R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(2148,77.405,82.145,315.4,2118,85.056,78.213,329.3,1603,90.088,88.099,308.2,2066,99.285,106.25,335.8,2095,80.428,80.487,343.7,2210,88.017,80.336,349.2,1609,93.489,90.065,312.4,1964,103.961,108.888,337.6,2114,82.591,82.747,360.2,2054,90.913,82.213,372.1,1424,96.787,93.41,341.8,2025,106.045,109.465,377.4,2003,84.752,84.373,337.2,2017,94.173,98.715,384.6,1528,97.733,99.646,358.6,2130,108.499,115.239,383.4,2017,87.972,89.082,384.4,2260,96.091,89.934,402.7,1805,101.846,99.957,372.1,2394,115.652,122.717,364.9,2586,91.269,95.895,314.9,2429,100.911,97.085,320.7,1910,105.248,109.414,308.6,2515,118.681,126.945,328.7),dim=c(4,24),dimnames=list(c('FallBelg','Brutoloonindex','wgbijdrage','brutoomzetindex
'),1:24))
> y <- array(NA,dim=c(4,24),dimnames=list(c('FallBelg','Brutoloonindex','wgbijdrage','brutoomzetindex
'),1:24))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
FallBelg Brutoloonindex wgbijdrage brutoomzetindex\r
1 2148 77.405 82.145 315.4
2 2118 85.056 78.213 329.3
3 1603 90.088 88.099 308.2
4 2066 99.285 106.250 335.8
5 2095 80.428 80.487 343.7
6 2210 88.017 80.336 349.2
7 1609 93.489 90.065 312.4
8 1964 103.961 108.888 337.6
9 2114 82.591 82.747 360.2
10 2054 90.913 82.213 372.1
11 1424 96.787 93.410 341.8
12 2025 106.045 109.465 377.4
13 2003 84.752 84.373 337.2
14 2017 94.173 98.715 384.6
15 1528 97.733 99.646 358.6
16 2130 108.499 115.239 383.4
17 2017 87.972 89.082 384.4
18 2260 96.091 89.934 402.7
19 1805 101.846 99.957 372.1
20 2394 115.652 122.717 364.9
21 2586 91.269 95.895 314.9
22 2429 100.911 97.085 320.7
23 1910 105.248 109.414 308.6
24 2515 118.681 126.945 328.7
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Brutoloonindex wgbijdrage
1702.1107 -11.3569 12.7424
`brutoomzetindex\r`
0.5761
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-566.10 -146.78 -17.59 204.49 517.07
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1702.1107 944.4516 1.802 0.0866 .
Brutoloonindex -11.3569 19.0160 -0.597 0.5571
wgbijdrage 12.7424 14.5645 0.875 0.3920
`brutoomzetindex\r` 0.5761 2.3458 0.246 0.8085
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 307 on 20 degrees of freedom
Multiple R-squared: 0.0628, Adjusted R-squared: -0.07778
F-statistic: 0.4467 on 3 and 20 DF, p-value: 0.7223
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
> postscript(file="/var/www/html/freestat/rcomp/tmp/1315j1291976423.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/2315j1291976423.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/3vanm1291976423.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/4vanm1291976423.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/5vanm1291976423.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 24
Frequency = 1
1 2 3 4 5 6
96.535638 195.522369 -376.145888 -55.885400 82.690622 282.633340
7 8 9 10 11 12
-358.992496 -139.432285 87.951862 122.412441 -566.098520 -85.045946
13 14 15 16 17 18
-5.974774 -95.041265 -540.495349 -29.207698 -42.601983 271.205062
19 20 21 22 23 24
-228.524909 231.397862 517.066594 451.064494 -168.811501 353.777728
> postscript(file="/var/www/html/freestat/rcomp/tmp/6o2m71291976423.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 24
Frequency = 1
lag(myerror, k = 1) myerror
0 96.535638 NA
1 195.522369 96.535638
2 -376.145888 195.522369
3 -55.885400 -376.145888
4 82.690622 -55.885400
5 282.633340 82.690622
6 -358.992496 282.633340
7 -139.432285 -358.992496
8 87.951862 -139.432285
9 122.412441 87.951862
10 -566.098520 122.412441
11 -85.045946 -566.098520
12 -5.974774 -85.045946
13 -95.041265 -5.974774
14 -540.495349 -95.041265
15 -29.207698 -540.495349
16 -42.601983 -29.207698
17 271.205062 -42.601983
18 -228.524909 271.205062
19 231.397862 -228.524909
20 517.066594 231.397862
21 451.064494 517.066594
22 -168.811501 451.064494
23 353.777728 -168.811501
24 NA 353.777728
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 195.522369 96.535638
[2,] -376.145888 195.522369
[3,] -55.885400 -376.145888
[4,] 82.690622 -55.885400
[5,] 282.633340 82.690622
[6,] -358.992496 282.633340
[7,] -139.432285 -358.992496
[8,] 87.951862 -139.432285
[9,] 122.412441 87.951862
[10,] -566.098520 122.412441
[11,] -85.045946 -566.098520
[12,] -5.974774 -85.045946
[13,] -95.041265 -5.974774
[14,] -540.495349 -95.041265
[15,] -29.207698 -540.495349
[16,] -42.601983 -29.207698
[17,] 271.205062 -42.601983
[18,] -228.524909 271.205062
[19,] 231.397862 -228.524909
[20,] 517.066594 231.397862
[21,] 451.064494 517.066594
[22,] -168.811501 451.064494
[23,] 353.777728 -168.811501
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 195.522369 96.535638
2 -376.145888 195.522369
3 -55.885400 -376.145888
4 82.690622 -55.885400
5 282.633340 82.690622
6 -358.992496 282.633340
7 -139.432285 -358.992496
8 87.951862 -139.432285
9 122.412441 87.951862
10 -566.098520 122.412441
11 -85.045946 -566.098520
12 -5.974774 -85.045946
13 -95.041265 -5.974774
14 -540.495349 -95.041265
15 -29.207698 -540.495349
16 -42.601983 -29.207698
17 271.205062 -42.601983
18 -228.524909 271.205062
19 231.397862 -228.524909
20 517.066594 231.397862
21 451.064494 517.066594
22 -168.811501 451.064494
23 353.777728 -168.811501
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/7zt3a1291976423.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/8zt3a1291976423.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/9zt3a1291976423.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10s22d1291976423.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
>
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/11v3j11291976423.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/12h3071291976423.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/13cdfy1291976423.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/14ydem1291976423.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/151wur1291976423.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/165xtx1291976423.tab")
+ }
>
> try(system("convert tmp/1315j1291976423.ps tmp/1315j1291976423.png",intern=TRUE))
character(0)
> try(system("convert tmp/2315j1291976423.ps tmp/2315j1291976423.png",intern=TRUE))
character(0)
> try(system("convert tmp/3vanm1291976423.ps tmp/3vanm1291976423.png",intern=TRUE))
character(0)
> try(system("convert tmp/4vanm1291976423.ps tmp/4vanm1291976423.png",intern=TRUE))
character(0)
> try(system("convert tmp/5vanm1291976423.ps tmp/5vanm1291976423.png",intern=TRUE))
character(0)
> try(system("convert tmp/6o2m71291976423.ps tmp/6o2m71291976423.png",intern=TRUE))
character(0)
> try(system("convert tmp/7zt3a1291976423.ps tmp/7zt3a1291976423.png",intern=TRUE))
character(0)
> try(system("convert tmp/8zt3a1291976423.ps tmp/8zt3a1291976423.png",intern=TRUE))
character(0)
> try(system("convert tmp/9zt3a1291976423.ps tmp/9zt3a1291976423.png",intern=TRUE))
character(0)
> try(system("convert tmp/10s22d1291976423.ps tmp/10s22d1291976423.png",intern=TRUE))
convert: unable to open image `tmp/10s22d1291976423.ps': No such file or directory.
convert: missing an image filename `tmp/10s22d1291976423.png'.
character(0)
>
>
> proc.time()
user system elapsed
3.159 2.211 3.478