R version 3.2.2 (2015-08-14) -- "Fire Safety"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-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(10
+ ,70.3
+ ,213
+ ,582
+ ,6
+ ,7.05
+ ,36
+ ,13
+ ,61
+ ,91
+ ,132
+ ,8.2
+ ,48.52
+ ,100
+ ,12
+ ,56.7
+ ,453
+ ,716
+ ,8.7
+ ,20.66
+ ,67
+ ,17
+ ,51.9
+ ,454
+ ,515
+ ,9
+ ,12.95
+ ,86
+ ,56
+ ,49.1
+ ,412
+ ,158
+ ,9
+ ,43.37
+ ,127
+ ,36
+ ,54
+ ,80
+ ,80
+ ,9
+ ,40.25
+ ,114
+ ,29
+ ,57.3
+ ,434
+ ,757
+ ,9.3
+ ,38.89
+ ,111
+ ,14
+ ,68.4
+ ,136
+ ,529
+ ,8.8
+ ,54.47
+ ,116
+ ,10
+ ,75.5
+ ,207
+ ,335
+ ,9
+ ,59.8
+ ,128
+ ,24
+ ,61.5
+ ,368
+ ,497
+ ,9.1
+ ,48.34
+ ,115
+ ,110
+ ,50.6
+ ,3344
+ ,3369
+ ,10.4
+ ,34.44
+ ,122
+ ,28
+ ,52.3
+ ,361
+ ,746
+ ,9.7
+ ,38.74
+ ,121
+ ,17
+ ,49
+ ,104
+ ,201
+ ,11.2
+ ,30.85
+ ,103
+ ,8
+ ,56.6
+ ,125
+ ,277
+ ,12.7
+ ,30.58
+ ,82
+ ,30
+ ,55.6
+ ,291
+ ,593
+ ,8.3
+ ,43.11
+ ,123
+ ,9
+ ,68.3
+ ,204
+ ,361
+ ,8.4
+ ,56.77
+ ,113
+ ,47
+ ,55
+ ,625
+ ,905
+ ,9.6
+ ,41.31
+ ,111
+ ,35
+ ,49.9
+ ,1064
+ ,1513
+ ,10.1
+ ,30.96
+ ,129
+ ,29
+ ,43.5
+ ,699
+ ,744
+ ,10.6
+ ,25.94
+ ,137
+ ,14
+ ,54.5
+ ,381
+ ,507
+ ,10
+ ,37
+ ,99
+ ,56
+ ,55.9
+ ,775
+ ,622
+ ,9.5
+ ,35.89
+ ,105
+ ,14
+ ,51.5
+ ,181
+ ,347
+ ,10.9
+ ,30.18
+ ,98
+ ,11
+ ,56.8
+ ,46
+ ,244
+ ,8.9
+ ,7.77
+ ,58
+ ,46
+ ,47.6
+ ,44
+ ,116
+ ,8.8
+ ,33.36
+ ,135
+ ,11
+ ,47.1
+ ,391
+ ,463
+ ,12.4
+ ,36.11
+ ,166
+ ,23
+ ,54
+ ,462
+ ,453
+ ,7.1
+ ,39.04
+ ,132
+ ,65
+ ,49.7
+ ,1007
+ ,751
+ ,10.9
+ ,34.99
+ ,155
+ ,26
+ ,51.5
+ ,266
+ ,540
+ ,8.6
+ ,37.01
+ ,134
+ ,69
+ ,54.6
+ ,1692
+ ,1950
+ ,9.6
+ ,39.93
+ ,115
+ ,61
+ ,50.4
+ ,347
+ ,520
+ ,9.4
+ ,36.22
+ ,147
+ ,94
+ ,50
+ ,343
+ ,179
+ ,10.6
+ ,42.75
+ ,125
+ ,10
+ ,61.6
+ ,337
+ ,624
+ ,9.2
+ ,49.1
+ ,105
+ ,18
+ ,59.4
+ ,275
+ ,448
+ ,7.9
+ ,46
+ ,119
+ ,9
+ ,66.2
+ ,641
+ ,844
+ ,10.9
+ ,35.94
+ ,78
+ ,10
+ ,68.9
+ ,721
+ ,1233
+ ,10.8
+ ,48.19
+ ,103
+ ,28
+ ,51
+ ,137
+ ,176
+ ,8.7
+ ,15.17
+ ,89
+ ,31
+ ,59.3
+ ,96
+ ,308
+ ,10.6
+ ,44.68
+ ,116
+ ,26
+ ,57.8
+ ,197
+ ,299
+ ,7.6
+ ,42.59
+ ,115
+ ,29
+ ,51.1
+ ,379
+ ,531
+ ,9.4
+ ,38.79
+ ,164
+ ,31
+ ,55.2
+ ,35
+ ,71
+ ,6.5
+ ,40.75
+ ,148
+ ,16
+ ,45.7
+ ,569
+ ,717
+ ,11.8
+ ,29.07
+ ,123)
+ ,dim=c(7
+ ,41)
+ ,dimnames=list(c('SO2'
+ ,'Temp'
+ ,'Man'
+ ,'Pop'
+ ,'Wind'
+ ,'Rain'
+ ,'RainDays')
+ ,1:41))
> y <- array(NA,dim=c(7,41),dimnames=list(c('SO2','Temp','Man','Pop','Wind','Rain','RainDays'),1:41))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par5 = '0'
> par4 = '0'
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> par5 <- '0'
> par4 <- '0'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal Dummies'
> par1 <- '1'
> #'GNU S' R Code compiled by R2WASP v. 1.2.327 (Mon, 30 Nov 2015 20:34:54 +0000)
> #Author: root
> #To cite this work: Wessa P., (2015), Multiple Regression (v1.0.36) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_multipleregression.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> mywarning <- ''
> par1 <- as.numeric(par1)
> if(is.na(par1)) {
+ par1 <- 1
+ mywarning = 'Warning: you did not specify the column number of the endogenous series! The first column was selected by default.'
+ }
> if (par4=='') par4 <- 0
> par4 <- as.numeric(par4)
> if (par5=='') par5 <- 0
> par5 <- as.numeric(par5)
> x <- na.omit(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(par4 > 0) {
+ x2 <- array(0, dim=c(n-par4,par4), dimnames=list(1:(n-par4), paste(colnames(x)[par1],'(t-',1:par4,')',sep='')))
+ for (i in 1:(n-par4)) {
+ for (j in 1:par4) {
+ x2[i,j] <- x[i+par4-j,par1]
+ }
+ }
+ x <- cbind(x[(par4+1):n,], x2)
+ n <- n - par4
+ }
> 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+par4,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
SO2 Temp Man Pop Wind Rain RainDays
1 10 70.3 213 582 6.0 7.05 36
2 13 61.0 91 132 8.2 48.52 100
3 12 56.7 453 716 8.7 20.66 67
4 17 51.9 454 515 9.0 12.95 86
5 56 49.1 412 158 9.0 43.37 127
6 36 54.0 80 80 9.0 40.25 114
7 29 57.3 434 757 9.3 38.89 111
8 14 68.4 136 529 8.8 54.47 116
9 10 75.5 207 335 9.0 59.80 128
10 24 61.5 368 497 9.1 48.34 115
11 110 50.6 3344 3369 10.4 34.44 122
12 28 52.3 361 746 9.7 38.74 121
13 17 49.0 104 201 11.2 30.85 103
14 8 56.6 125 277 12.7 30.58 82
15 30 55.6 291 593 8.3 43.11 123
16 9 68.3 204 361 8.4 56.77 113
17 47 55.0 625 905 9.6 41.31 111
18 35 49.9 1064 1513 10.1 30.96 129
19 29 43.5 699 744 10.6 25.94 137
20 14 54.5 381 507 10.0 37.00 99
21 56 55.9 775 622 9.5 35.89 105
22 14 51.5 181 347 10.9 30.18 98
23 11 56.8 46 244 8.9 7.77 58
24 46 47.6 44 116 8.8 33.36 135
25 11 47.1 391 463 12.4 36.11 166
26 23 54.0 462 453 7.1 39.04 132
27 65 49.7 1007 751 10.9 34.99 155
28 26 51.5 266 540 8.6 37.01 134
29 69 54.6 1692 1950 9.6 39.93 115
30 61 50.4 347 520 9.4 36.22 147
31 94 50.0 343 179 10.6 42.75 125
32 10 61.6 337 624 9.2 49.10 105
33 18 59.4 275 448 7.9 46.00 119
34 9 66.2 641 844 10.9 35.94 78
35 10 68.9 721 1233 10.8 48.19 103
36 28 51.0 137 176 8.7 15.17 89
37 31 59.3 96 308 10.6 44.68 116
38 26 57.8 197 299 7.6 42.59 115
39 29 51.1 379 531 9.4 38.79 164
40 31 55.2 35 71 6.5 40.75 148
41 16 45.7 569 717 11.8 29.07 123
> k <- length(x[1+par4,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Temp Man Pop Wind Rain
111.72848 -1.26794 0.06492 -0.03928 -3.18137 0.51236
RainDays
-0.05205
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-23.004 -8.542 -0.991 5.758 48.758
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 111.72848 47.31810 2.361 0.024087 *
Temp -1.26794 0.62118 -2.041 0.049056 *
Man 0.06492 0.01575 4.122 0.000228 ***
Pop -0.03928 0.01513 -2.595 0.013846 *
Wind -3.18137 1.81502 -1.753 0.088650 .
Rain 0.51236 0.36276 1.412 0.166918
RainDays -0.05205 0.16201 -0.321 0.749972
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.64 on 34 degrees of freedom
Multiple R-squared: 0.6695, Adjusted R-squared: 0.6112
F-statistic: 11.48 on 6 and 34 DF, p-value: 5.419e-07
> 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.140510597 0.281021193 0.85948940
[2,] 0.054423745 0.108847489 0.94557626
[3,] 0.026710397 0.053420794 0.97328960
[4,] 0.040572149 0.081144297 0.95942785
[5,] 0.039936179 0.079872358 0.96006382
[6,] 0.017667321 0.035334643 0.98233268
[7,] 0.012286407 0.024572813 0.98771359
[8,] 0.014750252 0.029500505 0.98524975
[9,] 0.011351848 0.022703696 0.98864815
[10,] 0.010237593 0.020475185 0.98976241
[11,] 0.012054577 0.024109154 0.98794542
[12,] 0.011284878 0.022569755 0.98871512
[13,] 0.006570141 0.013140282 0.99342986
[14,] 0.003711478 0.007422956 0.99628852
[15,] 0.003905165 0.007810330 0.99609484
[16,] 0.006653024 0.013306047 0.99334698
[17,] 0.017054233 0.034108466 0.98294577
[18,] 0.028104259 0.056208519 0.97189574
[19,] 0.015435843 0.030871686 0.98456416
[20,] 0.014825521 0.029651042 0.98517448
[21,] 0.178713129 0.357426259 0.82128687
[22,] 0.973384261 0.053231478 0.02661574
> postscript(file="/var/wessaorg/rcomp/tmp/12hko1449264906.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/20alb1449264906.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/33r231449264906.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/41c891449264906.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/545cx1449264906.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 = 41
Frequency = 1
1 2 3 4 5 6
13.7891430 -15.6745359 -8.5420949 -11.6941046 -0.9914755 4.6325896
7 8 9 10 11 12
6.9211854 7.0728641 -1.6236301 -3.9506805 -0.5429891 5.7583641
13 14 15 16 17 18
-6.2700392 -0.1946044 3.8886796 -11.6739808 15.1330095 -0.1218054
19 20 21 22 23 24
-16.1661852 -15.4368601 6.5676608 -6.2358612 5.2399001 14.2557175
25 26 27 28 29 30
-18.6183425 -23.0036633 5.2296552 -1.1162405 9.3175985 30.0716264
31 32 33 34 35 36
48.7575842 -11.4150350 -10.9110457 -6.9311354 2.2842257 3.4932963
37 38 39 40 41
17.1936194 -5.1438353 -3.1180614 -2.7179416 -17.5125721
> postscript(file="/var/wessaorg/rcomp/tmp/61kcq1449264906.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 = 41
Frequency = 1
lag(myerror, k = 1) myerror
0 13.7891430 NA
1 -15.6745359 13.7891430
2 -8.5420949 -15.6745359
3 -11.6941046 -8.5420949
4 -0.9914755 -11.6941046
5 4.6325896 -0.9914755
6 6.9211854 4.6325896
7 7.0728641 6.9211854
8 -1.6236301 7.0728641
9 -3.9506805 -1.6236301
10 -0.5429891 -3.9506805
11 5.7583641 -0.5429891
12 -6.2700392 5.7583641
13 -0.1946044 -6.2700392
14 3.8886796 -0.1946044
15 -11.6739808 3.8886796
16 15.1330095 -11.6739808
17 -0.1218054 15.1330095
18 -16.1661852 -0.1218054
19 -15.4368601 -16.1661852
20 6.5676608 -15.4368601
21 -6.2358612 6.5676608
22 5.2399001 -6.2358612
23 14.2557175 5.2399001
24 -18.6183425 14.2557175
25 -23.0036633 -18.6183425
26 5.2296552 -23.0036633
27 -1.1162405 5.2296552
28 9.3175985 -1.1162405
29 30.0716264 9.3175985
30 48.7575842 30.0716264
31 -11.4150350 48.7575842
32 -10.9110457 -11.4150350
33 -6.9311354 -10.9110457
34 2.2842257 -6.9311354
35 3.4932963 2.2842257
36 17.1936194 3.4932963
37 -5.1438353 17.1936194
38 -3.1180614 -5.1438353
39 -2.7179416 -3.1180614
40 -17.5125721 -2.7179416
41 NA -17.5125721
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -15.6745359 13.7891430
[2,] -8.5420949 -15.6745359
[3,] -11.6941046 -8.5420949
[4,] -0.9914755 -11.6941046
[5,] 4.6325896 -0.9914755
[6,] 6.9211854 4.6325896
[7,] 7.0728641 6.9211854
[8,] -1.6236301 7.0728641
[9,] -3.9506805 -1.6236301
[10,] -0.5429891 -3.9506805
[11,] 5.7583641 -0.5429891
[12,] -6.2700392 5.7583641
[13,] -0.1946044 -6.2700392
[14,] 3.8886796 -0.1946044
[15,] -11.6739808 3.8886796
[16,] 15.1330095 -11.6739808
[17,] -0.1218054 15.1330095
[18,] -16.1661852 -0.1218054
[19,] -15.4368601 -16.1661852
[20,] 6.5676608 -15.4368601
[21,] -6.2358612 6.5676608
[22,] 5.2399001 -6.2358612
[23,] 14.2557175 5.2399001
[24,] -18.6183425 14.2557175
[25,] -23.0036633 -18.6183425
[26,] 5.2296552 -23.0036633
[27,] -1.1162405 5.2296552
[28,] 9.3175985 -1.1162405
[29,] 30.0716264 9.3175985
[30,] 48.7575842 30.0716264
[31,] -11.4150350 48.7575842
[32,] -10.9110457 -11.4150350
[33,] -6.9311354 -10.9110457
[34,] 2.2842257 -6.9311354
[35,] 3.4932963 2.2842257
[36,] 17.1936194 3.4932963
[37,] -5.1438353 17.1936194
[38,] -3.1180614 -5.1438353
[39,] -2.7179416 -3.1180614
[40,] -17.5125721 -2.7179416
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -15.6745359 13.7891430
2 -8.5420949 -15.6745359
3 -11.6941046 -8.5420949
4 -0.9914755 -11.6941046
5 4.6325896 -0.9914755
6 6.9211854 4.6325896
7 7.0728641 6.9211854
8 -1.6236301 7.0728641
9 -3.9506805 -1.6236301
10 -0.5429891 -3.9506805
11 5.7583641 -0.5429891
12 -6.2700392 5.7583641
13 -0.1946044 -6.2700392
14 3.8886796 -0.1946044
15 -11.6739808 3.8886796
16 15.1330095 -11.6739808
17 -0.1218054 15.1330095
18 -16.1661852 -0.1218054
19 -15.4368601 -16.1661852
20 6.5676608 -15.4368601
21 -6.2358612 6.5676608
22 5.2399001 -6.2358612
23 14.2557175 5.2399001
24 -18.6183425 14.2557175
25 -23.0036633 -18.6183425
26 5.2296552 -23.0036633
27 -1.1162405 5.2296552
28 9.3175985 -1.1162405
29 30.0716264 9.3175985
30 48.7575842 30.0716264
31 -11.4150350 48.7575842
32 -10.9110457 -11.4150350
33 -6.9311354 -10.9110457
34 2.2842257 -6.9311354
35 3.4932963 2.2842257
36 17.1936194 3.4932963
37 -5.1438353 17.1936194
38 -3.1180614 -5.1438353
39 -2.7179416 -3.1180614
40 -17.5125721 -2.7179416
> 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/7ixac1449264906.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/8o7ja1449264906.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/9b6zg1449264906.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/10mr891449264906.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, signif(mysum$coefficients[i,1],6), 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.row.start(a)
> a<-table.element(a, mywarning)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/11sos11449264906.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,formatC(signif(mysum$coefficients[i,1],5),format='g',flag='+'))
+ a<-table.element(a,formatC(signif(mysum$coefficients[i,2],5),format='g',flag=' '))
+ a<-table.element(a,formatC(signif(mysum$coefficients[i,3],4),format='e',flag='+'))
+ a<-table.element(a,formatC(signif(mysum$coefficients[i,4],4),format='g',flag=' '))
+ a<-table.element(a,formatC(signif(mysum$coefficients[i,4]/2,4),format='g',flag=' '))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/127m2u1449264906.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,formatC(signif(sqrt(mysum$r.squared),6),format='g',flag=' '))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a,formatC(signif(mysum$r.squared,6),format='g',flag=' '))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a,formatC(signif(mysum$adj.r.squared,6),format='g',flag=' '))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a,formatC(signif(mysum$fstatistic[1],6),format='g',flag=' '))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, signif(mysum$fstatistic[2],6))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, signif(mysum$fstatistic[3],6))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a,formatC(signif(1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]),6),format='g',flag=' '))
> 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,formatC(signif(mysum$sigma,6),format='g',flag=' '))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a,formatC(signif(sum(myerror*myerror),6),format='g',flag=' '))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/13hwqp1449264906.tab")
> if(n < 200) {
+ 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,formatC(signif(x[i],6),format='g',flag=' '))
+ a<-table.element(a,formatC(signif(x[i]-mysum$resid[i],6),format='g',flag=' '))
+ a<-table.element(a,formatC(signif(mysum$resid[i],6),format='g',flag=' '))
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/wessaorg/rcomp/tmp/14l5do1449264906.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,formatC(signif(gqarr[mypoint-kp3+1,1],6),format='g',flag=' '))
+ a<-table.element(a,formatC(signif(gqarr[mypoint-kp3+1,2],6),format='g',flag=' '))
+ a<-table.element(a,formatC(signif(gqarr[mypoint-kp3+1,3],6),format='g',flag=' '))
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/wessaorg/rcomp/tmp/15ylkb1449264906.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,signif(numsignificant1,6))
+ a<-table.element(a,formatC(signif(numsignificant1/numgqtests,6),format='g',flag=' '))
+ 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,signif(numsignificant5,6))
+ a<-table.element(a,signif(numsignificant5/numgqtests,6))
+ 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,signif(numsignificant10,6))
+ a<-table.element(a,signif(numsignificant10/numgqtests,6))
+ 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/16vbma1449264906.tab")
+ }
+ }
>
> try(system("convert tmp/12hko1449264906.ps tmp/12hko1449264906.png",intern=TRUE))
character(0)
> try(system("convert tmp/20alb1449264906.ps tmp/20alb1449264906.png",intern=TRUE))
character(0)
> try(system("convert tmp/33r231449264906.ps tmp/33r231449264906.png",intern=TRUE))
character(0)
> try(system("convert tmp/41c891449264906.ps tmp/41c891449264906.png",intern=TRUE))
character(0)
> try(system("convert tmp/545cx1449264906.ps tmp/545cx1449264906.png",intern=TRUE))
character(0)
> try(system("convert tmp/61kcq1449264906.ps tmp/61kcq1449264906.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ixac1449264906.ps tmp/7ixac1449264906.png",intern=TRUE))
character(0)
> try(system("convert tmp/8o7ja1449264906.ps tmp/8o7ja1449264906.png",intern=TRUE))
character(0)
> try(system("convert tmp/9b6zg1449264906.ps tmp/9b6zg1449264906.png",intern=TRUE))
character(0)
> try(system("convert tmp/10mr891449264906.ps tmp/10mr891449264906.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.128 0.676 4.855