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(235.1 + ,1 + ,280.7 + ,1 + ,264.6 + ,2 + ,240.7 + ,0 + ,201.4 + ,1 + ,240.8 + ,0 + ,241.1 + ,-1 + ,223.8 + ,-3 + ,206.1 + ,-3 + ,174.7 + ,-3 + ,203.3 + ,-4 + ,220.5 + ,-8 + ,299.5 + ,-9 + ,347.4 + ,-13 + ,338.3 + ,-18 + ,327.7 + ,-11 + ,351.6 + ,-9 + ,396.6 + ,-10 + ,438.8 + ,-13 + ,395.6 + ,-11 + ,363.5 + ,-5 + ,378.8 + ,-15 + ,357 + ,-6 + ,369 + ,-6 + ,464.8 + ,-3 + ,479.1 + ,-1 + ,431.3 + ,-3 + ,366.5 + ,-4 + ,326.3 + ,-6 + ,355.1 + ,0 + ,331.6 + ,-4 + ,261.3 + ,-2 + ,249 + ,-2 + ,205.5 + ,-6 + ,235.6 + ,-7 + ,240.9 + ,-6 + ,264.9 + ,-6 + ,253.8 + ,-3 + ,232.3 + ,-2 + ,193.8 + ,-5 + ,177 + ,-11 + ,213.2 + ,-11 + ,207.2 + ,-11 + ,180.6 + ,-10 + ,188.6 + ,-14 + ,175.4 + ,-8 + ,199 + ,-9 + ,179.6 + ,-5 + ,225.8 + ,-1 + ,234 + ,-2 + ,200.2 + ,-5 + ,183.6 + ,-4 + ,178.2 + ,-6 + ,203.2 + ,-2 + ,208.5 + ,-2 + ,191.8 + ,-2 + ,172.8 + ,-2 + ,148 + ,2 + ,159.4 + ,1 + ,154.5 + ,-8 + ,213.2 + ,-1 + ,196.4 + ,1 + ,182.8 + ,-1 + ,176.4 + ,2 + ,153.6 + ,2 + ,173.2 + ,1 + ,171 + ,-1 + ,151.2 + ,-2 + ,161.9 + ,-2 + ,157.2 + ,-1 + ,201.7 + ,-8 + ,236.4 + ,-4 + ,356.1 + ,-6 + ,398.3 + ,-3 + ,403.7 + ,-3 + ,384.6 + ,-7 + ,365.8 + ,-9 + ,368.1 + ,-11 + ,367.9 + ,-13 + ,347 + ,-11 + ,343.3 + ,-9 + ,292.9 + ,-17 + ,311.5 + ,-22 + ,300.9 + ,-25 + ,366.9 + ,-20 + ,356.9 + ,-24 + ,329.7 + ,-24 + ,316.2 + ,-22 + ,269 + ,-19 + ,289.3 + ,-18 + ,266.2 + ,-17 + ,253.6 + ,-11 + ,233.8 + ,-11 + ,228.4 + ,-12 + ,253.6 + ,-10 + ,260.1 + ,-15 + ,306.6 + ,-15 + ,309.2 + ,-15 + ,309.5 + ,-13 + ,271 + ,-8 + ,279.9 + ,-13 + ,317.9 + ,-9 + ,298.4 + ,-7 + ,246.7 + ,-4 + ,227.3 + ,-4 + ,209.1 + ,-2) + ,dim=c(2 + ,106) + ,dimnames=list(c('Werkloosheid' + ,'Consumentenvertrouwen') + ,1:106)) > y <- array(NA,dim=c(2,106),dimnames=list(c('Werkloosheid','Consumentenvertrouwen'),1:106)) > 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 = '3' > 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] "Werkloosheid" > x[,par1] [1] 235.1 280.7 264.6 240.7 201.4 240.8 241.1 223.8 206.1 174.7 203.3 220.5 [13] 299.5 347.4 338.3 327.7 351.6 396.6 438.8 395.6 363.5 378.8 357.0 369.0 [25] 464.8 479.1 431.3 366.5 326.3 355.1 331.6 261.3 249.0 205.5 235.6 240.9 [37] 264.9 253.8 232.3 193.8 177.0 213.2 207.2 180.6 188.6 175.4 199.0 179.6 [49] 225.8 234.0 200.2 183.6 178.2 203.2 208.5 191.8 172.8 148.0 159.4 154.5 [61] 213.2 196.4 182.8 176.4 153.6 173.2 171.0 151.2 161.9 157.2 201.7 236.4 [73] 356.1 398.3 403.7 384.6 365.8 368.1 367.9 347.0 343.3 292.9 311.5 300.9 [85] 366.9 356.9 329.7 316.2 269.0 289.3 266.2 253.6 233.8 228.4 253.6 260.1 [97] 306.6 309.2 309.5 271.0 279.9 317.9 298.4 246.7 227.3 209.1 > 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]) 148 151.2 153.6 154.5 157.2 159.4 161.9 171 172.8 173.2 174.7 175.4 176.4 1 1 1 1 1 1 1 1 1 1 1 1 1 177 178.2 179.6 180.6 182.8 183.6 188.6 191.8 193.8 196.4 199 200.2 201.4 1 1 1 1 1 1 1 1 1 1 1 1 1 201.7 203.2 203.3 205.5 206.1 207.2 208.5 209.1 213.2 220.5 223.8 225.8 227.3 1 1 1 1 1 1 1 1 2 1 1 1 1 228.4 232.3 233.8 234 235.1 235.6 236.4 240.7 240.8 240.9 241.1 246.7 249 1 1 1 1 1 1 1 1 1 1 1 1 1 253.6 253.8 260.1 261.3 264.6 264.9 266.2 269 271 279.9 280.7 289.3 292.9 2 1 1 1 1 1 1 1 1 1 1 1 1 298.4 299.5 300.9 306.6 309.2 309.5 311.5 316.2 317.9 326.3 327.7 329.7 331.6 1 1 1 1 1 1 1 1 1 1 1 1 1 338.3 343.3 347 347.4 351.6 355.1 356.1 356.9 357 363.5 365.8 366.5 366.9 1 1 1 1 1 1 1 1 1 1 1 1 1 367.9 368.1 369 378.8 384.6 395.6 396.6 398.3 403.7 431.3 438.8 464.8 479.1 1 1 1 1 1 1 1 1 1 1 1 1 1 > colnames(x) [1] "Werkloosheid" "Consumentenvertrouwen" > colnames(x)[par1] [1] "Werkloosheid" > x[,par1] [1] 235.1 280.7 264.6 240.7 201.4 240.8 241.1 223.8 206.1 174.7 203.3 220.5 [13] 299.5 347.4 338.3 327.7 351.6 396.6 438.8 395.6 363.5 378.8 357.0 369.0 [25] 464.8 479.1 431.3 366.5 326.3 355.1 331.6 261.3 249.0 205.5 235.6 240.9 [37] 264.9 253.8 232.3 193.8 177.0 213.2 207.2 180.6 188.6 175.4 199.0 179.6 [49] 225.8 234.0 200.2 183.6 178.2 203.2 208.5 191.8 172.8 148.0 159.4 154.5 [61] 213.2 196.4 182.8 176.4 153.6 173.2 171.0 151.2 161.9 157.2 201.7 236.4 [73] 356.1 398.3 403.7 384.6 365.8 368.1 367.9 347.0 343.3 292.9 311.5 300.9 [85] 366.9 356.9 329.7 316.2 269.0 289.3 266.2 253.6 233.8 228.4 253.6 260.1 [97] 306.6 309.2 309.5 271.0 279.9 317.9 298.4 246.7 227.3 209.1 > 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/1lr0n1291988416.tab") + } + } > m Conditional inference tree with 2 terminal nodes Response: Werkloosheid Input: Consumentenvertrouwen Number of observations: 106 1) Consumentenvertrouwen <= -3; criterion = 1, statistic = 13.902 2)* weights = 75 1) Consumentenvertrouwen > -3 3)* weights = 31 > postscript(file="/var/www/html/freestat/rcomp/tmp/2lr0n1291988416.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/3lr0n1291988416.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 235.1 218.4097 16.690323 2 280.7 218.4097 62.290323 3 264.6 218.4097 46.190323 4 240.7 218.4097 22.290323 5 201.4 218.4097 -17.009677 6 240.8 218.4097 22.390323 7 241.1 218.4097 22.690323 8 223.8 288.1693 -64.369333 9 206.1 288.1693 -82.069333 10 174.7 288.1693 -113.469333 11 203.3 288.1693 -84.869333 12 220.5 288.1693 -67.669333 13 299.5 288.1693 11.330667 14 347.4 288.1693 59.230667 15 338.3 288.1693 50.130667 16 327.7 288.1693 39.530667 17 351.6 288.1693 63.430667 18 396.6 288.1693 108.430667 19 438.8 288.1693 150.630667 20 395.6 288.1693 107.430667 21 363.5 288.1693 75.330667 22 378.8 288.1693 90.630667 23 357.0 288.1693 68.830667 24 369.0 288.1693 80.830667 25 464.8 288.1693 176.630667 26 479.1 218.4097 260.690323 27 431.3 288.1693 143.130667 28 366.5 288.1693 78.330667 29 326.3 288.1693 38.130667 30 355.1 218.4097 136.690323 31 331.6 288.1693 43.430667 32 261.3 218.4097 42.890323 33 249.0 218.4097 30.590323 34 205.5 288.1693 -82.669333 35 235.6 288.1693 -52.569333 36 240.9 288.1693 -47.269333 37 264.9 288.1693 -23.269333 38 253.8 288.1693 -34.369333 39 232.3 218.4097 13.890323 40 193.8 288.1693 -94.369333 41 177.0 288.1693 -111.169333 42 213.2 288.1693 -74.969333 43 207.2 288.1693 -80.969333 44 180.6 288.1693 -107.569333 45 188.6 288.1693 -99.569333 46 175.4 288.1693 -112.769333 47 199.0 288.1693 -89.169333 48 179.6 288.1693 -108.569333 49 225.8 218.4097 7.390323 50 234.0 218.4097 15.590323 51 200.2 288.1693 -87.969333 52 183.6 288.1693 -104.569333 53 178.2 288.1693 -109.969333 54 203.2 218.4097 -15.209677 55 208.5 218.4097 -9.909677 56 191.8 218.4097 -26.609677 57 172.8 218.4097 -45.609677 58 148.0 218.4097 -70.409677 59 159.4 218.4097 -59.009677 60 154.5 288.1693 -133.669333 61 213.2 218.4097 -5.209677 62 196.4 218.4097 -22.009677 63 182.8 218.4097 -35.609677 64 176.4 218.4097 -42.009677 65 153.6 218.4097 -64.809677 66 173.2 218.4097 -45.209677 67 171.0 218.4097 -47.409677 68 151.2 218.4097 -67.209677 69 161.9 218.4097 -56.509677 70 157.2 218.4097 -61.209677 71 201.7 288.1693 -86.469333 72 236.4 288.1693 -51.769333 73 356.1 288.1693 67.930667 74 398.3 288.1693 110.130667 75 403.7 288.1693 115.530667 76 384.6 288.1693 96.430667 77 365.8 288.1693 77.630667 78 368.1 288.1693 79.930667 79 367.9 288.1693 79.730667 80 347.0 288.1693 58.830667 81 343.3 288.1693 55.130667 82 292.9 288.1693 4.730667 83 311.5 288.1693 23.330667 84 300.9 288.1693 12.730667 85 366.9 288.1693 78.730667 86 356.9 288.1693 68.730667 87 329.7 288.1693 41.530667 88 316.2 288.1693 28.030667 89 269.0 288.1693 -19.169333 90 289.3 288.1693 1.130667 91 266.2 288.1693 -21.969333 92 253.6 288.1693 -34.569333 93 233.8 288.1693 -54.369333 94 228.4 288.1693 -59.769333 95 253.6 288.1693 -34.569333 96 260.1 288.1693 -28.069333 97 306.6 288.1693 18.430667 98 309.2 288.1693 21.030667 99 309.5 288.1693 21.330667 100 271.0 288.1693 -17.169333 101 279.9 288.1693 -8.269333 102 317.9 288.1693 29.730667 103 298.4 288.1693 10.230667 104 246.7 288.1693 -41.469333 105 227.3 288.1693 -60.869333 106 209.1 218.4097 -9.309677 > 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/4w0z81291988416.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/5hjyw1291988416.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/621ek1291988416.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/7dsen1291988416.tab") + } > > try(system("convert tmp/2lr0n1291988416.ps tmp/2lr0n1291988416.png",intern=TRUE)) character(0) > try(system("convert tmp/3lr0n1291988416.ps tmp/3lr0n1291988416.png",intern=TRUE)) character(0) > try(system("convert tmp/4w0z81291988416.ps tmp/4w0z81291988416.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.695 0.733 3.857