R version 2.12.0 (2010-10-15) Copyright (C) 2010 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-pc-linux-gnu (32-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > x <- array(list(-5,-6,33,5,15,-1,-3,24,6,17,-2,-4,24,6,13,-5,-7,31,5,12,-4,-7,25,5,13,-6,-7,28,3,10,-2,-3,24,5,14,-2,0,25,5,13,-2,-5,16,5,10,-2,-3,17,3,11,2,3,11,6,12,1,2,12,6,7,-8,-7,39,4,11,-1,-1,19,6,9,1,0,14,5,13,-1,-3,15,4,12,2,4,7,5,5,2,2,12,5,13,1,3,12,4,11,-1,0,14,3,8,-2,-10,9,2,8,-2,-10,8,3,8,-1,-9,4,2,8,-8,-22,7,-1,0,-4,-16,3,0,3,-6,-18,5,-2,0,-3,-14,0,1,-1,-3,-12,-2,-2,-1,-7,-17,6,-2,-4,-9,-23,11,-2,1,-11,-28,9,-6,-1,-13,-31,17,-4,0,-11,-21,21,-2,-1,-9,-19,21,0,6,-17,-22,41,-5,0,-22,-22,57,-4,-3,-25,-25,65,-5,-3,-20,-16,68,-1,4,-24,-22,73,-2,1,-24,-21,71,-4,0,-22,-10,71,-1,-4,-19,-7,70,1,-2,-18,-5,69,1,3,-17,-4,65,-2,2,-11,7,57,1,5,-11,6,57,1,6,-12,3,57,3,6,-10,10,55,3,3,-15,0,65,1,4,-15,-2,65,1,7,-15,-1,64,0,5,-13,2,60,2,6,-8,8,43,2,1,-13,-6,47,-1,3,-9,-4,40,1,6,-7,4,31,0,0,-4,7,27,1,3,-4,3,24,1,4,-2,3,23,3,7,0,8,17,2,6),dim=c(5,60),dimnames=list(c('indicator','vooruitzichten','werkloosheid','financiën','spaarvermogen'),1:60)) > y <- array(NA,dim=c(5,60),dimnames=list(c('indicator','vooruitzichten','werkloosheid','financiën','spaarvermogen'),1:60)) > 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 = '4' > par2 = 'none' > par1 = '3' > #'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 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] 33 24 24 31 25 28 24 25 16 17 11 12 39 19 14 15 7 12 12 14 9 8 4 7 3 [26] 5 0 -2 6 11 9 17 21 21 41 57 65 68 73 71 71 70 69 65 57 57 57 55 65 65 [51] 64 60 43 47 40 31 27 24 23 17 > 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]) -2 0 3 4 5 6 7 8 9 11 12 14 15 16 17 19 21 23 24 25 27 28 31 33 39 40 1 1 1 1 1 1 2 1 2 2 3 2 1 1 3 1 2 1 4 2 1 1 2 1 1 1 41 43 47 55 57 60 64 65 68 69 70 71 73 1 1 1 1 4 1 1 4 1 1 1 2 1 > colnames(x) [1] "indicator" "vooruitzichten" "werkloosheid" "financi..n" [5] "spaarvermogen" > colnames(x)[par1] [1] "werkloosheid" > x[,par1] [1] 33 24 24 31 25 28 24 25 16 17 11 12 39 19 14 15 7 12 12 14 9 8 4 7 3 [26] 5 0 -2 6 11 9 17 21 21 41 57 65 68 73 71 71 70 69 65 57 57 57 55 65 65 [51] 64 60 43 47 40 31 27 24 23 17 > if (par2 == 'none') { + m <- ctree(as.formula(paste(colnames(x)[par1],' ~ .',sep='')),data = x) + } > > #Note: the /var/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/rcomp/tmp/14z5x1293271081.tab") + } + } > m Conditional inference tree with 6 terminal nodes Response: werkloosheid Inputs: indicator, vooruitzichten, financi..n, spaarvermogen Number of observations: 60 1) indicator <= -10; criterion = 1, statistic = 41.186 2) indicator <= -15; criterion = 0.965, statistic = 6.841 3)* weights = 13 2) indicator > -15 4)* weights = 9 1) indicator > -10 5) vooruitzichten <= -9; criterion = 0.96, statistic = 6.607 6)* weights = 11 5) vooruitzichten > -9 7) indicator <= -5; criterion = 1, statistic = 23.255 8)* weights = 7 7) indicator > -5 9) indicator <= -1; criterion = 0.999, statistic = 14.264 10)* weights = 13 9) indicator > -1 11)* weights = 7 > postscript(file="/var/www/rcomp/tmp/2x9m01293271081.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/rcomp/tmp/3x9m01293271081.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 33 35.000000 -2.00000000 2 24 21.307692 2.69230769 3 24 21.307692 2.69230769 4 31 35.000000 -4.00000000 5 25 21.307692 3.69230769 6 28 35.000000 -7.00000000 7 24 21.307692 2.69230769 8 25 21.307692 3.69230769 9 16 21.307692 -5.30769231 10 17 21.307692 -4.30769231 11 11 12.142857 -1.14285714 12 12 12.142857 -0.14285714 13 39 35.000000 4.00000000 14 19 21.307692 -2.30769231 15 14 12.142857 1.85714286 16 15 21.307692 -6.30769231 17 7 12.142857 -5.14285714 18 12 12.142857 -0.14285714 19 12 12.142857 -0.14285714 20 14 21.307692 -7.30769231 21 9 6.545455 2.45454545 22 8 6.545455 1.45454545 23 4 6.545455 -2.54545455 24 7 6.545455 0.45454545 25 3 6.545455 -3.54545455 26 5 6.545455 -1.54545455 27 0 6.545455 -6.54545455 28 -2 6.545455 -8.54545455 29 6 6.545455 -0.54545455 30 11 6.545455 4.45454545 31 9 42.222222 -33.22222222 32 17 42.222222 -25.22222222 33 21 42.222222 -21.22222222 34 21 6.545455 14.45454545 35 41 64.923077 -23.92307692 36 57 64.923077 -7.92307692 37 65 64.923077 0.07692308 38 68 64.923077 3.07692308 39 73 64.923077 8.07692308 40 71 64.923077 6.07692308 41 71 64.923077 6.07692308 42 70 64.923077 5.07692308 43 69 64.923077 4.07692308 44 65 64.923077 0.07692308 45 57 42.222222 14.77777778 46 57 42.222222 14.77777778 47 57 42.222222 14.77777778 48 55 42.222222 12.77777778 49 65 64.923077 0.07692308 50 65 64.923077 0.07692308 51 64 64.923077 -0.92307692 52 60 42.222222 17.77777778 53 43 35.000000 8.00000000 54 47 42.222222 4.77777778 55 40 35.000000 5.00000000 56 31 35.000000 -4.00000000 57 27 21.307692 5.69230769 58 24 21.307692 2.69230769 59 23 21.307692 1.69230769 60 17 12.142857 4.85714286 > 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/rcomp/tmp/47hll1293271081.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/rcomp/tmp/5bi2r1293271081.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/rcomp/tmp/6e0iw1293271081.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/rcomp/tmp/7zjh21293271081.tab") + } > > try(system("convert tmp/2x9m01293271081.ps tmp/2x9m01293271081.png",intern=TRUE)) character(0) > try(system("convert tmp/3x9m01293271081.ps tmp/3x9m01293271081.png",intern=TRUE)) character(0) > try(system("convert tmp/47hll1293271081.ps tmp/47hll1293271081.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.230 0.410 2.639