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