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(3.48,4143,3.6,4429,3.66,5219,3.45,4929,3.3,5761,3.14,5592,3.21,4163,3.12,4962,3.14,5208,3.4,4755,3.42,4491,3.29,5732,3.49,5731,3.52,5040,3.81,6102,4.03,4904,3.98,5369,4.1,5578,3.96,4619,3.83,4731,3.72,5011,3.82,5299,3.76,4146,3.98,4625,4.14,4736,4,4219,4.13,5116,4.28,4205,4.46,4121,4.63,5103,4.49,4300,4.41,4578,4.5,3809,4.39,5657,4.33,4248,4.45,3830,4.17,4736,4.13,4839,4.33,4411,4.47,4570,4.63,4104,4.9,4801,4.77,3953,4.51,3828,4.63,4440,4.36,4026,3.95,4109,3.74,4785,4.15,3224,4.14,3552,3.97,3940,3.81,3913,4.07,3681,3.84,4309,3.63,3830,3.55,4143,3.6,4087,3.63,3818,3.55,3380,3.69,3430,3.53,3458,3.43,3970,3.4,5260,3.41,5024,3.09,5634,3.35,6549,3.22,4676),dim=c(2,67),dimnames=list(c('leningen','nieuwbouw'),1:67)) > y <- array(NA,dim=c(2,67),dimnames=list(c('leningen','nieuwbouw'),1:67)) > for (i in 1:dim(x)[1]) + { + for (j in 1:dim(x)[2]) + { + y[i,j] <- as.numeric(x[i,j]) + } + } > par4 = 'yes' > 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 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] "leningen" > x[,par1] [1] 3.48 3.60 3.66 3.45 3.30 3.14 3.21 3.12 3.14 3.40 3.42 3.29 3.49 3.52 3.81 [16] 4.03 3.98 4.10 3.96 3.83 3.72 3.82 3.76 3.98 4.14 4.00 4.13 4.28 4.46 4.63 [31] 4.49 4.41 4.50 4.39 4.33 4.45 4.17 4.13 4.33 4.47 4.63 4.90 4.77 4.51 4.63 [46] 4.36 3.95 3.74 4.15 4.14 3.97 3.81 4.07 3.84 3.63 3.55 3.60 3.63 3.55 3.69 [61] 3.53 3.43 3.40 3.41 3.09 3.35 3.22 > 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]) 3.09 3.12 3.14 3.21 3.22 3.29 3.3 3.35 3.4 3.41 3.42 3.43 3.45 3.48 3.49 3.52 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 1 3.53 3.55 3.6 3.63 3.66 3.69 3.72 3.74 3.76 3.81 3.82 3.83 3.84 3.95 3.96 3.97 1 2 2 2 1 1 1 1 1 2 1 1 1 1 1 1 3.98 4 4.03 4.07 4.1 4.13 4.14 4.15 4.17 4.28 4.33 4.36 4.39 4.41 4.45 4.46 2 1 1 1 1 2 2 1 1 1 2 1 1 1 1 1 4.47 4.49 4.5 4.51 4.63 4.77 4.9 1 1 1 1 3 1 1 > colnames(x) [1] "leningen" "nieuwbouw" > colnames(x)[par1] [1] "leningen" > x[,par1] [1] 3.48 3.60 3.66 3.45 3.30 3.14 3.21 3.12 3.14 3.40 3.42 3.29 3.49 3.52 3.81 [16] 4.03 3.98 4.10 3.96 3.83 3.72 3.82 3.76 3.98 4.14 4.00 4.13 4.28 4.46 4.63 [31] 4.49 4.41 4.50 4.39 4.33 4.45 4.17 4.13 4.33 4.47 4.63 4.90 4.77 4.51 4.63 [46] 4.36 3.95 3.74 4.15 4.14 3.97 3.81 4.07 3.84 3.63 3.55 3.60 3.63 3.55 3.69 [61] 3.53 3.43 3.40 3.41 3.09 3.35 3.22 > 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/1e4q91293205885.tab") + } + } > m Conditional inference tree with 2 terminal nodes Response: leningen Input: nieuwbouw Number of observations: 67 1) nieuwbouw <= 4904; criterion = 0.986, statistic = 6.048 2)* weights = 46 1) nieuwbouw > 4904 3)* weights = 21 > postscript(file="/var/www/rcomp/tmp/2e4q91293205885.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/3e4q91293205885.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 3.48 4.002826 -0.522826087 2 3.60 4.002826 -0.402826087 3 3.66 3.616190 0.043809524 4 3.45 3.616190 -0.166190476 5 3.30 3.616190 -0.316190476 6 3.14 3.616190 -0.476190476 7 3.21 4.002826 -0.792826087 8 3.12 3.616190 -0.496190476 9 3.14 3.616190 -0.476190476 10 3.40 4.002826 -0.602826087 11 3.42 4.002826 -0.582826087 12 3.29 3.616190 -0.326190476 13 3.49 3.616190 -0.126190476 14 3.52 3.616190 -0.096190476 15 3.81 3.616190 0.193809524 16 4.03 4.002826 0.027173913 17 3.98 3.616190 0.363809524 18 4.10 3.616190 0.483809524 19 3.96 4.002826 -0.042826087 20 3.83 4.002826 -0.172826087 21 3.72 3.616190 0.103809524 22 3.82 3.616190 0.203809524 23 3.76 4.002826 -0.242826087 24 3.98 4.002826 -0.022826087 25 4.14 4.002826 0.137173913 26 4.00 4.002826 -0.002826087 27 4.13 3.616190 0.513809524 28 4.28 4.002826 0.277173913 29 4.46 4.002826 0.457173913 30 4.63 3.616190 1.013809524 31 4.49 4.002826 0.487173913 32 4.41 4.002826 0.407173913 33 4.50 4.002826 0.497173913 34 4.39 3.616190 0.773809524 35 4.33 4.002826 0.327173913 36 4.45 4.002826 0.447173913 37 4.17 4.002826 0.167173913 38 4.13 4.002826 0.127173913 39 4.33 4.002826 0.327173913 40 4.47 4.002826 0.467173913 41 4.63 4.002826 0.627173913 42 4.90 4.002826 0.897173913 43 4.77 4.002826 0.767173913 44 4.51 4.002826 0.507173913 45 4.63 4.002826 0.627173913 46 4.36 4.002826 0.357173913 47 3.95 4.002826 -0.052826087 48 3.74 4.002826 -0.262826087 49 4.15 4.002826 0.147173913 50 4.14 4.002826 0.137173913 51 3.97 4.002826 -0.032826087 52 3.81 4.002826 -0.192826087 53 4.07 4.002826 0.067173913 54 3.84 4.002826 -0.162826087 55 3.63 4.002826 -0.372826087 56 3.55 4.002826 -0.452826087 57 3.60 4.002826 -0.402826087 58 3.63 4.002826 -0.372826087 59 3.55 4.002826 -0.452826087 60 3.69 4.002826 -0.312826087 61 3.53 4.002826 -0.472826087 62 3.43 4.002826 -0.572826087 63 3.40 3.616190 -0.216190476 64 3.41 3.616190 -0.206190476 65 3.09 3.616190 -0.526190476 66 3.35 3.616190 -0.266190476 67 3.22 4.002826 -0.782826087 > 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/47w7c1293205885.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/5l6nl1293205885.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/6wf461293205885.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/7hf2c1293205885.tab") + } > > try(system("convert tmp/2e4q91293205885.ps tmp/2e4q91293205885.png",intern=TRUE)) character(0) > try(system("convert tmp/3e4q91293205885.ps tmp/3e4q91293205885.png",intern=TRUE)) character(0) > try(system("convert tmp/47w7c1293205885.ps tmp/47w7c1293205885.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.030 0.370 2.394