R version 2.9.0 (2009-04-17) Copyright (C) 2009 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. 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(113,14.3,15.89,110,14.2,16.93,107,15.9,20.28,103,15.3,22.52,98,15.5,23.51,98,15.1,22.59,137,15,23.51,148,12.1,24.76,147,15.8,26.08,139,16.9,25.29,130,15.1,23.38,128,13.7,25.29,127,14.8,28.42,123,14.7,31.85,118,16,30.1,114,15.4,25.45,108,15,24.95,111,15.5,26.84,151,15.1,27.52,159,11.7,27.94,158,16.3,25.23,148,16.7,26.53,138,15,27.21,137,14.9,28.53,136,14.6,30.35,133,15.3,31.21,126,17.9,32.86,120,16.4,33.2,114,15.4,35.73,116,17.9,34.53,153,15.9,36.54,162,13.9,40.1,161,17.8,40.56,149,17.9,46.14,139,17.4,42.85,135,16.7,38.22,130,16,40.18,127,16.6,42.19,122,19.1,47.56,117,17.8,47.26,112,17.2,44.03,113,18.6,49.83,149,16.3,53.35,157,15.1,58.9,157,19.2,59.64,147,17.7,56.99,137,19.1,53.2,132,18,53.24,125,17.5,57.85,123,17.8,55.69,117,21.1,55.64,114,17.2,62.52,111,19.4,64.4,112,19.8,64.65,144,17.6,67.71,150,16.2,67.21,149,19.5,59.37,134,19.9,53.26,123,20,52.42,116,17.3,55.03),dim=c(3,60),dimnames=list(c('werklozen-25','buitenlandse_handel','ruwe_aardolie'),1:60)) > y <- array(NA,dim=c(3,60),dimnames=list(c('werklozen-25','buitenlandse_handel','ruwe_aardolie'),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 = '0' > 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] "werklozen.25" > x[,par1] [1] 113 110 107 103 98 98 137 148 147 139 130 128 127 123 118 114 108 111 151 [20] 159 158 148 138 137 136 133 126 120 114 116 153 162 161 149 139 135 130 127 [39] 122 117 112 113 149 157 157 147 137 132 125 123 117 114 111 112 144 150 149 [58] 134 123 116 > 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]) 98 103 107 108 110 111 112 113 114 116 117 118 120 122 123 125 126 127 128 130 2 1 1 1 1 2 2 2 3 2 2 1 1 1 3 1 1 2 1 2 132 133 134 135 136 137 138 139 144 147 148 149 150 151 153 157 158 159 161 162 1 1 1 1 1 3 1 2 1 2 2 3 1 1 1 2 1 1 1 1 > colnames(x) [1] "werklozen.25" "buitenlandse_handel" "ruwe_aardolie" > colnames(x)[par1] [1] "werklozen.25" > x[,par1] [1] 113 110 107 103 98 98 137 148 147 139 130 128 127 123 118 114 108 111 151 [20] 159 158 148 138 137 136 133 126 120 114 116 153 162 161 149 139 135 130 127 [39] 122 117 112 113 149 157 157 147 137 132 125 123 117 114 111 112 144 150 149 [58] 134 123 116 > if (par2 == 'none') { + m <- ctree(as.formula(paste(colnames(x)[par1],' ~ .',sep='')),data = x) + } > > #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/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/rcomp/tmp/1ncgo1292360606.tab") + } + } > m Conditional inference tree with 1 terminal nodes Response: werklozen.25 Inputs: buitenlandse_handel, ruwe_aardolie Number of observations: 60 1)* weights = 60 > postscript(file="/var/www/html/rcomp/tmp/2ncgo1292360606.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/rcomp/tmp/3ncgo1292360606.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 113 130.2 -17.2 2 110 130.2 -20.2 3 107 130.2 -23.2 4 103 130.2 -27.2 5 98 130.2 -32.2 6 98 130.2 -32.2 7 137 130.2 6.8 8 148 130.2 17.8 9 147 130.2 16.8 10 139 130.2 8.8 11 130 130.2 -0.2 12 128 130.2 -2.2 13 127 130.2 -3.2 14 123 130.2 -7.2 15 118 130.2 -12.2 16 114 130.2 -16.2 17 108 130.2 -22.2 18 111 130.2 -19.2 19 151 130.2 20.8 20 159 130.2 28.8 21 158 130.2 27.8 22 148 130.2 17.8 23 138 130.2 7.8 24 137 130.2 6.8 25 136 130.2 5.8 26 133 130.2 2.8 27 126 130.2 -4.2 28 120 130.2 -10.2 29 114 130.2 -16.2 30 116 130.2 -14.2 31 153 130.2 22.8 32 162 130.2 31.8 33 161 130.2 30.8 34 149 130.2 18.8 35 139 130.2 8.8 36 135 130.2 4.8 37 130 130.2 -0.2 38 127 130.2 -3.2 39 122 130.2 -8.2 40 117 130.2 -13.2 41 112 130.2 -18.2 42 113 130.2 -17.2 43 149 130.2 18.8 44 157 130.2 26.8 45 157 130.2 26.8 46 147 130.2 16.8 47 137 130.2 6.8 48 132 130.2 1.8 49 125 130.2 -5.2 50 123 130.2 -7.2 51 117 130.2 -13.2 52 114 130.2 -16.2 53 111 130.2 -19.2 54 112 130.2 -18.2 55 144 130.2 13.8 56 150 130.2 19.8 57 149 130.2 18.8 58 134 130.2 3.8 59 123 130.2 -7.2 60 116 130.2 -14.2 > 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/rcomp/tmp/4xlf91292360606.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/rcomp/tmp/5tvvz1292360606.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/rcomp/tmp/6fwc51292360606.tab") + } Warning message: In cor(result$Forecasts, result$Actuals) : the standard deviation is zero > 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/rcomp/tmp/7p5b81292360606.tab") + } > > try(system("convert tmp/2ncgo1292360606.ps tmp/2ncgo1292360606.png",intern=TRUE)) character(0) > try(system("convert tmp/3ncgo1292360606.ps tmp/3ncgo1292360606.png",intern=TRUE)) character(0) > try(system("convert tmp/4xlf91292360606.ps tmp/4xlf91292360606.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.920 0.548 5.289