R version 2.8.0 (2008-10-20) Copyright (C) 2008 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. Natural language support but running in an English locale 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(99,94.6,106.3,95.9,128.9,104.7,111.1,102.8,102.9,98.1,130,113.9,87,80.9,87.5,95.7,117.6,113.2,103.4,105.9,110.8,108.8,112.6,102.3,102.5,99,112.4,100.7,135.6,115.5,105.1,100.7,127.7,109.9,137,114.6,91,85.4,90.5,100.5,122.4,114.8,123.3,116.5,124.3,112.9,120,102,118.1,106,119,105.3,142.7,118.8,123.6,106.1,129.6,109.3,151.6,117.2,110.4,92.5,99.2,104.2,130.5,112.5,136.2,122.4,129.7,113.3,128,100,121.6,110.7,135.8,112.8,143.8,109.8,147.5,117.3,136.2,109.1,156.6,115.9,123.3,96,104.5,99.8,139.8,116.8,136.5,115.7,112.1,99.4,118.5,94.3,94.4,91,102.3,93.2,111.4,103.1,99.2,94.1,87.8,91.8,115.8,102.7,79.7,82.6,72.7,89.1,104.5,104.5,103,105.1,95.1,95.1,104.2,88.7,78.3,86.3),dim=c(2,61),dimnames=list(c('tip','tipc '),1:61)) > y <- array(NA,dim=c(2,61),dimnames=list(c('tip','tipc '),1:61)) > 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 = '' > 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] "tip" > x[,par1] [1] 99.0 106.3 128.9 111.1 102.9 130.0 87.0 87.5 117.6 103.4 110.8 112.6 [13] 102.5 112.4 135.6 105.1 127.7 137.0 91.0 90.5 122.4 123.3 124.3 120.0 [25] 118.1 119.0 142.7 123.6 129.6 151.6 110.4 99.2 130.5 136.2 129.7 128.0 [37] 121.6 135.8 143.8 147.5 136.2 156.6 123.3 104.5 139.8 136.5 112.1 118.5 [49] 94.4 102.3 111.4 99.2 87.8 115.8 79.7 72.7 104.5 103.0 95.1 104.2 [61] 78.3 > 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]) 72.7 78.3 79.7 87 87.5 87.8 90.5 91 94.4 95.1 99 99.2 102.3 1 1 1 1 1 1 1 1 1 1 1 2 1 102.5 102.9 103 103.4 104.2 104.5 105.1 106.3 110.4 110.8 111.1 111.4 112.1 1 1 1 1 1 2 1 1 1 1 1 1 1 112.4 112.6 115.8 117.6 118.1 118.5 119 120 121.6 122.4 123.3 123.6 124.3 1 1 1 1 1 1 1 1 1 1 2 1 1 127.7 128 128.9 129.6 129.7 130 130.5 135.6 135.8 136.2 136.5 137 139.8 1 1 1 1 1 1 1 1 1 2 1 1 1 142.7 143.8 147.5 151.6 156.6 1 1 1 1 1 > colnames(x) [1] "tip" "tipc." > colnames(x)[par1] [1] "tip" > x[,par1] [1] 99.0 106.3 128.9 111.1 102.9 130.0 87.0 87.5 117.6 103.4 110.8 112.6 [13] 102.5 112.4 135.6 105.1 127.7 137.0 91.0 90.5 122.4 123.3 124.3 120.0 [25] 118.1 119.0 142.7 123.6 129.6 151.6 110.4 99.2 130.5 136.2 129.7 128.0 [37] 121.6 135.8 143.8 147.5 136.2 156.6 123.3 104.5 139.8 136.5 112.1 118.5 [49] 94.4 102.3 111.4 99.2 87.8 115.8 79.7 72.7 104.5 103.0 95.1 104.2 [61] 78.3 > if (par2 == 'none') { + m <- ctree(as.formula(paste(colnames(x)[par1],' ~ .',sep='')),data = x) + } > > #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/157601293577352.tab") + } + } > m Conditional inference tree with 3 terminal nodes Response: tip Input: tipc. Number of observations: 61 1) tipc. <= 108.8; criterion = 1, statistic = 43.066 2) tipc. <= 91.8; criterion = 1, statistic = 16.929 3)* weights = 8 2) tipc. > 91.8 4)* weights = 31 1) tipc. > 108.8 5)* weights = 22 > postscript(file="/var/www/html/freestat/rcomp/tmp/2yh631293577352.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/freestat/rcomp/tmp/3yh631293577352.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 99.0 109.0645 -10.064516 2 106.3 109.0645 -2.764516 3 128.9 109.0645 19.835484 4 111.1 109.0645 2.035484 5 102.9 109.0645 -6.164516 6 130.0 134.3636 -4.363636 7 87.0 86.8875 0.112500 8 87.5 109.0645 -21.564516 9 117.6 134.3636 -16.763636 10 103.4 109.0645 -5.664516 11 110.8 109.0645 1.735484 12 112.6 109.0645 3.535484 13 102.5 109.0645 -6.564516 14 112.4 109.0645 3.335484 15 135.6 134.3636 1.236364 16 105.1 109.0645 -3.964516 17 127.7 134.3636 -6.663636 18 137.0 134.3636 2.636364 19 91.0 86.8875 4.112500 20 90.5 109.0645 -18.564516 21 122.4 134.3636 -11.963636 22 123.3 134.3636 -11.063636 23 124.3 134.3636 -10.063636 24 120.0 109.0645 10.935484 25 118.1 109.0645 9.035484 26 119.0 109.0645 9.935484 27 142.7 134.3636 8.336364 28 123.6 109.0645 14.535484 29 129.6 134.3636 -4.763636 30 151.6 134.3636 17.236364 31 110.4 109.0645 1.335484 32 99.2 109.0645 -9.864516 33 130.5 134.3636 -3.863636 34 136.2 134.3636 1.836364 35 129.7 134.3636 -4.663636 36 128.0 109.0645 18.935484 37 121.6 134.3636 -12.763636 38 135.8 134.3636 1.436364 39 143.8 134.3636 9.436364 40 147.5 134.3636 13.136364 41 136.2 134.3636 1.836364 42 156.6 134.3636 22.236364 43 123.3 109.0645 14.235484 44 104.5 109.0645 -4.564516 45 139.8 134.3636 5.436364 46 136.5 134.3636 2.136364 47 112.1 109.0645 3.035484 48 118.5 109.0645 9.435484 49 94.4 86.8875 7.512500 50 102.3 109.0645 -6.764516 51 111.4 109.0645 2.335484 52 99.2 109.0645 -9.864516 53 87.8 86.8875 0.912500 54 115.8 109.0645 6.735484 55 79.7 86.8875 -7.187500 56 72.7 86.8875 -14.187500 57 104.5 109.0645 -4.564516 58 103.0 109.0645 -6.064516 59 95.1 109.0645 -13.964516 60 104.2 86.8875 17.312500 61 78.3 86.8875 -8.587500 > 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/freestat/rcomp/tmp/4885o1293577352.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/freestat/rcomp/tmp/5tq3u1293577352.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/freestat/rcomp/tmp/6fr2i1293577352.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/freestat/rcomp/tmp/7is4g1293577353.tab") + } > > try(system("convert tmp/2yh631293577352.ps tmp/2yh631293577352.png",intern=TRUE)) character(0) > try(system("convert tmp/3yh631293577352.ps tmp/3yh631293577352.png",intern=TRUE)) character(0) > try(system("convert tmp/4885o1293577352.ps tmp/4885o1293577352.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.576 0.723 3.729