R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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.
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(1000.00
+ ,6600.00
+ ,6.3
+ ,2.00
+ ,8.30
+ ,4.50
+ ,42.00
+ ,3.00
+ ,1.00
+ ,3.00
+ ,2547000.00
+ ,4603000.00
+ ,2.1
+ ,1.80
+ ,3.90
+ ,69.00
+ ,624.00
+ ,3.00
+ ,5.00
+ ,4.00
+ ,10550.00
+ ,179500.00
+ ,9.1
+ ,0.70
+ ,9.80
+ ,27.00
+ ,180.00
+ ,4.00
+ ,4.00
+ ,4.00
+ ,0.02
+ ,.300
+ ,15.8
+ ,3.90
+ ,19.70
+ ,19.00
+ ,35.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,160000.00
+ ,169000.00
+ ,5.2
+ ,1.00
+ ,6.20
+ ,30.40
+ ,392.00
+ ,4.00
+ ,5.00
+ ,4.00
+ ,3300.00
+ ,25600.00
+ ,10.9
+ ,3.60
+ ,14.50
+ ,28.00
+ ,63.00
+ ,1.00
+ ,2.00
+ ,1.00
+ ,52160.00
+ ,440000.00
+ ,8.3
+ ,1.40
+ ,9.70
+ ,50.00
+ ,230.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,0.43
+ ,6400.00
+ ,11.0
+ ,1.40
+ ,12.50
+ ,7.00
+ ,112.00
+ ,5.00
+ ,4.00
+ ,4.00
+ ,465000.00
+ ,423000.00
+ ,3.2
+ ,0.70
+ ,3.90
+ ,30.00
+ ,281.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,0.75
+ ,1200.00
+ ,6.3
+ ,2.10
+ ,8.40
+ ,3.50
+ ,42.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,0.79
+ ,3500.00
+ ,6.6
+ ,4.10
+ ,10.70
+ ,6.00
+ ,42.00
+ ,2.00
+ ,2.00
+ ,2.00
+ ,0.20
+ ,5000.00
+ ,9.5
+ ,1.20
+ ,10.70
+ ,10.40
+ ,120.00
+ ,2.00
+ ,2.00
+ ,2.00
+ ,27660.00
+ ,115000.00
+ ,3.3
+ ,0.50
+ ,3.80
+ ,20.00
+ ,148.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,0.12
+ ,1000.00
+ ,11.0
+ ,3.40
+ ,14.40
+ ,3.90
+ ,16.00
+ ,3.00
+ ,1.00
+ ,2.00
+ ,85000.00
+ ,325000.00
+ ,4.7
+ ,1.50
+ ,6.20
+ ,41.00
+ ,310.00
+ ,1.00
+ ,3.00
+ ,1.00
+ ,0.10
+ ,4000.00
+ ,10.4
+ ,3.40
+ ,13.80
+ ,9.00
+ ,28.00
+ ,5.00
+ ,1.00
+ ,3.00
+ ,1040.00
+ ,5500.00
+ ,7.4
+ ,0.80
+ ,8.20
+ ,7.60
+ ,68.00
+ ,5.00
+ ,3.00
+ ,4.00
+ ,521000.00
+ ,655000.00
+ ,2.1
+ ,0.80
+ ,2.90
+ ,46.00
+ ,336.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,0.10
+ ,0.25
+ ,17.9
+ ,2.00
+ ,19.90
+ ,24.00
+ ,50.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,62000.00
+ ,1320000.00
+ ,6.1
+ ,1.90
+ ,8.00
+ ,100.00
+ ,267.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,0.23
+ ,0.40
+ ,11.9
+ ,1.30
+ ,13.20
+ ,3.20
+ ,19.00
+ ,4.00
+ ,1.00
+ ,3.00
+ ,1700.00
+ ,6300.00
+ ,13.8
+ ,5.60
+ ,19.40
+ ,5.00
+ ,12.00
+ ,2.00
+ ,1.00
+ ,1.00
+ ,3500.00
+ ,10800.00
+ ,14.3
+ ,3.10
+ ,17.40
+ ,6.50
+ ,120.00
+ ,2.00
+ ,1.00
+ ,1.00
+ ,0.48
+ ,15500.00
+ ,15.2
+ ,1.80
+ ,17.00
+ ,12.00
+ ,140.00
+ ,2.00
+ ,2.00
+ ,2.00
+ ,10000.00
+ ,115000.00
+ ,10.0
+ ,0.90
+ ,10.90
+ ,20.20
+ ,170.00
+ ,4.00
+ ,4.00
+ ,4.00
+ ,1620.00
+ ,11400.00
+ ,11.9
+ ,1.80
+ ,13.70
+ ,13.00
+ ,17.00
+ ,2.00
+ ,1.00
+ ,2.00
+ ,192000.00
+ ,180000.00
+ ,6.5
+ ,1.90
+ ,8.40
+ ,27.00
+ ,115.00
+ ,4.00
+ ,4.00
+ ,4.00
+ ,2500.00
+ ,12100.00
+ ,7.5
+ ,0.90
+ ,8.40
+ ,18.00
+ ,31.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,0.28
+ ,1900.00
+ ,10.6
+ ,2.60
+ ,13.20
+ ,4.70
+ ,21.00
+ ,3.00
+ ,1.00
+ ,3.00
+ ,4235.00
+ ,50400.00
+ ,7.4
+ ,2.40
+ ,9.80
+ ,9.80
+ ,52.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,6800.00
+ ,179000.00
+ ,8.4
+ ,1.20
+ ,9.60
+ ,29.00
+ ,164.00
+ ,2.00
+ ,3.00
+ ,2.00
+ ,0.75
+ ,12300.00
+ ,5.7
+ ,0.90
+ ,6.60
+ ,7.00
+ ,225.00
+ ,2.00
+ ,2.00
+ ,2.00
+ ,3600.00
+ ,21000.00
+ ,4.9
+ ,0.50
+ ,5.40
+ ,6.00
+ ,225.00
+ ,3.00
+ ,2.00
+ ,3.00
+ ,55500.00
+ ,175000.00
+ ,3.2
+ ,0.60
+ ,3.80
+ ,20.00
+ ,151.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,0.90
+ ,2600.00
+ ,11.0
+ ,2.30
+ ,13.30
+ ,4.50
+ ,60.00
+ ,2.00
+ ,1.00
+ ,2.00
+ ,2000.00
+ ,12300.00
+ ,4.9
+ ,0.50
+ ,5.40
+ ,7.50
+ ,200.00
+ ,3.00
+ ,1.00
+ ,3.00
+ ,0.10
+ ,2500.00
+ ,13.2
+ ,2.60
+ ,15.80
+ ,2.30
+ ,46.00
+ ,3.00
+ ,2.00
+ ,2.00
+ ,4190.00
+ ,58000.00
+ ,9.7
+ ,0.60
+ ,10.30
+ ,24.00
+ ,210.00
+ ,4.00
+ ,3.00
+ ,4.00
+ ,3500.00
+ ,3900.00
+ ,12.8
+ ,6.60
+ ,19.40
+ ,3.00
+ ,14.00
+ ,2.00
+ ,1.00
+ ,1.00)
+ ,dim=c(10
+ ,39)
+ ,dimnames=list(c('body'
+ ,'brain'
+ ,'slowwave'
+ ,'paradoxical'
+ ,'total_sleep'
+ ,'lifespan'
+ ,'gestation'
+ ,'predation'
+ ,'sleepexp.'
+ ,'danger')
+ ,1:39))
> y <- array(NA,dim=c(10,39),dimnames=list(c('body','brain','slowwave','paradoxical','total_sleep','lifespan','gestation','predation','sleepexp.','danger'),1:39))
> 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
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
total_sleep body brain slowwave paradoxical lifespan gestation
1 8.3 1.000e+03 6.600e+03 6.3 2.0 4.5 42
2 3.9 2.547e+06 4.603e+06 2.1 1.8 69.0 624
3 9.8 1.055e+04 1.795e+05 9.1 0.7 27.0 180
4 19.7 2.000e-02 3.000e-01 15.8 3.9 19.0 35
5 6.2 1.600e+05 1.690e+05 5.2 1.0 30.4 392
6 14.5 3.300e+03 2.560e+04 10.9 3.6 28.0 63
7 9.7 5.216e+04 4.400e+05 8.3 1.4 50.0 230
8 12.5 4.300e-01 6.400e+03 11.0 1.4 7.0 112
9 3.9 4.650e+05 4.230e+05 3.2 0.7 30.0 281
10 8.4 7.500e-01 1.200e+03 6.3 2.1 3.5 42
11 10.7 7.900e-01 3.500e+03 6.6 4.1 6.0 42
12 10.7 2.000e-01 5.000e+03 9.5 1.2 10.4 120
13 3.8 2.766e+04 1.150e+05 3.3 0.5 20.0 148
14 14.4 1.200e-01 1.000e+03 11.0 3.4 3.9 16
15 6.2 8.500e+04 3.250e+05 4.7 1.5 41.0 310
16 13.8 1.000e-01 4.000e+03 10.4 3.4 9.0 28
17 8.2 1.040e+03 5.500e+03 7.4 0.8 7.6 68
18 2.9 5.210e+05 6.550e+05 2.1 0.8 46.0 336
19 19.9 1.000e-01 2.500e-01 17.9 2.0 24.0 50
20 8.0 6.200e+04 1.320e+06 6.1 1.9 100.0 267
21 13.2 2.300e-01 4.000e-01 11.9 1.3 3.2 19
22 19.4 1.700e+03 6.300e+03 13.8 5.6 5.0 12
23 17.4 3.500e+03 1.080e+04 14.3 3.1 6.5 120
24 17.0 4.800e-01 1.550e+04 15.2 1.8 12.0 140
25 10.9 1.000e+04 1.150e+05 10.0 0.9 20.2 170
26 13.7 1.620e+03 1.140e+04 11.9 1.8 13.0 17
27 8.4 1.920e+05 1.800e+05 6.5 1.9 27.0 115
28 8.4 2.500e+03 1.210e+04 7.5 0.9 18.0 31
29 13.2 2.800e-01 1.900e+03 10.6 2.6 4.7 21
30 9.8 4.235e+03 5.040e+04 7.4 2.4 9.8 52
31 9.6 6.800e+03 1.790e+05 8.4 1.2 29.0 164
32 6.6 7.500e-01 1.230e+04 5.7 0.9 7.0 225
33 5.4 3.600e+03 2.100e+04 4.9 0.5 6.0 225
34 3.8 5.550e+04 1.750e+05 3.2 0.6 20.0 151
35 13.3 9.000e-01 2.600e+03 11.0 2.3 4.5 60
36 5.4 2.000e+03 1.230e+04 4.9 0.5 7.5 200
37 15.8 1.000e-01 2.500e+03 13.2 2.6 2.3 46
38 10.3 4.190e+03 5.800e+04 9.7 0.6 24.0 210
39 19.4 3.500e+03 3.900e+03 12.8 6.6 3.0 14
predation sleepexp. danger
1 3 1 3
2 3 5 4
3 4 4 4
4 1 1 1
5 4 5 4
6 1 2 1
7 1 1 1
8 5 4 4
9 5 5 5
10 1 1 1
11 2 2 2
12 2 2 2
13 5 5 5
14 3 1 2
15 1 3 1
16 5 1 3
17 5 3 4
18 5 5 5
19 1 1 1
20 1 1 1
21 4 1 3
22 2 1 1
23 2 1 1
24 2 2 2
25 4 4 4
26 2 1 2
27 4 4 4
28 5 5 5
29 3 1 3
30 1 1 1
31 2 3 2
32 2 2 2
33 3 2 3
34 5 5 5
35 2 1 2
36 3 1 3
37 3 2 2
38 4 3 4
39 2 1 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) body brain slowwave paradoxical lifespan
-9.002e-03 -3.724e-08 2.462e-08 1.001e+00 9.980e-01 -3.258e-04
gestation predation sleepexp. danger
5.946e-06 9.547e-03 5.910e-03 -1.152e-02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.015027 -0.007351 -0.001416 0.004266 0.077772
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -9.002e-03 2.036e-02 -0.442 0.662
body -3.724e-08 3.763e-08 -0.990 0.330
brain 2.462e-08 2.149e-08 1.145 0.261
slowwave 1.001e+00 1.001e-03 1000.143 <2e-16 ***
paradoxical 9.980e-01 3.398e-03 293.696 <2e-16 ***
lifespan -3.258e-04 3.019e-04 -1.079 0.289
gestation 5.946e-06 5.041e-05 0.118 0.907
predation 9.547e-03 6.910e-03 1.382 0.178
sleepexp. 5.910e-03 3.919e-03 1.508 0.142
danger -1.152e-02 9.730e-03 -1.184 0.246
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.01638 on 29 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 3.721e+05 on 9 and 29 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,] 1 2.115945e-182 1.057972e-182
[2,] 1 6.269453e-191 3.134727e-191
[3,] 1 5.224817e-181 2.612409e-181
[4,] 1 7.318760e-163 3.659380e-163
[5,] 1 3.629650e-155 1.814825e-155
[6,] 1 4.302678e-139 2.151339e-139
[7,] 1 3.243008e-131 1.621504e-131
[8,] 1 1.633725e-115 8.168624e-116
[9,] 1 3.344468e-103 1.672234e-103
[10,] 1 2.658167e-90 1.329084e-90
[11,] 1 1.665265e-75 8.326327e-76
[12,] 1 6.381039e-63 3.190519e-63
[13,] 1 9.840384e-53 4.920192e-53
[14,] 1 1.431686e-40 7.158431e-41
> postscript(file="/var/www/html/rcomp/tmp/11pz91292173968.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/rcomp/tmp/21pz91292173968.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/rcomp/tmp/3chzu1292173968.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/rcomp/tmp/4chzu1292173968.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/rcomp/tmp/5chzu1292173968.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 = 39
Frequency = 1
1 2 3 4 5
8.279971e-03 -1.105590e-03 -1.008044e-02 4.227178e-03 -6.093206e-03
6 7 8 9 10
4.522348e-03 6.213195e-03 7.777236e-02 2.770280e-03 4.303858e-03
11 12 13 14 15
4.872104e-03 -2.723407e-03 -8.899281e-03 -4.727835e-03 -1.416421e-03
16 17 18 19 20
-1.022785e-02 -1.366495e-02 5.254062e-03 -6.128862e-06 4.038525e-03
21 22 23 24 25
-8.037170e-03 -4.558932e-03 -1.025122e-02 -6.664393e-03 -1.110263e-02
26 27 28 29 30
3.527591e-03 1.880593e-03 -1.035408e-02 5.763385e-03 4.826174e-03
31 32 33 34 35
-5.842511e-03 -1.710707e-03 -2.052084e-04 -9.063121e-03 2.500407e-03
36 37 38 39
6.496328e-03 -1.502662e-02 -3.331131e-03 -2.155528e-03
> postscript(file="/var/www/html/rcomp/tmp/64qgf1292173968.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 = 39
Frequency = 1
lag(myerror, k = 1) myerror
0 8.279971e-03 NA
1 -1.105590e-03 8.279971e-03
2 -1.008044e-02 -1.105590e-03
3 4.227178e-03 -1.008044e-02
4 -6.093206e-03 4.227178e-03
5 4.522348e-03 -6.093206e-03
6 6.213195e-03 4.522348e-03
7 7.777236e-02 6.213195e-03
8 2.770280e-03 7.777236e-02
9 4.303858e-03 2.770280e-03
10 4.872104e-03 4.303858e-03
11 -2.723407e-03 4.872104e-03
12 -8.899281e-03 -2.723407e-03
13 -4.727835e-03 -8.899281e-03
14 -1.416421e-03 -4.727835e-03
15 -1.022785e-02 -1.416421e-03
16 -1.366495e-02 -1.022785e-02
17 5.254062e-03 -1.366495e-02
18 -6.128862e-06 5.254062e-03
19 4.038525e-03 -6.128862e-06
20 -8.037170e-03 4.038525e-03
21 -4.558932e-03 -8.037170e-03
22 -1.025122e-02 -4.558932e-03
23 -6.664393e-03 -1.025122e-02
24 -1.110263e-02 -6.664393e-03
25 3.527591e-03 -1.110263e-02
26 1.880593e-03 3.527591e-03
27 -1.035408e-02 1.880593e-03
28 5.763385e-03 -1.035408e-02
29 4.826174e-03 5.763385e-03
30 -5.842511e-03 4.826174e-03
31 -1.710707e-03 -5.842511e-03
32 -2.052084e-04 -1.710707e-03
33 -9.063121e-03 -2.052084e-04
34 2.500407e-03 -9.063121e-03
35 6.496328e-03 2.500407e-03
36 -1.502662e-02 6.496328e-03
37 -3.331131e-03 -1.502662e-02
38 -2.155528e-03 -3.331131e-03
39 NA -2.155528e-03
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.105590e-03 8.279971e-03
[2,] -1.008044e-02 -1.105590e-03
[3,] 4.227178e-03 -1.008044e-02
[4,] -6.093206e-03 4.227178e-03
[5,] 4.522348e-03 -6.093206e-03
[6,] 6.213195e-03 4.522348e-03
[7,] 7.777236e-02 6.213195e-03
[8,] 2.770280e-03 7.777236e-02
[9,] 4.303858e-03 2.770280e-03
[10,] 4.872104e-03 4.303858e-03
[11,] -2.723407e-03 4.872104e-03
[12,] -8.899281e-03 -2.723407e-03
[13,] -4.727835e-03 -8.899281e-03
[14,] -1.416421e-03 -4.727835e-03
[15,] -1.022785e-02 -1.416421e-03
[16,] -1.366495e-02 -1.022785e-02
[17,] 5.254062e-03 -1.366495e-02
[18,] -6.128862e-06 5.254062e-03
[19,] 4.038525e-03 -6.128862e-06
[20,] -8.037170e-03 4.038525e-03
[21,] -4.558932e-03 -8.037170e-03
[22,] -1.025122e-02 -4.558932e-03
[23,] -6.664393e-03 -1.025122e-02
[24,] -1.110263e-02 -6.664393e-03
[25,] 3.527591e-03 -1.110263e-02
[26,] 1.880593e-03 3.527591e-03
[27,] -1.035408e-02 1.880593e-03
[28,] 5.763385e-03 -1.035408e-02
[29,] 4.826174e-03 5.763385e-03
[30,] -5.842511e-03 4.826174e-03
[31,] -1.710707e-03 -5.842511e-03
[32,] -2.052084e-04 -1.710707e-03
[33,] -9.063121e-03 -2.052084e-04
[34,] 2.500407e-03 -9.063121e-03
[35,] 6.496328e-03 2.500407e-03
[36,] -1.502662e-02 6.496328e-03
[37,] -3.331131e-03 -1.502662e-02
[38,] -2.155528e-03 -3.331131e-03
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.105590e-03 8.279971e-03
2 -1.008044e-02 -1.105590e-03
3 4.227178e-03 -1.008044e-02
4 -6.093206e-03 4.227178e-03
5 4.522348e-03 -6.093206e-03
6 6.213195e-03 4.522348e-03
7 7.777236e-02 6.213195e-03
8 2.770280e-03 7.777236e-02
9 4.303858e-03 2.770280e-03
10 4.872104e-03 4.303858e-03
11 -2.723407e-03 4.872104e-03
12 -8.899281e-03 -2.723407e-03
13 -4.727835e-03 -8.899281e-03
14 -1.416421e-03 -4.727835e-03
15 -1.022785e-02 -1.416421e-03
16 -1.366495e-02 -1.022785e-02
17 5.254062e-03 -1.366495e-02
18 -6.128862e-06 5.254062e-03
19 4.038525e-03 -6.128862e-06
20 -8.037170e-03 4.038525e-03
21 -4.558932e-03 -8.037170e-03
22 -1.025122e-02 -4.558932e-03
23 -6.664393e-03 -1.025122e-02
24 -1.110263e-02 -6.664393e-03
25 3.527591e-03 -1.110263e-02
26 1.880593e-03 3.527591e-03
27 -1.035408e-02 1.880593e-03
28 5.763385e-03 -1.035408e-02
29 4.826174e-03 5.763385e-03
30 -5.842511e-03 4.826174e-03
31 -1.710707e-03 -5.842511e-03
32 -2.052084e-04 -1.710707e-03
33 -9.063121e-03 -2.052084e-04
34 2.500407e-03 -9.063121e-03
35 6.496328e-03 2.500407e-03
36 -1.502662e-02 6.496328e-03
37 -3.331131e-03 -1.502662e-02
38 -2.155528e-03 -3.331131e-03
> 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/rcomp/tmp/7fhf01292173968.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/rcomp/tmp/8fhf01292173968.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/rcomp/tmp/9fhf01292173968.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')
Warning messages:
1: In sqrt(crit * p * (1 - hh)/hh) : NaNs produced
2: In sqrt(crit * p * (1 - hh)/hh) : NaNs produced
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/10qrw31292173968.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/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/rcomp/tmp/11brv81292173968.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/rcomp/tmp/12erte1292173968.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/rcomp/tmp/13bj9n1292173968.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/rcomp/tmp/14ekqt1292173968.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/rcomp/tmp/15i3oz1292173968.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/rcomp/tmp/16vum71292173968.tab")
+ }
>
> try(system("convert tmp/11pz91292173968.ps tmp/11pz91292173968.png",intern=TRUE))
character(0)
> try(system("convert tmp/21pz91292173968.ps tmp/21pz91292173968.png",intern=TRUE))
character(0)
> try(system("convert tmp/3chzu1292173968.ps tmp/3chzu1292173968.png",intern=TRUE))
character(0)
> try(system("convert tmp/4chzu1292173968.ps tmp/4chzu1292173968.png",intern=TRUE))
character(0)
> try(system("convert tmp/5chzu1292173968.ps tmp/5chzu1292173968.png",intern=TRUE))
character(0)
> try(system("convert tmp/64qgf1292173968.ps tmp/64qgf1292173968.png",intern=TRUE))
character(0)
> try(system("convert tmp/7fhf01292173968.ps tmp/7fhf01292173968.png",intern=TRUE))
character(0)
> try(system("convert tmp/8fhf01292173968.ps tmp/8fhf01292173968.png",intern=TRUE))
character(0)
> try(system("convert tmp/9fhf01292173968.ps tmp/9fhf01292173968.png",intern=TRUE))
character(0)
> try(system("convert tmp/10qrw31292173968.ps tmp/10qrw31292173968.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.340 1.658 5.513