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(3.18 + ,0.22 + ,6.62 + ,3.64 + ,3.14 + ,0.22 + ,6.56 + ,3.62 + ,3.02 + ,0.23 + ,6.59 + ,3.61 + ,3.02 + ,0.24 + ,6.56 + ,3.6 + ,3.03 + ,0.25 + ,6.57 + ,3.6 + ,3.04 + ,0.25 + ,6.62 + ,3.63 + ,3.09 + ,0.24 + ,6.69 + ,3.59 + ,3.06 + ,0.24 + ,6.69 + ,3.55 + ,3.06 + ,0.22 + ,6.64 + ,3.54 + ,3.09 + ,0.21 + ,6.6 + ,3.53 + ,3.11 + ,0.21 + ,6.66 + ,3.53 + ,3.1 + ,0.21 + ,6.62 + ,3.53 + ,3.09 + ,0.2 + ,6.64 + ,3.52 + ,3.19 + ,0.2 + ,6.64 + ,3.52 + ,3.22 + ,0.2 + ,6.73 + ,3.48 + ,3.22 + ,0.2 + ,6.73 + ,3.49 + ,3.25 + ,0.2 + ,6.69 + ,3.47 + ,3.25 + ,0.2 + ,6.78 + ,3.46 + ,3.27 + ,0.2 + ,6.77 + ,3.4 + ,3.28 + ,0.2 + ,6.8 + ,3.36 + ,3.24 + ,0.2 + ,6.8 + ,3.3 + ,3.23 + ,0.2 + ,6.74 + ,3.28 + ,3.2 + ,0.2 + ,6.84 + ,3.28 + ,3.19 + ,0.2 + ,6.83 + ,3.24 + ,3.23 + ,0.2 + ,6.89 + ,3.23 + ,3.19 + ,0.2 + ,6.9 + ,3.2 + ,3.16 + ,0.2 + ,6.86 + ,3.15 + ,3.11 + ,0.2 + ,6.78 + ,3.1 + ,3.11 + ,0.2 + ,6.82 + ,3.07 + ,3.07 + ,0.2 + ,6.81 + ,3.03 + ,3.05 + ,0.21 + ,6.81 + ,2.96 + ,3 + ,0.2 + ,6.78 + ,2.88 + ,2.95 + ,0.2 + ,6.79 + ,2.83 + ,2.9 + ,0.19 + ,6.83 + ,2.8 + ,2.88 + ,0.18 + ,6.9 + ,2.8 + ,2.9 + ,0.18 + ,6.79 + ,2.79 + ,2.89 + ,0.17 + ,6.88 + ,2.79 + ,2.89 + ,0.17 + ,6.89 + ,2.78 + ,2.91 + ,0.17 + ,6.91 + ,2.79 + ,2.9 + ,0.17 + ,6.93 + ,2.78 + ,2.9 + ,0.17 + ,6.89 + ,2.78 + ,2.88 + ,0.16 + ,7 + ,2.74 + ,2.83 + ,0.16 + ,7.01 + ,2.71 + ,2.8 + ,0.16 + ,7.15 + ,2.69 + ,2.77 + ,0.16 + ,7.25 + ,2.68 + ,2.78 + ,0.16 + ,7.33 + ,2.68 + ,2.75 + ,0.16 + ,7.39 + ,2.68 + ,2.74 + ,0.15 + ,7.38 + ,2.69 + ,2.73 + ,0.15 + ,7.38 + ,2.68 + ,2.69 + ,0.15 + ,7.35 + ,2.69 + ,2.67 + ,0.15 + ,7.38 + ,2.68 + ,2.66 + ,0.15 + ,7.34 + ,2.68 + ,2.67 + ,0.16 + ,7.25 + ,2.63 + ,2.65 + ,0.15 + ,7.07 + ,2.58 + ,2.64 + ,0.15 + ,6.73 + ,2.52 + ,2.63 + ,0.15 + ,6.56 + ,2.5) + ,dim=c(4 + ,56) + ,dimnames=list(c('Mayonaise' + ,'Eieren' + ,'Olijfolie' + ,'Mosterd') + ,1:56)) > y <- array(NA,dim=c(4,56),dimnames=list(c('Mayonaise','Eieren','Olijfolie','Mosterd'),1:56)) > 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] "Mayonaise" > x[,par1] [1] 3.18 3.14 3.02 3.02 3.03 3.04 3.09 3.06 3.06 3.09 3.11 3.10 3.09 3.19 3.22 [16] 3.22 3.25 3.25 3.27 3.28 3.24 3.23 3.20 3.19 3.23 3.19 3.16 3.11 3.11 3.07 [31] 3.05 3.00 2.95 2.90 2.88 2.90 2.89 2.89 2.91 2.90 2.90 2.88 2.83 2.80 2.77 [46] 2.78 2.75 2.74 2.73 2.69 2.67 2.66 2.67 2.65 2.64 2.63 > 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.63 2.64 2.65 2.66 2.67 2.69 2.73 2.74 2.75 2.77 2.78 2.8 2.83 2.88 2.89 2.9 1 1 1 1 2 1 1 1 1 1 1 1 1 2 2 4 2.91 2.95 3 3.02 3.03 3.04 3.05 3.06 3.07 3.09 3.1 3.11 3.14 3.16 3.18 3.19 1 1 1 2 1 1 1 2 1 3 1 3 1 1 1 3 3.2 3.22 3.23 3.24 3.25 3.27 3.28 1 2 2 1 2 1 1 > colnames(x) [1] "Mayonaise" "Eieren" "Olijfolie" "Mosterd" > colnames(x)[par1] [1] "Mayonaise" > x[,par1] [1] 3.18 3.14 3.02 3.02 3.03 3.04 3.09 3.06 3.06 3.09 3.11 3.10 3.09 3.19 3.22 [16] 3.22 3.25 3.25 3.27 3.28 3.24 3.23 3.20 3.19 3.23 3.19 3.16 3.11 3.11 3.07 [31] 3.05 3.00 2.95 2.90 2.88 2.90 2.89 2.89 2.91 2.90 2.90 2.88 2.83 2.80 2.77 [46] 2.78 2.75 2.74 2.73 2.69 2.67 2.66 2.67 2.65 2.64 2.63 > 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/1ejwe1292347258.tab") + } + } > m Conditional inference tree with 4 terminal nodes Response: Mayonaise Inputs: Eieren, Olijfolie, Mosterd Number of observations: 56 1) Mosterd <= 2.83; criterion = 1, statistic = 37.696 2) Mosterd <= 2.69; criterion = 1, statistic = 18.947 3)* weights = 13 2) Mosterd > 2.69 4)* weights = 11 1) Mosterd > 2.83 5) Eieren <= 0.2; criterion = 0.999, statistic = 12.174 6)* weights = 19 5) Eieren > 0.2 7)* weights = 13 > postscript(file="/var/www/html/rcomp/tmp/2oadz1292347258.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/3oadz1292347258.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.18 3.076154 0.103846154 2 3.14 3.076154 0.063846154 3 3.02 3.076154 -0.056153846 4 3.02 3.076154 -0.056153846 5 3.03 3.076154 -0.046153846 6 3.04 3.076154 -0.036153846 7 3.09 3.076154 0.013846154 8 3.06 3.076154 -0.016153846 9 3.06 3.076154 -0.016153846 10 3.09 3.076154 0.013846154 11 3.11 3.076154 0.033846154 12 3.10 3.076154 0.023846154 13 3.09 3.184211 -0.094210526 14 3.19 3.184211 0.005789474 15 3.22 3.184211 0.035789474 16 3.22 3.184211 0.035789474 17 3.25 3.184211 0.065789474 18 3.25 3.184211 0.065789474 19 3.27 3.184211 0.085789474 20 3.28 3.184211 0.095789474 21 3.24 3.184211 0.055789474 22 3.23 3.184211 0.045789474 23 3.20 3.184211 0.015789474 24 3.19 3.184211 0.005789474 25 3.23 3.184211 0.045789474 26 3.19 3.184211 0.005789474 27 3.16 3.184211 -0.024210526 28 3.11 3.184211 -0.074210526 29 3.11 3.184211 -0.074210526 30 3.07 3.184211 -0.114210526 31 3.05 3.076154 -0.026153846 32 3.00 3.184211 -0.184210526 33 2.95 2.893636 0.056363636 34 2.90 2.893636 0.006363636 35 2.88 2.893636 -0.013636364 36 2.90 2.893636 0.006363636 37 2.89 2.893636 -0.003636364 38 2.89 2.893636 -0.003636364 39 2.91 2.893636 0.016363636 40 2.90 2.893636 0.006363636 41 2.90 2.893636 0.006363636 42 2.88 2.893636 -0.013636364 43 2.83 2.893636 -0.063636364 44 2.80 2.706154 0.093846154 45 2.77 2.706154 0.063846154 46 2.78 2.706154 0.073846154 47 2.75 2.706154 0.043846154 48 2.74 2.706154 0.033846154 49 2.73 2.706154 0.023846154 50 2.69 2.706154 -0.016153846 51 2.67 2.706154 -0.036153846 52 2.66 2.706154 -0.046153846 53 2.67 2.706154 -0.036153846 54 2.65 2.706154 -0.056153846 55 2.64 2.706154 -0.066153846 56 2.63 2.706154 -0.076153846 > 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/4h1ck1292347258.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/5dbsa1292347258.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/6o2rd1292347258.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/rcomp/tmp/7938j1292347258.tab") + } > > try(system("convert tmp/2oadz1292347258.ps tmp/2oadz1292347258.png",intern=TRUE)) character(0) > try(system("convert tmp/3oadz1292347258.ps tmp/3oadz1292347258.png",intern=TRUE)) character(0) > try(system("convert tmp/4h1ck1292347258.ps tmp/4h1ck1292347258.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.247 0.570 5.964