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(78
+ ,284
+ ,9.100000381
+ ,109
+ ,8
+ ,68
+ ,433
+ ,8.699999809
+ ,144
+ ,9.300000191
+ ,70
+ ,739
+ ,7.199999809
+ ,113
+ ,7.5
+ ,96
+ ,1792
+ ,8.899999619
+ ,97
+ ,8.899999619
+ ,74
+ ,477
+ ,8.300000191
+ ,206
+ ,10.19999981
+ ,111
+ ,362
+ ,10.89999962
+ ,124
+ ,8.300000191
+ ,77
+ ,671
+ ,10
+ ,152
+ ,8.800000191
+ ,168
+ ,636
+ ,9.100000381
+ ,162
+ ,8.800000191
+ ,82
+ ,329
+ ,8.699999809
+ ,150
+ ,10.69999981
+ ,89
+ ,634
+ ,7.599999905
+ ,134
+ ,11.69999981
+ ,149
+ ,631
+ ,10.80000019
+ ,292
+ ,8.5
+ ,60
+ ,257
+ ,9.5
+ ,108
+ ,8.300000191
+ ,96
+ ,284
+ ,8.800000191
+ ,111
+ ,8.199999809
+ ,83
+ ,603
+ ,9.5
+ ,182
+ ,7.900000095
+ ,130
+ ,686
+ ,8.699999809
+ ,129
+ ,10.30000019
+ ,145
+ ,345
+ ,11.19999981
+ ,158
+ ,7.400000095
+ ,112
+ ,1357
+ ,9.699999809
+ ,186
+ ,9.600000381
+ ,131
+ ,544
+ ,9.600000381
+ ,177
+ ,9.300000191
+ ,80
+ ,205
+ ,9.100000381
+ ,127
+ ,10.60000038
+ ,130
+ ,1264
+ ,9.199999809
+ ,179
+ ,9.699999809
+ ,140
+ ,688
+ ,8.300000191
+ ,80
+ ,11.60000038
+ ,154
+ ,354
+ ,8.399999619
+ ,103
+ ,8.100000381
+ ,118
+ ,1632
+ ,9.399999619
+ ,101
+ ,9.800000191
+ ,94
+ ,348
+ ,9.800000191
+ ,117
+ ,7.400000095
+ ,119
+ ,370
+ ,10.39999962
+ ,88
+ ,9.399999619
+ ,153
+ ,648
+ ,9.899999619
+ ,78
+ ,11.19999981
+ ,116
+ ,366
+ ,9.199999809
+ ,102
+ ,9.100000381
+ ,97
+ ,540
+ ,10.30000019
+ ,95
+ ,10.5
+ ,176
+ ,680
+ ,8.899999619
+ ,80
+ ,11.89999962
+ ,75
+ ,345
+ ,9.600000381
+ ,92
+ ,8.399999619
+ ,134
+ ,525
+ ,10.30000019
+ ,126
+ ,5
+ ,161
+ ,870
+ ,10.39999962
+ ,108
+ ,9.800000191
+ ,111
+ ,669
+ ,9.699999809
+ ,77
+ ,9.800000191
+ ,114
+ ,452
+ ,9.600000381
+ ,60
+ ,10.80000019
+ ,142
+ ,430
+ ,10.69999981
+ ,71
+ ,10.10000038
+ ,238
+ ,822
+ ,10.30000019
+ ,86
+ ,10.89999962
+ ,78
+ ,190
+ ,10.69999981
+ ,93
+ ,9.199999809
+ ,196
+ ,867
+ ,9.600000381
+ ,106
+ ,8.300000191
+ ,125
+ ,969
+ ,10.5
+ ,162
+ ,7.300000191
+ ,82
+ ,499
+ ,7.699999809
+ ,95
+ ,9.399999619
+ ,125
+ ,925
+ ,10.19999981
+ ,91
+ ,9.399999619
+ ,129
+ ,353
+ ,9.899999619
+ ,52
+ ,9.800000191
+ ,84
+ ,288
+ ,8.399999619
+ ,110
+ ,3.599999905
+ ,183
+ ,718
+ ,10.39999962
+ ,69
+ ,8.399999619
+ ,119
+ ,540
+ ,9.199999809
+ ,57
+ ,10.80000019
+ ,180
+ ,668
+ ,13
+ ,106
+ ,10.10000038
+ ,82
+ ,347
+ ,8.800000191
+ ,40
+ ,9
+ ,71
+ ,345
+ ,9.199999809
+ ,50
+ ,10
+ ,118
+ ,463
+ ,7.800000191
+ ,35
+ ,11.30000019
+ ,121
+ ,728
+ ,8.199999809
+ ,86
+ ,11.30000019
+ ,68
+ ,383
+ ,7.400000095
+ ,57
+ ,12.80000019
+ ,112
+ ,316
+ ,10.39999962
+ ,57
+ ,10
+ ,109
+ ,388
+ ,8.899999619
+ ,94
+ ,6.699999809)
+ ,dim=c(5
+ ,53)
+ ,dimnames=list(c('Beschikbaarheiddokter'
+ ,'beschikbaarheidziekenhuis'
+ ,'inkomen'
+ ,'bevolkingsdichtheid'
+ ,'sterftecijfer')
+ ,1:53))
> y <- array(NA,dim=c(5,53),dimnames=list(c('Beschikbaarheiddokter','beschikbaarheidziekenhuis','inkomen','bevolkingsdichtheid','sterftecijfer'),1:53))
> 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 = '5'
> #'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
> 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
sterftecijfer Beschikbaarheiddokter beschikbaarheidziekenhuis inkomen
1 8.0 78 284 9.1
2 9.3 68 433 8.7
3 7.5 70 739 7.2
4 8.9 96 1792 8.9
5 10.2 74 477 8.3
6 8.3 111 362 10.9
7 8.8 77 671 10.0
8 8.8 168 636 9.1
9 10.7 82 329 8.7
10 11.7 89 634 7.6
11 8.5 149 631 10.8
12 8.3 60 257 9.5
13 8.2 96 284 8.8
14 7.9 83 603 9.5
15 10.3 130 686 8.7
16 7.4 145 345 11.2
17 9.6 112 1357 9.7
18 9.3 131 544 9.6
19 10.6 80 205 9.1
20 9.7 130 1264 9.2
21 11.6 140 688 8.3
22 8.1 154 354 8.4
23 9.8 118 1632 9.4
24 7.4 94 348 9.8
25 9.4 119 370 10.4
26 11.2 153 648 9.9
27 9.1 116 366 9.2
28 10.5 97 540 10.3
29 11.9 176 680 8.9
30 8.4 75 345 9.6
31 5.0 134 525 10.3
32 9.8 161 870 10.4
33 9.8 111 669 9.7
34 10.8 114 452 9.6
35 10.1 142 430 10.7
36 10.9 238 822 10.3
37 9.2 78 190 10.7
38 8.3 196 867 9.6
39 7.3 125 969 10.5
40 9.4 82 499 7.7
41 9.4 125 925 10.2
42 9.8 129 353 9.9
43 3.6 84 288 8.4
44 8.4 183 718 10.4
45 10.8 119 540 9.2
46 10.1 180 668 13.0
47 9.0 82 347 8.8
48 10.0 71 345 9.2
49 11.3 118 463 7.8
50 11.3 121 728 8.2
51 12.8 68 383 7.4
52 10.0 112 316 10.4
53 6.7 109 388 8.9
bevolkingsdichtheid
1 109
2 144
3 113
4 97
5 206
6 124
7 152
8 162
9 150
10 134
11 292
12 108
13 111
14 182
15 129
16 158
17 186
18 177
19 127
20 179
21 80
22 103
23 101
24 117
25 88
26 78
27 102
28 95
29 80
30 92
31 126
32 108
33 77
34 60
35 71
36 86
37 93
38 106
39 162
40 95
41 91
42 52
43 110
44 69
45 57
46 106
47 40
48 50
49 35
50 86
51 57
52 57
53 94
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Beschikbaarheiddokter
12.2662552 0.0073916
beschikbaarheidziekenhuis inkomen
0.0005837 -0.3302302
bevolkingsdichtheid
-0.0094629
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-5.6404 -0.7904 0.3053 0.9164 2.7906
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.2662552 2.0201467 6.072 1.95e-07 ***
Beschikbaarheiddokter 0.0073916 0.0069336 1.066 0.2917
beschikbaarheidziekenhuis 0.0005837 0.0007219 0.809 0.4228
inkomen -0.3302302 0.2345518 -1.408 0.1656
bevolkingsdichtheid -0.0094629 0.0048868 -1.936 0.0587 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.601 on 48 degrees of freedom
Multiple R-squared: 0.1437, Adjusted R-squared: 0.07235
F-statistic: 2.014 on 4 and 48 DF, p-value: 0.1075
> 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.04735538 0.09471077 0.9526446
[2,] 0.18431886 0.36863772 0.8156811
[3,] 0.45654232 0.91308464 0.5434577
[4,] 0.44680065 0.89360131 0.5531993
[5,] 0.33208288 0.66416577 0.6679171
[6,] 0.24800068 0.49600136 0.7519993
[7,] 0.20672543 0.41345085 0.7932746
[8,] 0.17879923 0.35759846 0.8212008
[9,] 0.12315396 0.24630792 0.8768460
[10,] 0.08844399 0.17688797 0.9115560
[11,] 0.07078761 0.14157521 0.9292124
[12,] 0.13920769 0.27841538 0.8607923
[13,] 0.13791671 0.27583342 0.8620833
[14,] 0.18165772 0.36331544 0.8183423
[15,] 0.19592614 0.39185227 0.8040739
[16,] 0.18234504 0.36469009 0.8176550
[17,] 0.15320961 0.30641921 0.8467904
[18,] 0.13469996 0.26939992 0.8653000
[19,] 0.16569258 0.33138517 0.8343074
[20,] 0.14393095 0.28786191 0.8560690
[21,] 0.16867180 0.33734359 0.8313282
[22,] 0.22046957 0.44093914 0.7795304
[23,] 0.16535669 0.33071337 0.8346433
[24,] 0.39131485 0.78262971 0.6086851
[25,] 0.31706157 0.63412315 0.6829384
[26,] 0.25087979 0.50175957 0.7491202
[27,] 0.21140037 0.42280074 0.7885996
[28,] 0.16518466 0.33036932 0.8348153
[29,] 0.15909670 0.31819340 0.8409033
[30,] 0.15886440 0.31772879 0.8411356
[31,] 0.13129618 0.26259236 0.8687038
[32,] 0.09541635 0.19083270 0.9045837
[33,] 0.07491765 0.14983529 0.9250824
[34,] 0.12493212 0.24986424 0.8750679
[35,] 0.09577038 0.19154075 0.9042296
[36,] 0.31268277 0.62536555 0.6873172
[37,] 0.33328742 0.66657484 0.6667126
[38,] 0.21176195 0.42352390 0.7882381
> postscript(file="/var/wessaorg/rcomp/tmp/1hyb81322067364.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/2j69e1322067364.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/3cmty1322067364.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/4l7xb1322067364.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/5yhry1322067364.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 = 53
Frequency = 1
1 2 3 4 5 6 7
-0.9720267 0.5140247 -2.2680705 -1.2649202 1.7985981 -0.2251221 0.3135783
8 9 10 11 12 13 14
-0.5412068 1.9280255 2.1835915 1.0937186 -0.4005880 -0.9852193 -0.4723074
15 16 17 18 19 20 21
0.7661213 -0.9457068 0.6771120 0.5930432 1.8296360 0.4669926 1.3952646
22 23 24 25 26 27 28
-1.7625878 -0.2311740 -1.4207860 0.3052957 1.4319640 -0.2339897 1.5018973
29 30 31 32 33 34 35
1.6319733 -0.5812128 -3.9694873 0.2922483 0.2546453 1.1652449 0.7384664
36 37 38 39 40 41 42
0.6099054 0.6598045 -1.7478168 -1.4554225 -0.3218953 -0.1006735 0.1355243
43 44 45 46 47 48 49
-5.6404098 -1.5506955 0.9164389 1.4093913 -0.7903754 0.5188206 0.7982710
50 51 52 53
1.2361106 2.7906403 0.6952086 -2.7698629
> postscript(file="/var/wessaorg/rcomp/tmp/6o4021322067364.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 = 53
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.9720267 NA
1 0.5140247 -0.9720267
2 -2.2680705 0.5140247
3 -1.2649202 -2.2680705
4 1.7985981 -1.2649202
5 -0.2251221 1.7985981
6 0.3135783 -0.2251221
7 -0.5412068 0.3135783
8 1.9280255 -0.5412068
9 2.1835915 1.9280255
10 1.0937186 2.1835915
11 -0.4005880 1.0937186
12 -0.9852193 -0.4005880
13 -0.4723074 -0.9852193
14 0.7661213 -0.4723074
15 -0.9457068 0.7661213
16 0.6771120 -0.9457068
17 0.5930432 0.6771120
18 1.8296360 0.5930432
19 0.4669926 1.8296360
20 1.3952646 0.4669926
21 -1.7625878 1.3952646
22 -0.2311740 -1.7625878
23 -1.4207860 -0.2311740
24 0.3052957 -1.4207860
25 1.4319640 0.3052957
26 -0.2339897 1.4319640
27 1.5018973 -0.2339897
28 1.6319733 1.5018973
29 -0.5812128 1.6319733
30 -3.9694873 -0.5812128
31 0.2922483 -3.9694873
32 0.2546453 0.2922483
33 1.1652449 0.2546453
34 0.7384664 1.1652449
35 0.6099054 0.7384664
36 0.6598045 0.6099054
37 -1.7478168 0.6598045
38 -1.4554225 -1.7478168
39 -0.3218953 -1.4554225
40 -0.1006735 -0.3218953
41 0.1355243 -0.1006735
42 -5.6404098 0.1355243
43 -1.5506955 -5.6404098
44 0.9164389 -1.5506955
45 1.4093913 0.9164389
46 -0.7903754 1.4093913
47 0.5188206 -0.7903754
48 0.7982710 0.5188206
49 1.2361106 0.7982710
50 2.7906403 1.2361106
51 0.6952086 2.7906403
52 -2.7698629 0.6952086
53 NA -2.7698629
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.5140247 -0.9720267
[2,] -2.2680705 0.5140247
[3,] -1.2649202 -2.2680705
[4,] 1.7985981 -1.2649202
[5,] -0.2251221 1.7985981
[6,] 0.3135783 -0.2251221
[7,] -0.5412068 0.3135783
[8,] 1.9280255 -0.5412068
[9,] 2.1835915 1.9280255
[10,] 1.0937186 2.1835915
[11,] -0.4005880 1.0937186
[12,] -0.9852193 -0.4005880
[13,] -0.4723074 -0.9852193
[14,] 0.7661213 -0.4723074
[15,] -0.9457068 0.7661213
[16,] 0.6771120 -0.9457068
[17,] 0.5930432 0.6771120
[18,] 1.8296360 0.5930432
[19,] 0.4669926 1.8296360
[20,] 1.3952646 0.4669926
[21,] -1.7625878 1.3952646
[22,] -0.2311740 -1.7625878
[23,] -1.4207860 -0.2311740
[24,] 0.3052957 -1.4207860
[25,] 1.4319640 0.3052957
[26,] -0.2339897 1.4319640
[27,] 1.5018973 -0.2339897
[28,] 1.6319733 1.5018973
[29,] -0.5812128 1.6319733
[30,] -3.9694873 -0.5812128
[31,] 0.2922483 -3.9694873
[32,] 0.2546453 0.2922483
[33,] 1.1652449 0.2546453
[34,] 0.7384664 1.1652449
[35,] 0.6099054 0.7384664
[36,] 0.6598045 0.6099054
[37,] -1.7478168 0.6598045
[38,] -1.4554225 -1.7478168
[39,] -0.3218953 -1.4554225
[40,] -0.1006735 -0.3218953
[41,] 0.1355243 -0.1006735
[42,] -5.6404098 0.1355243
[43,] -1.5506955 -5.6404098
[44,] 0.9164389 -1.5506955
[45,] 1.4093913 0.9164389
[46,] -0.7903754 1.4093913
[47,] 0.5188206 -0.7903754
[48,] 0.7982710 0.5188206
[49,] 1.2361106 0.7982710
[50,] 2.7906403 1.2361106
[51,] 0.6952086 2.7906403
[52,] -2.7698629 0.6952086
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.5140247 -0.9720267
2 -2.2680705 0.5140247
3 -1.2649202 -2.2680705
4 1.7985981 -1.2649202
5 -0.2251221 1.7985981
6 0.3135783 -0.2251221
7 -0.5412068 0.3135783
8 1.9280255 -0.5412068
9 2.1835915 1.9280255
10 1.0937186 2.1835915
11 -0.4005880 1.0937186
12 -0.9852193 -0.4005880
13 -0.4723074 -0.9852193
14 0.7661213 -0.4723074
15 -0.9457068 0.7661213
16 0.6771120 -0.9457068
17 0.5930432 0.6771120
18 1.8296360 0.5930432
19 0.4669926 1.8296360
20 1.3952646 0.4669926
21 -1.7625878 1.3952646
22 -0.2311740 -1.7625878
23 -1.4207860 -0.2311740
24 0.3052957 -1.4207860
25 1.4319640 0.3052957
26 -0.2339897 1.4319640
27 1.5018973 -0.2339897
28 1.6319733 1.5018973
29 -0.5812128 1.6319733
30 -3.9694873 -0.5812128
31 0.2922483 -3.9694873
32 0.2546453 0.2922483
33 1.1652449 0.2546453
34 0.7384664 1.1652449
35 0.6099054 0.7384664
36 0.6598045 0.6099054
37 -1.7478168 0.6598045
38 -1.4554225 -1.7478168
39 -0.3218953 -1.4554225
40 -0.1006735 -0.3218953
41 0.1355243 -0.1006735
42 -5.6404098 0.1355243
43 -1.5506955 -5.6404098
44 0.9164389 -1.5506955
45 1.4093913 0.9164389
46 -0.7903754 1.4093913
47 0.5188206 -0.7903754
48 0.7982710 0.5188206
49 1.2361106 0.7982710
50 2.7906403 1.2361106
51 0.6952086 2.7906403
52 -2.7698629 0.6952086
> 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/79wdd1322067364.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/8xi5f1322067364.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/9hhn11322067364.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/10nwg01322067364.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/11j9z51322067364.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/12hwgf1322067364.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/13d8e11322067364.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/14cf4m1322067364.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/15bhz91322067364.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/1605oi1322067364.tab")
+ }
>
> try(system("convert tmp/1hyb81322067364.ps tmp/1hyb81322067364.png",intern=TRUE))
character(0)
> try(system("convert tmp/2j69e1322067364.ps tmp/2j69e1322067364.png",intern=TRUE))
character(0)
> try(system("convert tmp/3cmty1322067364.ps tmp/3cmty1322067364.png",intern=TRUE))
character(0)
> try(system("convert tmp/4l7xb1322067364.ps tmp/4l7xb1322067364.png",intern=TRUE))
character(0)
> try(system("convert tmp/5yhry1322067364.ps tmp/5yhry1322067364.png",intern=TRUE))
character(0)
> try(system("convert tmp/6o4021322067364.ps tmp/6o4021322067364.png",intern=TRUE))
character(0)
> try(system("convert tmp/79wdd1322067364.ps tmp/79wdd1322067364.png",intern=TRUE))
character(0)
> try(system("convert tmp/8xi5f1322067364.ps tmp/8xi5f1322067364.png",intern=TRUE))
character(0)
> try(system("convert tmp/9hhn11322067364.ps tmp/9hhn11322067364.png",intern=TRUE))
character(0)
> try(system("convert tmp/10nwg01322067364.ps tmp/10nwg01322067364.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.170 0.474 3.695