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(2649.2,31077,2579.4,31293,2504.6,30236,2462.3,30160,2467.4,32436,2446.7,30695,2656.3,27525,2626.2,26434,2482.6,25739,2539.9,25204,2502.7,24977,2466.9,24320,2513.2,22680,2443.3,22052,2293.4,21467,2070.8,21383,2029.6,21777,2052,21928,1864.4,21814,1670.1,22937,1811,23595,1905.4,20830,1862.8,19650,2014.5,19195,2197.8,19644,2962.3,18483,3047,18079,3032.6,19178,3504.4,18391,3801.1,18441,3857.6,18584,3674.4,20108,3721,20148,3844.5,19394,4116.7,17745,4105.2,17696,4435.2,17032,4296.5,16438,4202.5,15683,4562.8,15594,4621.4,15713,4697,15937,4591.3,16171,4357,15928,4502.6,16348,4443.9,15579,4290.9,15305,4199.8,15648,4138.5,14954,3970.1,15137,3862.3,15839,3701.6,16050,3570.12,15168,3801.06,17064,3895.51,16005,3917.96,14886,3813.06,14931,3667.03,14544,3494.17,13812,3364,13031,3295.3,12574,3277.0,11964,3257.2,11451,3161.7,11346,3097.3,11353,3061.3,10702,3119.3,10646,3106.22,10556,3080.58,10463,2981.85,10407),dim=c(2,70),dimnames=list(c('Bel20','Goudprijs'),1:70)) > y <- array(NA,dim=c(2,70),dimnames=list(c('Bel20','Goudprijs'),1:70)) > for (i in 1:dim(x)[1]) + { + for (j in 1:dim(x)[2]) + { + y[i,j] <- as.numeric(x[i,j]) + } + } > par4 = 'no' > par3 = '' > par2 = 'none' > par1 = '1' > #'GNU S' R Code compiled by R2WASP v. 1.0.44 () > #Author: Dr. Ian E. Holliday > #To cite this work: Ian E. Holliday, 2009, 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: > #Technical description: > library(party) Loading required package: survival Loading required package: splines Loading required package: grid Loading required package: modeltools Loading required package: stats4 Loading required package: coin Loading required package: mvtnorm Loading required package: zoo Attaching package: 'zoo' The following object(s) are masked from package:base : as.Date.numeric Loading required package: sandwich Loading required package: strucchange Loading required package: vcd Loading required package: MASS Loading required package: colorspace > library(Hmisc) Attaching package: 'Hmisc' The following object(s) are masked from package:survival : untangle.specials The following object(s) are masked from package:base : format.pval, round.POSIXt, trunc.POSIXt, units > par1 <- as.numeric(par1) > par3 <- as.numeric(par3) > x <- data.frame(t(y)) > is.data.frame(x) [1] TRUE > x <- x[!is.na(x[,par1]),] > k <- length(x[1,]) > n <- length(x[,1]) > colnames(x)[par1] [1] "Bel20" > x[,par1] [1] 2649.20 2579.40 2504.60 2462.30 2467.40 2446.70 2656.30 2626.20 2482.60 [10] 2539.90 2502.70 2466.90 2513.20 2443.30 2293.40 2070.80 2029.60 2052.00 [19] 1864.40 1670.10 1811.00 1905.40 1862.80 2014.50 2197.80 2962.30 3047.00 [28] 3032.60 3504.40 3801.10 3857.60 3674.40 3721.00 3844.50 4116.70 4105.20 [37] 4435.20 4296.50 4202.50 4562.80 4621.40 4697.00 4591.30 4357.00 4502.60 [46] 4443.90 4290.90 4199.80 4138.50 3970.10 3862.30 3701.60 3570.12 3801.06 [55] 3895.51 3917.96 3813.06 3667.03 3494.17 3364.00 3295.30 3277.00 3257.20 [64] 3161.70 3097.30 3061.30 3119.30 3106.22 3080.58 2981.85 > if (par2 == 'kmeans') { + cl <- kmeans(x[,par1], par3) + print(cl) + clm <- matrix(cbind(cl$centers,1:par3),ncol=2) + clm <- clm[sort.list(clm[,1]),] + for (i in 1:par3) { + cl$cluster[cl$cluster==clm[i,2]] <- paste('C',i,sep='') + } + cl$cluster <- as.factor(cl$cluster) + print(cl$cluster) + x[,par1] <- cl$cluster + } > if (par2 == 'quantiles') { + x[,par1] <- cut2(x[,par1],g=par3) + } > if (par2 == 'hclust') { + hc <- hclust(dist(x[,par1])^2, 'cen') + print(hc) + memb <- cutree(hc, k = par3) + dum <- c(mean(x[memb==1,par1])) + for (i in 2:par3) { + dum <- c(dum, mean(x[memb==i,par1])) + } + hcm <- matrix(cbind(dum,1:par3),ncol=2) + hcm <- hcm[sort.list(hcm[,1]),] + for (i in 1:par3) { + memb[memb==hcm[i,2]] <- paste('C',i,sep='') + } + memb <- as.factor(memb) + print(memb) + x[,par1] <- memb + } > if (par2=='equal') { + ed <- cut(as.numeric(x[,par1]),par3,labels=paste('C',1:par3,sep='')) + x[,par1] <- as.factor(ed) + } > table(x[,par1]) 1670.1 1811 1862.8 1864.4 1905.4 2014.5 2029.6 2052 2070.8 2197.8 1 1 1 1 1 1 1 1 1 1 2293.4 2443.3 2446.7 2462.3 2466.9 2467.4 2482.6 2502.7 2504.6 2513.2 1 1 1 1 1 1 1 1 1 1 2539.9 2579.4 2626.2 2649.2 2656.3 2962.3 2981.85 3032.6 3047 3061.3 1 1 1 1 1 1 1 1 1 1 3080.58 3097.3 3106.22 3119.3 3161.7 3257.2 3277 3295.3 3364 3494.17 1 1 1 1 1 1 1 1 1 1 3504.4 3570.12 3667.03 3674.4 3701.6 3721 3801.06 3801.1 3813.06 3844.5 1 1 1 1 1 1 1 1 1 1 3857.6 3862.3 3895.51 3917.96 3970.1 4105.2 4116.7 4138.5 4199.8 4202.5 1 1 1 1 1 1 1 1 1 1 4290.9 4296.5 4357 4435.2 4443.9 4502.6 4562.8 4591.3 4621.4 4697 1 1 1 1 1 1 1 1 1 1 > colnames(x) [1] "Bel20" "Goudprijs" > colnames(x)[par1] [1] "Bel20" > x[,par1] [1] 2649.20 2579.40 2504.60 2462.30 2467.40 2446.70 2656.30 2626.20 2482.60 [10] 2539.90 2502.70 2466.90 2513.20 2443.30 2293.40 2070.80 2029.60 2052.00 [19] 1864.40 1670.10 1811.00 1905.40 1862.80 2014.50 2197.80 2962.30 3047.00 [28] 3032.60 3504.40 3801.10 3857.60 3674.40 3721.00 3844.50 4116.70 4105.20 [37] 4435.20 4296.50 4202.50 4562.80 4621.40 4697.00 4591.30 4357.00 4502.60 [46] 4443.90 4290.90 4199.80 4138.50 3970.10 3862.30 3701.60 3570.12 3801.06 [55] 3895.51 3917.96 3813.06 3667.03 3494.17 3364.00 3295.30 3277.00 3257.20 [64] 3161.70 3097.30 3061.30 3119.30 3106.22 3080.58 2981.85 > if (par2 == 'none') { + m <- ctree(as.formula(paste(colnames(x)[par1],' ~ .',sep='')),data = x) + } > > #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") > > if (par2 != 'none') { + m <- ctree(as.formula(paste('as.factor(',colnames(x)[par1],') ~ .',sep='')),data = x) + if (par4=='yes') { + a<-table.start() + a<-table.row.start(a) + a<-table.element(a,'10-Fold Cross Validation',3+2*par3,TRUE) + a<-table.row.end(a) + a<-table.row.start(a) + a<-table.element(a,'',1,TRUE) + a<-table.element(a,'Prediction (training)',par3+1,TRUE) + a<-table.element(a,'Prediction (testing)',par3+1,TRUE) + a<-table.row.end(a) + a<-table.row.start(a) + a<-table.element(a,'Actual',1,TRUE) + for (jjj in 1:par3) a<-table.element(a,paste('C',jjj,sep=''),1,TRUE) + a<-table.element(a,'CV',1,TRUE) + for (jjj in 1:par3) a<-table.element(a,paste('C',jjj,sep=''),1,TRUE) + a<-table.element(a,'CV',1,TRUE) + a<-table.row.end(a) + for (i in 1:10) { + ind <- sample(2, nrow(x), replace=T, prob=c(0.9,0.1)) + m.ct <- ctree(as.formula(paste('as.factor(',colnames(x)[par1],') ~ .',sep='')),data =x[ind==1,]) + if (i==1) { + m.ct.i.pred <- predict(m.ct, newdata=x[ind==1,]) + m.ct.i.actu <- x[ind==1,par1] + m.ct.x.pred <- predict(m.ct, newdata=x[ind==2,]) + m.ct.x.actu <- x[ind==2,par1] + } else { + m.ct.i.pred <- c(m.ct.i.pred,predict(m.ct, newdata=x[ind==1,])) + m.ct.i.actu <- c(m.ct.i.actu,x[ind==1,par1]) + m.ct.x.pred <- c(m.ct.x.pred,predict(m.ct, newdata=x[ind==2,])) + m.ct.x.actu <- c(m.ct.x.actu,x[ind==2,par1]) + } + } + print(m.ct.i.tab <- table(m.ct.i.actu,m.ct.i.pred)) + numer <- 0 + for (i in 1:par3) { + print(m.ct.i.tab[i,i] / sum(m.ct.i.tab[i,])) + numer <- numer + m.ct.i.tab[i,i] + } + print(m.ct.i.cp <- numer / sum(m.ct.i.tab)) + print(m.ct.x.tab <- table(m.ct.x.actu,m.ct.x.pred)) + numer <- 0 + for (i in 1:par3) { + print(m.ct.x.tab[i,i] / sum(m.ct.x.tab[i,])) + numer <- numer + m.ct.x.tab[i,i] + } + print(m.ct.x.cp <- numer / sum(m.ct.x.tab)) + for (i in 1:par3) { + a<-table.row.start(a) + a<-table.element(a,paste('C',i,sep=''),1,TRUE) + for (jjj in 1:par3) a<-table.element(a,m.ct.i.tab[i,jjj]) + a<-table.element(a,round(m.ct.i.tab[i,i]/sum(m.ct.i.tab[i,]),4)) + for (jjj in 1:par3) a<-table.element(a,m.ct.x.tab[i,jjj]) + a<-table.element(a,round(m.ct.x.tab[i,i]/sum(m.ct.x.tab[i,]),4)) + a<-table.row.end(a) + } + a<-table.row.start(a) + a<-table.element(a,'Overall',1,TRUE) + for (jjj in 1:par3) a<-table.element(a,'-') + a<-table.element(a,round(m.ct.i.cp,4)) + for (jjj in 1:par3) a<-table.element(a,'-') + a<-table.element(a,round(m.ct.x.cp,4)) + a<-table.row.end(a) + a<-table.end(a) + table.save(a,file="/var/www/html/freestat/rcomp/tmp/1962x1292274563.tab") + } + } > m Conditional inference tree with 5 terminal nodes Response: Bel20 Input: Goudprijs Number of observations: 70 1) Goudprijs <= 18584; criterion = 1, statistic = 20.113 2) Goudprijs <= 13812; criterion = 1, statistic = 12.432 3)* weights = 12 2) Goudprijs > 13812 4) Goudprijs <= 17032; criterion = 0.96, statistic = 4.211 5) Goudprijs <= 15305; criterion = 0.98, statistic = 5.373 6)* weights = 7 5) Goudprijs > 15305 7)* weights = 14 4) Goudprijs > 17032 8)* weights = 8 1) Goudprijs > 18584 9)* weights = 29 > postscript(file="/var/www/html/freestat/rcomp/tmp/2962x1292274563.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(m) > dev.off() null device 1 > postscript(file="/var/www/html/freestat/rcomp/tmp/3962x1292274563.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(x[,par1] ~ as.factor(where(m)),main='Response by Terminal Node',xlab='Terminal Node',ylab='Response') > dev.off() null device 1 > if (par2 == 'none') { + forec <- predict(m) + result <- as.data.frame(cbind(x[,par1],forec,x[,par1]-forec)) + colnames(result) <- c('Actuals','Forecasts','Residuals') + print(result) + } Actuals Forecasts Residuals 1 2649.20 2461.552 187.6482759 2 2579.40 2461.552 117.8482759 3 2504.60 2461.552 43.0482759 4 2462.30 2461.552 0.7482759 5 2467.40 2461.552 5.8482759 6 2446.70 2461.552 -14.8517241 7 2656.30 2461.552 194.7482759 8 2626.20 2461.552 164.6482759 9 2482.60 2461.552 21.0482759 10 2539.90 2461.552 78.3482759 11 2502.70 2461.552 41.1482759 12 2466.90 2461.552 5.3482759 13 2513.20 2461.552 51.6482759 14 2443.30 2461.552 -18.2517241 15 2293.40 2461.552 -168.1517241 16 2070.80 2461.552 -390.7517241 17 2029.60 2461.552 -431.9517241 18 2052.00 2461.552 -409.5517241 19 1864.40 2461.552 -597.1517241 20 1670.10 2461.552 -791.4517241 21 1811.00 2461.552 -650.5517241 22 1905.40 2461.552 -556.1517241 23 1862.80 2461.552 -598.7517241 24 2014.50 2461.552 -447.0517241 25 2197.80 2461.552 -263.7517241 26 2962.30 3649.420 -687.1200000 27 3047.00 3649.420 -602.4200000 28 3032.60 2461.552 571.0482759 29 3504.40 3649.420 -145.0200000 30 3801.10 3649.420 151.6800000 31 3857.60 3649.420 208.1800000 32 3674.40 2461.552 1212.8482759 33 3721.00 2461.552 1259.4482759 34 3844.50 2461.552 1382.9482759 35 4116.70 3649.420 467.2800000 36 4105.20 3649.420 455.7800000 37 4435.20 4312.101 123.0992857 38 4296.50 4312.101 -15.6007143 39 4202.50 4312.101 -109.6007143 40 4562.80 4312.101 250.6992857 41 4621.40 4312.101 309.2992857 42 4697.00 4312.101 384.8992857 43 4591.30 4312.101 279.1992857 44 4357.00 4312.101 44.8992857 45 4502.60 4312.101 190.4992857 46 4443.90 4312.101 131.7992857 47 4290.90 3909.667 381.2328571 48 4199.80 4312.101 -112.3007143 49 4138.50 3909.667 228.8328571 50 3970.10 3909.667 60.4328571 51 3862.30 4312.101 -449.8007143 52 3701.60 4312.101 -610.5007143 53 3570.12 3909.667 -339.5471429 54 3801.06 3649.420 151.6400000 55 3895.51 4312.101 -416.5907143 56 3917.96 3909.667 8.2928571 57 3813.06 3909.667 -96.6071429 58 3667.03 3909.667 -242.6371429 59 3494.17 3191.327 302.8433333 60 3364.00 3191.327 172.6733333 61 3295.30 3191.327 103.9733333 62 3277.00 3191.327 85.6733333 63 3257.20 3191.327 65.8733333 64 3161.70 3191.327 -29.6266667 65 3097.30 3191.327 -94.0266667 66 3061.30 3191.327 -130.0266667 67 3119.30 3191.327 -72.0266667 68 3106.22 3191.327 -85.1066667 69 3080.58 3191.327 -110.7466667 70 2981.85 3191.327 -209.4766667 > if (par2 != 'none') { + print(cbind(as.factor(x[,par1]),predict(m))) + myt <- table(as.factor(x[,par1]),predict(m)) + print(myt) + } > postscript(file="/var/www/html/freestat/rcomp/tmp/4ky101292274563.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > if(par2=='none') { + op <- par(mfrow=c(2,2)) + plot(density(result$Actuals),main='Kernel Density Plot of Actuals') + plot(density(result$Residuals),main='Kernel Density Plot of Residuals') + plot(result$Forecasts,result$Actuals,main='Actuals versus Predictions',xlab='Predictions',ylab='Actuals') + plot(density(result$Forecasts),main='Kernel Density Plot of Predictions') + par(op) + } > if(par2!='none') { + plot(myt,main='Confusion Matrix',xlab='Actual',ylab='Predicted') + } > dev.off() null device 1 > if (par2 == 'none') { + detcoef <- cor(result$Forecasts,result$Actuals) + a<-table.start() + a<-table.row.start(a) + a<-table.element(a,'Goodness of Fit',2,TRUE) + a<-table.row.end(a) + a<-table.row.start(a) + a<-table.element(a,'Correlation',1,TRUE) + a<-table.element(a,round(detcoef,4)) + a<-table.row.end(a) + a<-table.row.start(a) + a<-table.element(a,'R-squared',1,TRUE) + a<-table.element(a,round(detcoef*detcoef,4)) + a<-table.row.end(a) + a<-table.row.start(a) + a<-table.element(a,'RMSE',1,TRUE) + a<-table.element(a,round(sqrt(mean((result$Residuals)^2)),4)) + a<-table.row.end(a) + a<-table.end(a) + table.save(a,file="/var/www/html/freestat/rcomp/tmp/5y8zr1292274563.tab") + a<-table.start() + a<-table.row.start(a) + a<-table.element(a,'Actuals, Predictions, and Residuals',4,TRUE) + a<-table.row.end(a) + a<-table.row.start(a) + a<-table.element(a,'#',header=TRUE) + a<-table.element(a,'Actuals',header=TRUE) + a<-table.element(a,'Forecasts',header=TRUE) + a<-table.element(a,'Residuals',header=TRUE) + a<-table.row.end(a) + for (i in 1:length(result$Actuals)) { + a<-table.row.start(a) + a<-table.element(a,i,header=TRUE) + a<-table.element(a,result$Actuals[i]) + a<-table.element(a,result$Forecasts[i]) + a<-table.element(a,result$Residuals[i]) + a<-table.row.end(a) + } + a<-table.end(a) + table.save(a,file="/var/www/html/freestat/rcomp/tmp/6jqxf1292274563.tab") + } > if (par2 != 'none') { + a<-table.start() + a<-table.row.start(a) + a<-table.element(a,'Confusion Matrix (predicted in columns / actuals in rows)',par3+1,TRUE) + a<-table.row.end(a) + a<-table.row.start(a) + a<-table.element(a,'',1,TRUE) + for (i in 1:par3) { + a<-table.element(a,paste('C',i,sep=''),1,TRUE) + } + a<-table.row.end(a) + for (i in 1:par3) { + a<-table.row.start(a) + a<-table.element(a,paste('C',i,sep=''),1,TRUE) + for (j in 1:par3) { + a<-table.element(a,myt[i,j]) + } + a<-table.row.end(a) + } + a<-table.end(a) + table.save(a,file="/var/www/html/freestat/rcomp/tmp/7cze01292274563.tab") + } > > try(system("convert tmp/2962x1292274563.ps tmp/2962x1292274563.png",intern=TRUE)) character(0) > try(system("convert tmp/3962x1292274563.ps tmp/3962x1292274563.png",intern=TRUE)) character(0) > try(system("convert tmp/4ky101292274563.ps tmp/4ky101292274563.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.836 0.749 4.013