R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(6654000
+ ,5712000
+ ,-999.0
+ ,3.3
+ ,38.6
+ ,645.0
+ ,3
+ ,5
+ ,3
+ ,1000
+ ,6600
+ ,2.0
+ ,8.3
+ ,4.5
+ ,42.0
+ ,3
+ ,1
+ ,3
+ ,3385
+ ,44500
+ ,-999.0
+ ,12.5
+ ,14.0
+ ,60.0
+ ,1
+ ,1
+ ,1
+ ,.920
+ ,5700
+ ,-999.0
+ ,16.5
+ ,-999.0
+ ,25.0
+ ,5
+ ,2
+ ,3
+ ,2547000
+ ,4603000
+ ,1.8
+ ,3.9
+ ,69.0
+ ,624.0
+ ,3
+ ,5
+ ,4
+ ,10550
+ ,179500
+ ,.7
+ ,9.8
+ ,27.0
+ ,180.0
+ ,4
+ ,4
+ ,4
+ ,.023
+ ,.300
+ ,3.9
+ ,19.7
+ ,19.0
+ ,35.0
+ ,1
+ ,1
+ ,1
+ ,160000
+ ,169000
+ ,1.0
+ ,6.2
+ ,30.4
+ ,392.0
+ ,4
+ ,5
+ ,4
+ ,3300
+ ,25600
+ ,3.6
+ ,14.5
+ ,28.0
+ ,63.0
+ ,1
+ ,2
+ ,1
+ ,52160
+ ,440000
+ ,1.4
+ ,9.7
+ ,50.0
+ ,230.0
+ ,1
+ ,1
+ ,1
+ ,.425
+ ,6400
+ ,1.5
+ ,12.5
+ ,7.0
+ ,112.0
+ ,5
+ ,4
+ ,4
+ ,465000
+ ,423000
+ ,.7
+ ,3.9
+ ,30.0
+ ,281.0
+ ,5
+ ,5
+ ,5
+ ,.550
+ ,2400
+ ,2.7
+ ,10.3
+ ,-999.0
+ ,-999.0
+ ,2
+ ,1
+ ,2
+ ,187100
+ ,419000
+ ,-999.0
+ ,3.1
+ ,40.0
+ ,365.0
+ ,5
+ ,5
+ ,5
+ ,.075
+ ,1200
+ ,2.1
+ ,8.4
+ ,3.5
+ ,42.0
+ ,1
+ ,1
+ ,1
+ ,3000
+ ,25000
+ ,.0
+ ,8.6
+ ,50.0
+ ,28.0
+ ,2
+ ,2
+ ,2
+ ,.785
+ ,3500
+ ,4.1
+ ,10.7
+ ,6.0
+ ,42.0
+ ,2
+ ,2
+ ,2
+ ,.200
+ ,5000
+ ,1.2
+ ,10.7
+ ,10.4
+ ,120.0
+ ,2
+ ,2
+ ,2
+ ,1410
+ ,17500
+ ,1.3
+ ,6.1
+ ,34.0
+ ,-999.0
+ ,1
+ ,2
+ ,1
+ ,60000
+ ,81000
+ ,6.1
+ ,18.1
+ ,7.0
+ ,-999.0
+ ,1
+ ,1
+ ,1
+ ,529000
+ ,680000
+ ,.3
+ ,-999.0
+ ,28.0
+ ,400.0
+ ,5
+ ,5
+ ,5
+ ,27660
+ ,115000
+ ,.5
+ ,3.8
+ ,20.0
+ ,148.0
+ ,5
+ ,5
+ ,5
+ ,.120
+ ,1000
+ ,3.4
+ ,14.4
+ ,3.9
+ ,16.0
+ ,3
+ ,1
+ ,2
+ ,207000
+ ,406000
+ ,-999.0
+ ,12.0
+ ,39.3
+ ,252.0
+ ,1
+ ,4
+ ,1
+ ,85000
+ ,325000
+ ,1.5
+ ,6.2
+ ,41.0
+ ,310.0
+ ,1
+ ,3
+ ,1
+ ,36330
+ ,119500
+ ,-999.0
+ ,13.0
+ ,16.2
+ ,63.0
+ ,1
+ ,1
+ ,1
+ ,.101
+ ,4000
+ ,3.4
+ ,13.8
+ ,9.0
+ ,28.0
+ ,5
+ ,1
+ ,3
+ ,1040
+ ,5500
+ ,.8
+ ,8.2
+ ,7.6
+ ,68.0
+ ,5
+ ,3
+ ,4
+ ,521000
+ ,655000
+ ,.8
+ ,2.9
+ ,46.0
+ ,336.0
+ ,5
+ ,5
+ ,5
+ ,100000
+ ,157000
+ ,-999.0
+ ,10.8
+ ,22.4
+ ,100.0
+ ,1
+ ,1
+ ,1
+ ,35000
+ ,56000
+ ,-999.0
+ ,-999.0
+ ,16.3
+ ,33.0
+ ,3
+ ,5
+ ,4
+ ,.005
+ ,.140
+ ,1.4
+ ,9.1
+ ,2.6
+ ,21.5
+ ,5
+ ,2
+ ,4
+ ,.010
+ ,.250
+ ,2.0
+ ,19.9
+ ,24.0
+ ,50.0
+ ,1
+ ,1
+ ,1
+ ,62000
+ ,1320000
+ ,1.9
+ ,8.0
+ ,100.0
+ ,267.0
+ ,1
+ ,1
+ ,1
+ ,.122
+ ,3000
+ ,2.4
+ ,10.6
+ ,-999.0
+ ,30.0
+ ,2
+ ,1
+ ,1
+ ,1350
+ ,8100
+ ,2.8
+ ,11.2
+ ,-999.0
+ ,45.0
+ ,3
+ ,1
+ ,3
+ ,.023
+ ,.400
+ ,1.3
+ ,13.2
+ ,3.2
+ ,19.0
+ ,4
+ ,1
+ ,3
+ ,.048
+ ,.330
+ ,2.0
+ ,12.8
+ ,2.0
+ ,30.0
+ ,4
+ ,1
+ ,3
+ ,1700
+ ,6300
+ ,5.6
+ ,19.4
+ ,5.0
+ ,12.0
+ ,2
+ ,1
+ ,1
+ ,3500
+ ,10800
+ ,3.1
+ ,17.4
+ ,6.5
+ ,120.0
+ ,2
+ ,1
+ ,1
+ ,250000
+ ,490000
+ ,1.0
+ ,-999.0
+ ,23.6
+ ,440.0
+ ,5
+ ,5
+ ,5
+ ,.480
+ ,15500
+ ,1.8
+ ,17.0
+ ,12.0
+ ,140.0
+ ,2
+ ,2
+ ,2
+ ,10000
+ ,115000
+ ,.9
+ ,10.9
+ ,20.2
+ ,170.0
+ ,4
+ ,4
+ ,4
+ ,1620
+ ,11400
+ ,1.8
+ ,13.7
+ ,13.0
+ ,17.0
+ ,2
+ ,1
+ ,2
+ ,192000
+ ,180000
+ ,1.9
+ ,8.4
+ ,27.0
+ ,115.0
+ ,4
+ ,4
+ ,4
+ ,2500
+ ,12100
+ ,.9
+ ,8.4
+ ,18.0
+ ,31.0
+ ,5
+ ,5
+ ,5
+ ,4288
+ ,39200
+ ,-999.0
+ ,12.5
+ ,13.7
+ ,63.0
+ ,2
+ ,2
+ ,2
+ ,.280
+ ,1900
+ ,2.6
+ ,13.2
+ ,4.7
+ ,21.0
+ ,3
+ ,1
+ ,3
+ ,4235
+ ,50400
+ ,2.4
+ ,9.8
+ ,9.8
+ ,52.0
+ ,1
+ ,1
+ ,1
+ ,6800
+ ,179000
+ ,1.2
+ ,9.6
+ ,29.0
+ ,164.0
+ ,2
+ ,3
+ ,2
+ ,.750
+ ,12300
+ ,.9
+ ,6.6
+ ,7.0
+ ,225.0
+ ,2
+ ,2
+ ,2
+ ,3600
+ ,21000
+ ,.5
+ ,5.4
+ ,6.0
+ ,225.0
+ ,3
+ ,2
+ ,3
+ ,14830
+ ,98200
+ ,-999.0
+ ,2.6
+ ,17.0
+ ,150.0
+ ,5
+ ,5
+ ,5
+ ,55500
+ ,175000
+ ,.6
+ ,3.8
+ ,20.0
+ ,151.0
+ ,5
+ ,5
+ ,5
+ ,1400
+ ,12500
+ ,-999.0
+ ,11.0
+ ,12.7
+ ,90.0
+ ,2
+ ,2
+ ,2
+ ,.060
+ ,1000
+ ,2.2
+ ,10.3
+ ,3.5
+ ,-999.0
+ ,3
+ ,1
+ ,2
+ ,.900
+ ,2600
+ ,2.3
+ ,13.3
+ ,4.5
+ ,60.0
+ ,2
+ ,1
+ ,2
+ ,2000
+ ,12300
+ ,.5
+ ,5.4
+ ,7.5
+ ,200.0
+ ,3
+ ,1
+ ,3
+ ,.104
+ ,2500
+ ,2.6
+ ,15.8
+ ,2.3
+ ,46.0
+ ,3
+ ,2
+ ,2
+ ,4190
+ ,58000
+ ,.6
+ ,10.3
+ ,24.0
+ ,210.0
+ ,4
+ ,3
+ ,4
+ ,3500
+ ,3900
+ ,6.6
+ ,19.4
+ ,3.0
+ ,14.0
+ ,2
+ ,1
+ ,1
+ ,4050
+ ,17000
+ ,-999.0
+ ,-999.0
+ ,13.0
+ ,38.0
+ ,3
+ ,1
+ ,1)
+ ,dim=c(9
+ ,62)
+ ,dimnames=list(c('bowgth'
+ ,'brwght'
+ ,'PS'
+ ,'TS'
+ ,'LIFESPAN'
+ ,'DRAAGTIJD'
+ ,'PRED'
+ ,'Exposure'
+ ,'OverallD')
+ ,1:62))
> y <- array(NA,dim=c(9,62),dimnames=list(c('bowgth','brwght','PS','TS','LIFESPAN','DRAAGTIJD','PRED','Exposure','OverallD'),1:62))
> 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 = '3'
> #'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
PS bowgth brwght TS LIFESPAN DRAAGTIJD PRED Exposure OverallD
1 -999.0 6.654e+06 5.712e+06 3.3 38.6 645.0 3 5 3
2 2.0 1.000e+03 6.600e+03 8.3 4.5 42.0 3 1 3
3 -999.0 3.385e+03 4.450e+04 12.5 14.0 60.0 1 1 1
4 -999.0 9.200e-01 5.700e+03 16.5 -999.0 25.0 5 2 3
5 1.8 2.547e+06 4.603e+06 3.9 69.0 624.0 3 5 4
6 0.7 1.055e+04 1.795e+05 9.8 27.0 180.0 4 4 4
7 3.9 2.300e-02 3.000e-01 19.7 19.0 35.0 1 1 1
8 1.0 1.600e+05 1.690e+05 6.2 30.4 392.0 4 5 4
9 3.6 3.300e+03 2.560e+04 14.5 28.0 63.0 1 2 1
10 1.4 5.216e+04 4.400e+05 9.7 50.0 230.0 1 1 1
11 1.5 4.250e-01 6.400e+03 12.5 7.0 112.0 5 4 4
12 0.7 4.650e+05 4.230e+05 3.9 30.0 281.0 5 5 5
13 2.7 5.500e-01 2.400e+03 10.3 -999.0 -999.0 2 1 2
14 -999.0 1.871e+05 4.190e+05 3.1 40.0 365.0 5 5 5
15 2.1 7.500e-02 1.200e+03 8.4 3.5 42.0 1 1 1
16 0.0 3.000e+03 2.500e+04 8.6 50.0 28.0 2 2 2
17 4.1 7.850e-01 3.500e+03 10.7 6.0 42.0 2 2 2
18 1.2 2.000e-01 5.000e+03 10.7 10.4 120.0 2 2 2
19 1.3 1.410e+03 1.750e+04 6.1 34.0 -999.0 1 2 1
20 6.1 6.000e+04 8.100e+04 18.1 7.0 -999.0 1 1 1
21 0.3 5.290e+05 6.800e+05 -999.0 28.0 400.0 5 5 5
22 0.5 2.766e+04 1.150e+05 3.8 20.0 148.0 5 5 5
23 3.4 1.200e-01 1.000e+03 14.4 3.9 16.0 3 1 2
24 -999.0 2.070e+05 4.060e+05 12.0 39.3 252.0 1 4 1
25 1.5 8.500e+04 3.250e+05 6.2 41.0 310.0 1 3 1
26 -999.0 3.633e+04 1.195e+05 13.0 16.2 63.0 1 1 1
27 3.4 1.010e-01 4.000e+03 13.8 9.0 28.0 5 1 3
28 0.8 1.040e+03 5.500e+03 8.2 7.6 68.0 5 3 4
29 0.8 5.210e+05 6.550e+05 2.9 46.0 336.0 5 5 5
30 -999.0 1.000e+05 1.570e+05 10.8 22.4 100.0 1 1 1
31 -999.0 3.500e+04 5.600e+04 -999.0 16.3 33.0 3 5 4
32 1.4 5.000e-03 1.400e-01 9.1 2.6 21.5 5 2 4
33 2.0 1.000e-02 2.500e-01 19.9 24.0 50.0 1 1 1
34 1.9 6.200e+04 1.320e+06 8.0 100.0 267.0 1 1 1
35 2.4 1.220e-01 3.000e+03 10.6 -999.0 30.0 2 1 1
36 2.8 1.350e+03 8.100e+03 11.2 -999.0 45.0 3 1 3
37 1.3 2.300e-02 4.000e-01 13.2 3.2 19.0 4 1 3
38 2.0 4.800e-02 3.300e-01 12.8 2.0 30.0 4 1 3
39 5.6 1.700e+03 6.300e+03 19.4 5.0 12.0 2 1 1
40 3.1 3.500e+03 1.080e+04 17.4 6.5 120.0 2 1 1
41 1.0 2.500e+05 4.900e+05 -999.0 23.6 440.0 5 5 5
42 1.8 4.800e-01 1.550e+04 17.0 12.0 140.0 2 2 2
43 0.9 1.000e+04 1.150e+05 10.9 20.2 170.0 4 4 4
44 1.8 1.620e+03 1.140e+04 13.7 13.0 17.0 2 1 2
45 1.9 1.920e+05 1.800e+05 8.4 27.0 115.0 4 4 4
46 0.9 2.500e+03 1.210e+04 8.4 18.0 31.0 5 5 5
47 -999.0 4.288e+03 3.920e+04 12.5 13.7 63.0 2 2 2
48 2.6 2.800e-01 1.900e+03 13.2 4.7 21.0 3 1 3
49 2.4 4.235e+03 5.040e+04 9.8 9.8 52.0 1 1 1
50 1.2 6.800e+03 1.790e+05 9.6 29.0 164.0 2 3 2
51 0.9 7.500e-01 1.230e+04 6.6 7.0 225.0 2 2 2
52 0.5 3.600e+03 2.100e+04 5.4 6.0 225.0 3 2 3
53 -999.0 1.483e+04 9.820e+04 2.6 17.0 150.0 5 5 5
54 0.6 5.550e+04 1.750e+05 3.8 20.0 151.0 5 5 5
55 -999.0 1.400e+03 1.250e+04 11.0 12.7 90.0 2 2 2
56 2.2 6.000e-02 1.000e+03 10.3 3.5 -999.0 3 1 2
57 2.3 9.000e-01 2.600e+03 13.3 4.5 60.0 2 1 2
58 0.5 2.000e+03 1.230e+04 5.4 7.5 200.0 3 1 3
59 2.6 1.040e-01 2.500e+03 15.8 2.3 46.0 3 2 2
60 0.6 4.190e+03 5.800e+04 10.3 24.0 210.0 4 3 4
61 6.6 3.500e+03 3.900e+03 19.4 3.0 14.0 2 1 1
62 -999.0 4.050e+03 1.700e+04 -999.0 13.0 38.0 3 1 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) bowgth brwght TS LIFESPAN DRAAGTIJD
-2.445e+02 -2.426e-04 1.970e-04 3.015e-01 1.841e-01 -1.574e-01
PRED Exposure OverallD
-3.062e+01 -1.037e+02 1.606e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-882.10 -30.60 114.98 212.82 466.70
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.445e+02 1.161e+02 -2.106 0.0400 *
bowgth -2.426e-04 1.588e-04 -1.528 0.1325
brwght 1.970e-04 1.595e-04 1.235 0.2222
TS 3.015e-01 2.069e-01 1.458 0.1508
LIFESPAN 1.841e-01 2.092e-01 0.880 0.3828
DRAAGTIJD -1.574e-01 1.876e-01 -0.839 0.4052
PRED -3.062e+01 9.593e+01 -0.319 0.7509
Exposure -1.037e+02 6.146e+01 -1.687 0.0975 .
OverallD 1.606e+02 1.226e+02 1.310 0.1958
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 384.7 on 53 degrees of freedom
Multiple R-squared: 0.1911, Adjusted R-squared: 0.06901
F-statistic: 1.565 on 8 and 53 DF, p-value: 0.1577
> 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.7383381 0.5233238 0.26166192
[2,] 0.7061522 0.5876955 0.29384777
[3,] 0.9315396 0.1369209 0.06846044
[4,] 0.9050518 0.1898964 0.09494819
[5,] 0.8527027 0.2945946 0.14729729
[6,] 0.7906649 0.4186701 0.20933506
[7,] 0.7286534 0.5426933 0.27134664
[8,] 0.7051035 0.5897930 0.29489651
[9,] 0.6639531 0.6720939 0.33604694
[10,] 0.6072212 0.7855575 0.39277876
[11,] 0.5184376 0.9631248 0.48156239
[12,] 0.4635001 0.9270003 0.53649987
[13,] 0.5215836 0.9568327 0.47841636
[14,] 0.5419155 0.9161690 0.45808450
[15,] 0.7732644 0.4534712 0.22673561
[16,] 0.7121848 0.5756303 0.28781516
[17,] 0.6349682 0.7300635 0.36503176
[18,] 0.5700932 0.8598136 0.42990678
[19,] 0.8545142 0.2909716 0.14548580
[20,] 0.8724250 0.2551499 0.12757497
[21,] 0.8221442 0.3557116 0.17785580
[22,] 0.7828587 0.4342826 0.21714132
[23,] 0.7852611 0.4294777 0.21473887
[24,] 0.7904161 0.4191679 0.20958394
[25,] 0.7331899 0.5336202 0.26681012
[26,] 0.6623726 0.6752547 0.33762736
[27,] 0.6013232 0.7973536 0.39867680
[28,] 0.5218455 0.9563089 0.47815447
[29,] 0.4389334 0.8778668 0.56106659
[30,] 0.3942323 0.7884645 0.60576774
[31,] 0.3502266 0.7004533 0.64977336
[32,] 0.2782231 0.5564462 0.72177690
[33,] 0.2046691 0.4093382 0.79533089
[34,] 0.1395209 0.2790417 0.86047913
[35,] 0.3995833 0.7991667 0.60041665
[36,] 0.5730130 0.8539740 0.42698700
[37,] 0.4498349 0.8996699 0.55016506
[38,] 0.3144339 0.6288677 0.68556614
[39,] 0.2058304 0.4116608 0.79416959
> postscript(file="/var/www/html/freestat/rcomp/tmp/109ki1292352907.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/209ki1292352907.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/3t1jl1292352907.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/4t1jl1292352907.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/5t1jl1292352907.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 = 62
Frequency = 1
1 2 3 4 5 6
-43.384190 -37.536808 -785.632834 -694.108985 9.688109 127.573490
7 8 9 10 11 12
218.184921 303.708790 321.631648 169.668934 182.693527 180.681550
13 14 15 16 17 18
113.985538 -874.035024 223.511446 180.309460 197.588118 205.861979
19 20 21 22 23 24
154.706812 58.895675 466.697332 115.974920 119.505084 -470.713597
25 26 27 28 29 30
423.038443 -792.497505 20.675553 73.012380 154.680299 -779.089315
31 32 33 34 35 36
-485.791662 -35.897834 217.665595 -3.669238 436.075306 147.388685
37 38 39 40 41 42
-11.421590 -8.648316 248.720103 263.099698 444.238366 205.347913
43 44 45 46 47 48
139.692282 84.325757 162.890784 111.103130 -810.158754 -41.073867
49 50 51 52 53 54
215.139030 280.738545 222.516367 91.835164 -882.099559 111.482336
55 56 57 58 59 60
-800.712632 -40.176091 94.621377 -14.720994 226.674809 51.320075
61 62
251.312693 -417.093229
> postscript(file="/var/www/html/freestat/rcomp/tmp/6la061292352907.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 = 62
Frequency = 1
lag(myerror, k = 1) myerror
0 -43.384190 NA
1 -37.536808 -43.384190
2 -785.632834 -37.536808
3 -694.108985 -785.632834
4 9.688109 -694.108985
5 127.573490 9.688109
6 218.184921 127.573490
7 303.708790 218.184921
8 321.631648 303.708790
9 169.668934 321.631648
10 182.693527 169.668934
11 180.681550 182.693527
12 113.985538 180.681550
13 -874.035024 113.985538
14 223.511446 -874.035024
15 180.309460 223.511446
16 197.588118 180.309460
17 205.861979 197.588118
18 154.706812 205.861979
19 58.895675 154.706812
20 466.697332 58.895675
21 115.974920 466.697332
22 119.505084 115.974920
23 -470.713597 119.505084
24 423.038443 -470.713597
25 -792.497505 423.038443
26 20.675553 -792.497505
27 73.012380 20.675553
28 154.680299 73.012380
29 -779.089315 154.680299
30 -485.791662 -779.089315
31 -35.897834 -485.791662
32 217.665595 -35.897834
33 -3.669238 217.665595
34 436.075306 -3.669238
35 147.388685 436.075306
36 -11.421590 147.388685
37 -8.648316 -11.421590
38 248.720103 -8.648316
39 263.099698 248.720103
40 444.238366 263.099698
41 205.347913 444.238366
42 139.692282 205.347913
43 84.325757 139.692282
44 162.890784 84.325757
45 111.103130 162.890784
46 -810.158754 111.103130
47 -41.073867 -810.158754
48 215.139030 -41.073867
49 280.738545 215.139030
50 222.516367 280.738545
51 91.835164 222.516367
52 -882.099559 91.835164
53 111.482336 -882.099559
54 -800.712632 111.482336
55 -40.176091 -800.712632
56 94.621377 -40.176091
57 -14.720994 94.621377
58 226.674809 -14.720994
59 51.320075 226.674809
60 251.312693 51.320075
61 -417.093229 251.312693
62 NA -417.093229
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -37.536808 -43.384190
[2,] -785.632834 -37.536808
[3,] -694.108985 -785.632834
[4,] 9.688109 -694.108985
[5,] 127.573490 9.688109
[6,] 218.184921 127.573490
[7,] 303.708790 218.184921
[8,] 321.631648 303.708790
[9,] 169.668934 321.631648
[10,] 182.693527 169.668934
[11,] 180.681550 182.693527
[12,] 113.985538 180.681550
[13,] -874.035024 113.985538
[14,] 223.511446 -874.035024
[15,] 180.309460 223.511446
[16,] 197.588118 180.309460
[17,] 205.861979 197.588118
[18,] 154.706812 205.861979
[19,] 58.895675 154.706812
[20,] 466.697332 58.895675
[21,] 115.974920 466.697332
[22,] 119.505084 115.974920
[23,] -470.713597 119.505084
[24,] 423.038443 -470.713597
[25,] -792.497505 423.038443
[26,] 20.675553 -792.497505
[27,] 73.012380 20.675553
[28,] 154.680299 73.012380
[29,] -779.089315 154.680299
[30,] -485.791662 -779.089315
[31,] -35.897834 -485.791662
[32,] 217.665595 -35.897834
[33,] -3.669238 217.665595
[34,] 436.075306 -3.669238
[35,] 147.388685 436.075306
[36,] -11.421590 147.388685
[37,] -8.648316 -11.421590
[38,] 248.720103 -8.648316
[39,] 263.099698 248.720103
[40,] 444.238366 263.099698
[41,] 205.347913 444.238366
[42,] 139.692282 205.347913
[43,] 84.325757 139.692282
[44,] 162.890784 84.325757
[45,] 111.103130 162.890784
[46,] -810.158754 111.103130
[47,] -41.073867 -810.158754
[48,] 215.139030 -41.073867
[49,] 280.738545 215.139030
[50,] 222.516367 280.738545
[51,] 91.835164 222.516367
[52,] -882.099559 91.835164
[53,] 111.482336 -882.099559
[54,] -800.712632 111.482336
[55,] -40.176091 -800.712632
[56,] 94.621377 -40.176091
[57,] -14.720994 94.621377
[58,] 226.674809 -14.720994
[59,] 51.320075 226.674809
[60,] 251.312693 51.320075
[61,] -417.093229 251.312693
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -37.536808 -43.384190
2 -785.632834 -37.536808
3 -694.108985 -785.632834
4 9.688109 -694.108985
5 127.573490 9.688109
6 218.184921 127.573490
7 303.708790 218.184921
8 321.631648 303.708790
9 169.668934 321.631648
10 182.693527 169.668934
11 180.681550 182.693527
12 113.985538 180.681550
13 -874.035024 113.985538
14 223.511446 -874.035024
15 180.309460 223.511446
16 197.588118 180.309460
17 205.861979 197.588118
18 154.706812 205.861979
19 58.895675 154.706812
20 466.697332 58.895675
21 115.974920 466.697332
22 119.505084 115.974920
23 -470.713597 119.505084
24 423.038443 -470.713597
25 -792.497505 423.038443
26 20.675553 -792.497505
27 73.012380 20.675553
28 154.680299 73.012380
29 -779.089315 154.680299
30 -485.791662 -779.089315
31 -35.897834 -485.791662
32 217.665595 -35.897834
33 -3.669238 217.665595
34 436.075306 -3.669238
35 147.388685 436.075306
36 -11.421590 147.388685
37 -8.648316 -11.421590
38 248.720103 -8.648316
39 263.099698 248.720103
40 444.238366 263.099698
41 205.347913 444.238366
42 139.692282 205.347913
43 84.325757 139.692282
44 162.890784 84.325757
45 111.103130 162.890784
46 -810.158754 111.103130
47 -41.073867 -810.158754
48 215.139030 -41.073867
49 280.738545 215.139030
50 222.516367 280.738545
51 91.835164 222.516367
52 -882.099559 91.835164
53 111.482336 -882.099559
54 -800.712632 111.482336
55 -40.176091 -800.712632
56 94.621377 -40.176091
57 -14.720994 94.621377
58 226.674809 -14.720994
59 51.320075 226.674809
60 251.312693 51.320075
61 -417.093229 251.312693
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/7ejzr1292352907.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/8ejzr1292352907.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/9ejzr1292352907.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/freestat/rcomp/tmp/10pthc1292352907.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/11abfi1292352907.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/12vue61292352907.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/132vtz1292352907.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/14v4s21292352907.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/15gm881292352907.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/16k5pw1292352907.tab")
+ }
>
> try(system("convert tmp/109ki1292352907.ps tmp/109ki1292352907.png",intern=TRUE))
character(0)
> try(system("convert tmp/209ki1292352907.ps tmp/209ki1292352907.png",intern=TRUE))
character(0)
> try(system("convert tmp/3t1jl1292352907.ps tmp/3t1jl1292352907.png",intern=TRUE))
character(0)
> try(system("convert tmp/4t1jl1292352907.ps tmp/4t1jl1292352907.png",intern=TRUE))
character(0)
> try(system("convert tmp/5t1jl1292352907.ps tmp/5t1jl1292352907.png",intern=TRUE))
character(0)
> try(system("convert tmp/6la061292352907.ps tmp/6la061292352907.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ejzr1292352907.ps tmp/7ejzr1292352907.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ejzr1292352907.ps tmp/8ejzr1292352907.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ejzr1292352907.ps tmp/9ejzr1292352907.png",intern=TRUE))
character(0)
> try(system("convert tmp/10pthc1292352907.ps tmp/10pthc1292352907.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.986 2.495 4.339