R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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.
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(1962
+ ,9.5
+ ,5.569
+ ,1.933
+ ,0.226
+ ,1963
+ ,9.6
+ ,5.634
+ ,1.947
+ ,0.231
+ ,1964
+ ,9.4
+ ,5.433
+ ,1.936
+ ,0.225
+ ,1965
+ ,9.4
+ ,5.425
+ ,1.956
+ ,0.229
+ ,1966
+ ,9.5
+ ,5.412
+ ,1.965
+ ,0.236
+ ,1967
+ ,9.4
+ ,5.247
+ ,1.973
+ ,0.234
+ ,1968
+ ,9.7
+ ,5.31
+ ,1.988
+ ,0.253
+ ,1969
+ ,9.5
+ ,5.168
+ ,1.985
+ ,0.251
+ ,1970
+ ,9.5
+ ,4.927
+ ,1.986
+ ,0.243
+ ,1971
+ ,9.3
+ ,4.929
+ ,1.993
+ ,0.239
+ ,1972
+ ,9.4
+ ,4.902
+ ,2.003
+ ,0.237
+ ,1973
+ ,9.3
+ ,4.82
+ ,2
+ ,0.23
+ ,1974
+ ,9.1
+ ,4.588
+ ,2.015
+ ,0.221
+ ,1975
+ ,8.8
+ ,4.312
+ ,2.001
+ ,0.203
+ ,1976
+ ,8.8
+ ,4.269
+ ,2.025
+ ,0.195
+ ,1977
+ ,8.6
+ ,4.137
+ ,2.035
+ ,0.182
+ ,1978
+ ,8.7
+ ,4.099
+ ,2.049
+ ,0.183
+ ,1979
+ ,8.5
+ ,4.016
+ ,2.04
+ ,0.175
+ ,1980
+ ,8.7
+ ,4.121
+ ,2.079
+ ,0.181
+ ,1981
+ ,8.6
+ ,3.97
+ ,2.064
+ ,0.176
+ ,1982
+ ,8.5
+ ,3.89
+ ,2.083
+ ,0.172
+ ,1983
+ ,8.6
+ ,3.889
+ ,2.091
+ ,0.176
+ ,1984
+ ,8.6
+ ,3.788
+ ,2.108
+ ,0.172
+ ,1985
+ ,8.7
+ ,3.75
+ ,2.113
+ ,0.174
+ ,1986
+ ,8.7
+ ,3.651
+ ,2.115
+ ,0.172
+ ,1987
+ ,8.7
+ ,3.559
+ ,2.117
+ ,0.174
+ ,1988
+ ,8.8
+ ,3.525
+ ,2.125
+ ,0.18
+ ,1989
+ ,8.7
+ ,3.32
+ ,2.142
+ ,0.205
+ ,1990
+ ,8.6
+ ,3.218
+ ,2.16
+ ,0.207
+ ,1991
+ ,8.5
+ ,3.138
+ ,2.158
+ ,0.207
+ ,1992
+ ,8.5
+ ,3.061
+ ,2.143
+ ,0.208
+ ,1993
+ ,8.8
+ ,3.099
+ ,2.146
+ ,0.22
+ ,1994
+ ,8.8
+ ,2.997
+ ,2.131
+ ,0.227
+ ,1995
+ ,8.8
+ ,2.963
+ ,2.117
+ ,0.234
+ ,1996
+ ,8.8
+ ,2.883
+ ,2.087
+ ,0.24
+ ,1997
+ ,8.6
+ ,2.804
+ ,2.057
+ ,0.24
+ ,1998
+ ,8.6
+ ,2.724
+ ,2.024
+ ,0.242
+ ,1999
+ ,8.8
+ ,2.678
+ ,2.027
+ ,0.252
+ ,2000
+ ,8.7
+ ,2.576
+ ,1.996
+ ,0.25
+ ,2001
+ ,8.5
+ ,2.478
+ ,1.96
+ ,0.253)
+ ,dim=c(5
+ ,40)
+ ,dimnames=list(c('Year'
+ ,'Rate'
+ ,'Heart_disease'
+ ,'Cancer'
+ ,'Diabetes
')
+ ,1:40))
> y <- array(NA,dim=c(5,40),dimnames=list(c('Year','Rate','Heart_disease','Cancer','Diabetes
'),1:40))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par20 = ''
> par19 = ''
> par18 = ''
> par17 = ''
> par16 = ''
> par15 = ''
> par14 = ''
> par13 = ''
> par12 = ''
> par11 = ''
> par10 = ''
> par9 = ''
> par8 = ''
> par7 = ''
> par6 = ''
> par5 = ''
> par4 = ''
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '2'
> library(lattice)
> library(lmtest)
Loading required package: zoo
> 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
Rate Year Heart_disease Cancer Diabetes\r
1 9.5 1962 5.569 1.933 0.226
2 9.6 1963 5.634 1.947 0.231
3 9.4 1964 5.433 1.936 0.225
4 9.4 1965 5.425 1.956 0.229
5 9.5 1966 5.412 1.965 0.236
6 9.4 1967 5.247 1.973 0.234
7 9.7 1968 5.310 1.988 0.253
8 9.5 1969 5.168 1.985 0.251
9 9.5 1970 4.927 1.986 0.243
10 9.3 1971 4.929 1.993 0.239
11 9.4 1972 4.902 2.003 0.237
12 9.3 1973 4.820 2.000 0.230
13 9.1 1974 4.588 2.015 0.221
14 8.8 1975 4.312 2.001 0.203
15 8.8 1976 4.269 2.025 0.195
16 8.6 1977 4.137 2.035 0.182
17 8.7 1978 4.099 2.049 0.183
18 8.5 1979 4.016 2.040 0.175
19 8.7 1980 4.121 2.079 0.181
20 8.6 1981 3.970 2.064 0.176
21 8.5 1982 3.890 2.083 0.172
22 8.6 1983 3.889 2.091 0.176
23 8.6 1984 3.788 2.108 0.172
24 8.7 1985 3.750 2.113 0.174
25 8.7 1986 3.651 2.115 0.172
26 8.7 1987 3.559 2.117 0.174
27 8.8 1988 3.525 2.125 0.180
28 8.7 1989 3.320 2.142 0.205
29 8.6 1990 3.218 2.160 0.207
30 8.5 1991 3.138 2.158 0.207
31 8.5 1992 3.061 2.143 0.208
32 8.8 1993 3.099 2.146 0.220
33 8.8 1994 2.997 2.131 0.227
34 8.8 1995 2.963 2.117 0.234
35 8.8 1996 2.883 2.087 0.240
36 8.6 1997 2.804 2.057 0.240
37 8.6 1998 2.724 2.024 0.242
38 8.8 1999 2.678 2.027 0.252
39 8.7 2000 2.576 1.996 0.250
40 8.5 2001 2.478 1.960 0.253
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Year Heart_disease Cancer `Diabetes\r`
-101.63575 0.05212 0.98985 0.98565 6.05031
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.117957 -0.054671 -0.003596 0.048395 0.157415
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -101.63575 25.24499 -4.026 0.000290 ***
Year 0.05212 0.01256 4.151 0.000201 ***
Heart_disease 0.98985 0.15026 6.588 1.31e-07 ***
Cancer 0.98565 0.31602 3.119 0.003622 **
`Diabetes\r` 6.05031 0.70730 8.554 4.27e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.07624 on 35 degrees of freedom
Multiple R-squared: 0.9644, Adjusted R-squared: 0.9603
F-statistic: 236.9 on 4 and 35 DF, p-value: < 2.2e-16
> 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
+ }
[,1] [,2] [,3]
[1,] 0.1192790 0.2385581 0.8807210
[2,] 0.4331268 0.8662536 0.5668732
[3,] 0.2974580 0.5949161 0.7025420
[4,] 0.4223472 0.8446945 0.5776528
[5,] 0.3953979 0.7907959 0.6046021
[6,] 0.3075487 0.6150973 0.6924513
[7,] 0.2208356 0.4416711 0.7791644
[8,] 0.1700438 0.3400875 0.8299562
[9,] 0.1269822 0.2539645 0.8730178
[10,] 0.2276232 0.4552465 0.7723768
[11,] 0.3129933 0.6259865 0.6870067
[12,] 0.2273131 0.4546262 0.7726869
[13,] 0.2923051 0.5846103 0.7076949
[14,] 0.2351293 0.4702585 0.7648707
[15,] 0.2016941 0.4033883 0.7983059
[16,] 0.2070536 0.4141071 0.7929464
[17,] 0.4573920 0.9147841 0.5426080
[18,] 0.6237182 0.7525636 0.3762818
[19,] 0.6618441 0.6763119 0.3381559
[20,] 0.6808063 0.6383875 0.3191937
[21,] 0.6882918 0.6234165 0.3117082
[22,] 0.8300602 0.3398795 0.1699398
[23,] 0.7730667 0.4538667 0.2269333
[24,] 0.6490020 0.7019961 0.3509980
[25,] 0.5644639 0.8710723 0.4355361
> postscript(file="/var/wessaorg/rcomp/tmp/1qwr41321954596.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/wessaorg/rcomp/tmp/2z7841321954596.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/wessaorg/rcomp/tmp/3citn1321954596.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/wessaorg/rcomp/tmp/4pj9f1321954596.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/wessaorg/rcomp/tmp/5dg0q1321954596.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 = 40
Frequency = 1
1 2 3 4 5 6
0.093976320 0.033466892 0.027451887 -0.060662206 -0.051135804 -0.035713994
7 8 9 10 11 12
0.020066296 -0.076436230 0.157415491 -0.079381130 -0.002529720 -0.028171611
13 14 15 16 17 18
-0.010977356 0.032806944 0.047998647 -0.004662449 0.060983757 -0.051704061
19 20 21 22 23 24
-0.082499055 -0.040114185 -0.107571060 -0.090786260 -0.035485011 0.032981760
25 26 27 28 29 30
0.088987485 0.113863034 0.151212216 0.033998748 -0.046997662 -0.117957066
31 32 33 34 35 36
-0.085122840 0.049583659 0.070862243 0.023845476 0.044182493 -0.100168484
37 38 39 40
-0.052673273 0.077281119 0.068782928 -0.068997938
> postscript(file="/var/wessaorg/rcomp/tmp/6l7301321954596.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 = 40
Frequency = 1
lag(myerror, k = 1) myerror
0 0.093976320 NA
1 0.033466892 0.093976320
2 0.027451887 0.033466892
3 -0.060662206 0.027451887
4 -0.051135804 -0.060662206
5 -0.035713994 -0.051135804
6 0.020066296 -0.035713994
7 -0.076436230 0.020066296
8 0.157415491 -0.076436230
9 -0.079381130 0.157415491
10 -0.002529720 -0.079381130
11 -0.028171611 -0.002529720
12 -0.010977356 -0.028171611
13 0.032806944 -0.010977356
14 0.047998647 0.032806944
15 -0.004662449 0.047998647
16 0.060983757 -0.004662449
17 -0.051704061 0.060983757
18 -0.082499055 -0.051704061
19 -0.040114185 -0.082499055
20 -0.107571060 -0.040114185
21 -0.090786260 -0.107571060
22 -0.035485011 -0.090786260
23 0.032981760 -0.035485011
24 0.088987485 0.032981760
25 0.113863034 0.088987485
26 0.151212216 0.113863034
27 0.033998748 0.151212216
28 -0.046997662 0.033998748
29 -0.117957066 -0.046997662
30 -0.085122840 -0.117957066
31 0.049583659 -0.085122840
32 0.070862243 0.049583659
33 0.023845476 0.070862243
34 0.044182493 0.023845476
35 -0.100168484 0.044182493
36 -0.052673273 -0.100168484
37 0.077281119 -0.052673273
38 0.068782928 0.077281119
39 -0.068997938 0.068782928
40 NA -0.068997938
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.033466892 0.093976320
[2,] 0.027451887 0.033466892
[3,] -0.060662206 0.027451887
[4,] -0.051135804 -0.060662206
[5,] -0.035713994 -0.051135804
[6,] 0.020066296 -0.035713994
[7,] -0.076436230 0.020066296
[8,] 0.157415491 -0.076436230
[9,] -0.079381130 0.157415491
[10,] -0.002529720 -0.079381130
[11,] -0.028171611 -0.002529720
[12,] -0.010977356 -0.028171611
[13,] 0.032806944 -0.010977356
[14,] 0.047998647 0.032806944
[15,] -0.004662449 0.047998647
[16,] 0.060983757 -0.004662449
[17,] -0.051704061 0.060983757
[18,] -0.082499055 -0.051704061
[19,] -0.040114185 -0.082499055
[20,] -0.107571060 -0.040114185
[21,] -0.090786260 -0.107571060
[22,] -0.035485011 -0.090786260
[23,] 0.032981760 -0.035485011
[24,] 0.088987485 0.032981760
[25,] 0.113863034 0.088987485
[26,] 0.151212216 0.113863034
[27,] 0.033998748 0.151212216
[28,] -0.046997662 0.033998748
[29,] -0.117957066 -0.046997662
[30,] -0.085122840 -0.117957066
[31,] 0.049583659 -0.085122840
[32,] 0.070862243 0.049583659
[33,] 0.023845476 0.070862243
[34,] 0.044182493 0.023845476
[35,] -0.100168484 0.044182493
[36,] -0.052673273 -0.100168484
[37,] 0.077281119 -0.052673273
[38,] 0.068782928 0.077281119
[39,] -0.068997938 0.068782928
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.033466892 0.093976320
2 0.027451887 0.033466892
3 -0.060662206 0.027451887
4 -0.051135804 -0.060662206
5 -0.035713994 -0.051135804
6 0.020066296 -0.035713994
7 -0.076436230 0.020066296
8 0.157415491 -0.076436230
9 -0.079381130 0.157415491
10 -0.002529720 -0.079381130
11 -0.028171611 -0.002529720
12 -0.010977356 -0.028171611
13 0.032806944 -0.010977356
14 0.047998647 0.032806944
15 -0.004662449 0.047998647
16 0.060983757 -0.004662449
17 -0.051704061 0.060983757
18 -0.082499055 -0.051704061
19 -0.040114185 -0.082499055
20 -0.107571060 -0.040114185
21 -0.090786260 -0.107571060
22 -0.035485011 -0.090786260
23 0.032981760 -0.035485011
24 0.088987485 0.032981760
25 0.113863034 0.088987485
26 0.151212216 0.113863034
27 0.033998748 0.151212216
28 -0.046997662 0.033998748
29 -0.117957066 -0.046997662
30 -0.085122840 -0.117957066
31 0.049583659 -0.085122840
32 0.070862243 0.049583659
33 0.023845476 0.070862243
34 0.044182493 0.023845476
35 -0.100168484 0.044182493
36 -0.052673273 -0.100168484
37 0.077281119 -0.052673273
38 0.068782928 0.077281119
39 -0.068997938 0.068782928
> 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/wessaorg/rcomp/tmp/768c41321954596.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/wessaorg/rcomp/tmp/8m6mc1321954596.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/wessaorg/rcomp/tmp/93eh31321954596.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/wessaorg/rcomp/tmp/10lp7x1321954596.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()
+ }
null device
1
>
> #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/1180lx1321954596.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/wessaorg/rcomp/tmp/12judw1321954596.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/wessaorg/rcomp/tmp/13aemv1321954596.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/wessaorg/rcomp/tmp/14lqwr1321954596.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/wessaorg/rcomp/tmp/15y9vk1321954596.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/wessaorg/rcomp/tmp/16gixj1321954596.tab")
+ }
>
> try(system("convert tmp/1qwr41321954596.ps tmp/1qwr41321954596.png",intern=TRUE))
character(0)
> try(system("convert tmp/2z7841321954596.ps tmp/2z7841321954596.png",intern=TRUE))
character(0)
> try(system("convert tmp/3citn1321954596.ps tmp/3citn1321954596.png",intern=TRUE))
character(0)
> try(system("convert tmp/4pj9f1321954596.ps tmp/4pj9f1321954596.png",intern=TRUE))
character(0)
> try(system("convert tmp/5dg0q1321954596.ps tmp/5dg0q1321954596.png",intern=TRUE))
character(0)
> try(system("convert tmp/6l7301321954596.ps tmp/6l7301321954596.png",intern=TRUE))
character(0)
> try(system("convert tmp/768c41321954596.ps tmp/768c41321954596.png",intern=TRUE))
character(0)
> try(system("convert tmp/8m6mc1321954596.ps tmp/8m6mc1321954596.png",intern=TRUE))
character(0)
> try(system("convert tmp/93eh31321954596.ps tmp/93eh31321954596.png",intern=TRUE))
character(0)
> try(system("convert tmp/10lp7x1321954596.ps tmp/10lp7x1321954596.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.055 0.739 3.814