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(6.3,0,3,2.1,3.406028945,4,9.1,1.02325246,4,15.8,-1.698970004,1,5.2,2.204119983,4,10.9,0.51851394,1,8.3,1.717337583,1,11,-0.366531544,4,3.2,2.667452953,5,6.3,-1.096910013,1,6.6,-0.102372909,2,9.5,-0.698970004,2,3.3,1.441852176,5,11,-0.920818754,2,4.7,1.929418926,1,10.4,-1,3,7.4,0.017033339,4,2.1,2.716837723,5,17.9,-2,1,6.1,1.792391689,1,11.9,-1.698970004,3,13.8,0.230448921,1,14.3,0.544068044,1,15.2,-0.318758763,2,10,1,4,6.5,0.209515015,4,7.5,2.283301229,5,10.6,0.397940009,3,7.4,-0.552841969,1,8.4,0.627365857,2,5.7,0.832508913,2,4.9,-0.124938737,3,3.2,0.556302501,5,11,1.744292983,2,4.9,-0.045757491,3,13.2,0.301029996,2,9.7,-1,4,12.8,0.622214023,1,11.9,0.544068044,2),dim=c(3,39),dimnames=list(c('SWS','logWb','ODI '),1:39)) > y <- array(NA,dim=c(3,39),dimnames=list(c('SWS','logWb','ODI '),1:39)) > 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] "SWS" > x[,par1] [1] 6.3 2.1 9.1 15.8 5.2 10.9 8.3 11.0 3.2 6.3 6.6 9.5 3.3 11.0 4.7 [16] 10.4 7.4 2.1 17.9 6.1 11.9 13.8 14.3 15.2 10.0 6.5 7.5 10.6 7.4 8.4 [31] 5.7 4.9 3.2 11.0 4.9 13.2 9.7 12.8 11.9 > 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.1 3.2 3.3 4.7 4.9 5.2 5.7 6.1 6.3 6.5 6.6 7.4 7.5 8.3 8.4 9.1 2 2 1 1 2 1 1 1 2 1 1 2 1 1 1 1 9.5 9.7 10 10.4 10.6 10.9 11 11.9 12.8 13.2 13.8 14.3 15.2 15.8 17.9 1 1 1 1 1 1 3 2 1 1 1 1 1 1 1 > colnames(x) [1] "SWS" "logWb" "ODI." > colnames(x)[par1] [1] "SWS" > x[,par1] [1] 6.3 2.1 9.1 15.8 5.2 10.9 8.3 11.0 3.2 6.3 6.6 9.5 3.3 11.0 4.7 [16] 10.4 7.4 2.1 17.9 6.1 11.9 13.8 14.3 15.2 10.0 6.5 7.5 10.6 7.4 8.4 [31] 5.7 4.9 3.2 11.0 4.9 13.2 9.7 12.8 11.9 > 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/134y21292848381.tab") + } + } > m Conditional inference tree with 3 terminal nodes Response: SWS Inputs: logWb, ODI. Number of observations: 39 1) logWb <= 1.744293; criterion = 1, statistic = 13.426 2) ODI. <= 2; criterion = 0.994, statistic = 8.758 3)* weights = 18 2) ODI. > 2 4)* weights = 14 1) logWb > 1.744293 5)* weights = 7 > postscript(file="/var/www/html/rcomp/tmp/234y21292848381.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/334y21292848381.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 6.3 7.800000 -1.5000000 2 2.1 4.414286 -2.3142857 3 9.1 7.800000 1.3000000 4 15.8 11.111111 4.6888889 5 5.2 4.414286 0.7857143 6 10.9 11.111111 -0.2111111 7 8.3 11.111111 -2.8111111 8 11.0 7.800000 3.2000000 9 3.2 4.414286 -1.2142857 10 6.3 11.111111 -4.8111111 11 6.6 11.111111 -4.5111111 12 9.5 11.111111 -1.6111111 13 3.3 7.800000 -4.5000000 14 11.0 11.111111 -0.1111111 15 4.7 4.414286 0.2857143 16 10.4 7.800000 2.6000000 17 7.4 7.800000 -0.4000000 18 2.1 4.414286 -2.3142857 19 17.9 11.111111 6.7888889 20 6.1 4.414286 1.6857143 21 11.9 7.800000 4.1000000 22 13.8 11.111111 2.6888889 23 14.3 11.111111 3.1888889 24 15.2 11.111111 4.0888889 25 10.0 7.800000 2.2000000 26 6.5 7.800000 -1.3000000 27 7.5 4.414286 3.0857143 28 10.6 7.800000 2.8000000 29 7.4 11.111111 -3.7111111 30 8.4 11.111111 -2.7111111 31 5.7 11.111111 -5.4111111 32 4.9 7.800000 -2.9000000 33 3.2 7.800000 -4.6000000 34 11.0 11.111111 -0.1111111 35 4.9 7.800000 -2.9000000 36 13.2 11.111111 2.0888889 37 9.7 7.800000 1.9000000 38 12.8 11.111111 1.6888889 39 11.9 11.111111 0.7888889 > 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/4hxi21292848382.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/5kxy81292848382.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/6nyfe1292848382.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/7ypez1292848382.tab") + } > > try(system("convert tmp/234y21292848381.ps tmp/234y21292848381.png",intern=TRUE)) character(0) > try(system("convert tmp/334y21292848381.ps tmp/334y21292848381.png",intern=TRUE)) character(0) > try(system("convert tmp/4hxi21292848382.ps tmp/4hxi21292848382.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.066 0.573 6.165