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(1.8 + ,0.8 + ,2.9 + ,1.8 + ,2.3 + ,0.8 + ,2.6 + ,1.7 + ,-0.1 + ,2.9 + ,1.7 + ,2.2 + ,1 + ,2.2 + ,1.4 + ,-1.5 + ,2.9 + ,1.6 + ,2.1 + ,0.6 + ,2.3 + ,1.2 + ,-4.4 + ,1.4 + ,1.8 + ,2.4 + ,0.9 + ,2.4 + ,1 + ,-4.2 + ,1.1 + ,1.6 + ,2.5 + ,0.6 + ,2.1 + ,1.7 + ,3.5 + ,1.9 + ,1.5 + ,2.4 + ,0.6 + ,1.9 + ,2.4 + ,10 + ,2.8 + ,1.5 + ,2.3 + ,0.4 + ,2.2 + ,2 + ,8.6 + ,1.4 + ,1.3 + ,2.1 + ,0.3 + ,1.9 + ,2.1 + ,9.5 + ,0.7 + ,1.4 + ,2.3 + ,0 + ,2.3 + ,2 + ,9.9 + ,-0.8 + ,1.4 + ,2.2 + ,0.3 + ,2.1 + ,1.8 + ,10.4 + ,-3.1 + ,1.3 + ,2.1 + ,0.1 + ,2.2 + ,2.7 + ,16 + ,0.1 + ,1.3 + ,2 + ,0 + ,2.3 + ,2.3 + ,12.7 + ,1 + ,1.2 + ,2.1 + ,0 + ,1.9 + ,1.9 + ,10.2 + ,1.9 + ,1.1 + ,2.1 + ,0 + ,1.7 + ,2 + ,8.9 + ,-0.5 + ,1.4 + ,2.5 + ,-0.2 + ,2.5 + ,2.3 + ,12.6 + ,1.5 + ,1.2 + ,2.2 + ,-0.3 + ,2.1 + ,2.8 + ,13.6 + ,3.9 + ,1.5 + ,2.3 + ,0.1 + ,2.4 + ,2.4 + ,14.8 + ,1.9 + ,1.1 + ,2.3 + ,0.1 + ,1.5 + ,2.3 + ,9.5 + ,2.6 + ,1.3 + ,2.2 + ,0.4 + ,1.9 + ,2.7 + ,13.7 + ,1.7 + ,1.5 + ,2.2 + ,0.4 + ,2.1 + ,2.7 + ,17 + ,1.4 + ,1.1 + ,1.6 + ,-0.5 + ,2.2 + ,2.9 + ,14.7 + ,2.8 + ,1.4 + ,1.8 + ,0.5 + ,2 + ,3 + ,17.4 + ,0.5 + ,1.3 + ,1.7 + ,0.4 + ,2 + ,2.2 + ,9 + ,1 + ,1.5 + ,1.9 + ,0.7 + ,2.2 + ,2.3 + ,9.1 + ,1.5 + ,1.6 + ,1.8 + ,0.8 + ,2.3 + ,2.8 + ,12.2 + ,1.8 + ,1.7 + ,1.9 + ,0.8 + ,2.3 + ,2.8 + ,15.9 + ,2.7 + ,1.1 + ,1.5 + ,0 + ,2 + ,2.8 + ,12.9 + ,3 + ,1.6 + ,1 + ,1.1 + ,2.2 + ,2.2 + ,10.9 + ,-0.3 + ,1.3 + ,0.8 + ,0.9 + ,1.9 + ,2.6 + ,10.6 + ,1.1 + ,1.7 + ,1.1 + ,1.1 + ,2.3 + ,2.8 + ,13.2 + ,1.7 + ,1.6 + ,1.5 + ,1 + ,2.2 + ,2.5 + ,9.6 + ,1.6 + ,1.7 + ,1.7 + ,1.1 + ,2.3 + ,2.4 + ,6.4 + ,3 + ,1.9 + ,2.3 + ,1.5 + ,2.1 + ,2.3 + ,5.8 + ,3.3 + ,1.8 + ,2.4 + ,1 + ,2.4 + ,1.9 + ,-1 + ,6.7 + ,1.9 + ,3 + ,1 + ,2.3 + ,1.7 + ,-0.2 + ,5.6 + ,1.6 + ,3 + ,0.9 + ,1.9 + ,2 + ,2.7 + ,6 + ,1.5 + ,3.2 + ,0.8 + ,1.6 + ,2.1 + ,3.6 + ,4.8 + ,1.6 + ,3.2 + ,0.8 + ,1.8 + ,1.7 + ,-0.9 + ,5.9 + ,1.6 + ,3.2 + ,0.8 + ,1.8 + ,1.8 + ,0.3 + ,4.3 + ,1.7 + ,3.5 + ,0.8 + ,2 + ,1.8 + ,-1.1 + ,3.7 + ,2 + ,4 + ,0.9 + ,2.3 + ,1.8 + ,-2.5 + ,5.6 + ,2 + ,4.3 + ,0.8 + ,2.2 + ,1.3 + ,-3.4 + ,1.7 + ,1.9 + ,4.1 + ,0.7 + ,2.2 + ,1.3 + ,-3.5 + ,3.2 + ,1.7 + ,4 + ,0.6 + ,2 + ,1.3 + ,-3.9 + ,3.6 + ,1.8 + ,4.1 + ,0.6 + ,2 + ,1.2 + ,-4.6 + ,1.7 + ,1.9 + ,4.2 + ,1 + ,1.9 + ,1.4 + ,-0.1 + ,0.5 + ,1.7 + ,4.5 + ,1 + ,1.5 + ,2.2 + ,4.3 + ,2.1 + ,2 + ,5.6 + ,1 + ,1.6 + ,2.9 + ,10.2 + ,1.5 + ,2.1 + ,6.5 + ,1.1 + ,1.5 + ,3.1 + ,8.7 + ,2.7 + ,2.4 + ,7.6 + ,1.1 + ,2 + ,3.5 + ,13.3 + ,1.4 + ,2.5 + ,8.5 + ,1.4 + ,1.5 + ,3.6 + ,15 + ,1.2 + ,2.5 + ,8.7 + ,1.2 + ,1.5 + ,4.4 + ,20.7 + ,2.3 + ,2.6 + ,8.3 + ,1.2 + ,1.9 + ,4.1 + ,20.7 + ,1.6 + ,2.2 + ,8.3 + ,1.3 + ,1.1 + ,5.1 + ,26.4 + ,4.7 + ,2.5 + ,8.5 + ,1.4 + ,1.5 + ,5.8 + ,31.2 + ,3.5 + ,2.8 + ,8.7 + ,1.4 + ,2.1 + ,5.9 + ,31.4 + ,4.4 + ,2.8 + ,8.7 + ,1.1 + ,2.3 + ,5.4 + ,26.6 + ,3.9 + ,2.9 + ,8.5 + ,1.1 + ,2.6 + ,5.5 + ,26.6 + ,3.5 + ,3 + ,7.9 + ,1.3 + ,2.9 + ,4.8 + ,19.2 + ,3 + ,3.1 + ,7 + ,1.5 + ,3.2 + ,3.2 + ,6.5 + ,1.6 + ,2.9 + ,5.8 + ,1.5 + ,3.2) + ,dim=c(7 + ,61) + ,dimnames=list(c('HICP' + ,'ED' + ,'NBL' + ,'IT' + ,'BL' + ,'NEI' + ,'D') + ,1:61)) > y <- array(NA,dim=c(7,61),dimnames=list(c('HICP','ED','NBL','IT','BL','NEI','D'),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 = '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 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] "HICP" > x[,par1] [1] 1.8 1.7 1.4 1.2 1.0 1.7 2.4 2.0 2.1 2.0 1.8 2.7 2.3 1.9 2.0 2.3 2.8 2.4 2.3 [20] 2.7 2.7 2.9 3.0 2.2 2.3 2.8 2.8 2.8 2.2 2.6 2.8 2.5 2.4 2.3 1.9 1.7 2.0 2.1 [39] 1.7 1.8 1.8 1.8 1.3 1.3 1.3 1.2 1.4 2.2 2.9 3.1 3.5 3.6 4.4 4.1 5.1 5.8 5.9 [58] 5.4 5.5 4.8 3.2 > 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]) 1 1.2 1.3 1.4 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3 3.1 3.2 1 2 3 2 4 5 2 4 2 3 5 3 1 1 3 5 2 1 1 1 3.5 3.6 4.1 4.4 4.8 5.1 5.4 5.5 5.8 5.9 1 1 1 1 1 1 1 1 1 1 > colnames(x) [1] "HICP" "ED" "NBL" "IT" "BL" "NEI" "D" > colnames(x)[par1] [1] "HICP" > x[,par1] [1] 1.8 1.7 1.4 1.2 1.0 1.7 2.4 2.0 2.1 2.0 1.8 2.7 2.3 1.9 2.0 2.3 2.8 2.4 2.3 [20] 2.7 2.7 2.9 3.0 2.2 2.3 2.8 2.8 2.8 2.2 2.6 2.8 2.5 2.4 2.3 1.9 1.7 2.0 2.1 [39] 1.7 1.8 1.8 1.8 1.3 1.3 1.3 1.2 1.4 2.2 2.9 3.1 3.5 3.6 4.4 4.1 5.1 5.8 5.9 [58] 5.4 5.5 4.8 3.2 > 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/153921293188589.tab") + } + } > m Conditional inference tree with 5 terminal nodes Response: HICP Inputs: ED, NBL, IT, BL, NEI, D Number of observations: 61 1) ED <= 17.4; criterion = 1, statistic = 49.489 2) ED <= 3.5; criterion = 1, statistic = 36.147 3)* weights = 18 2) ED > 3.5 4) IT <= 1.8; criterion = 0.998, statistic = 13.084 5) ED <= 10.9; criterion = 0.999, statistic = 13.44 6)* weights = 15 5) ED > 10.9 7)* weights = 13 4) IT > 1.8 8)* weights = 7 1) ED > 17.4 9)* weights = 8 > postscript(file="/var/www/rcomp/tmp/2gu8n1293188589.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/3gu8n1293188589.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 1.8 1.555556 0.244444444 2 1.7 1.555556 0.144444444 3 1.4 1.555556 -0.155555556 4 1.2 1.555556 -0.355555556 5 1.0 1.555556 -0.555555556 6 1.7 1.555556 0.144444444 7 2.4 2.180000 0.220000000 8 2.0 2.180000 -0.180000000 9 2.1 2.180000 -0.080000000 10 2.0 2.180000 -0.180000000 11 1.8 2.180000 -0.380000000 12 2.7 2.692308 0.007692308 13 2.3 2.692308 -0.392307692 14 1.9 2.180000 -0.280000000 15 2.0 2.180000 -0.180000000 16 2.3 2.692308 -0.392307692 17 2.8 2.692308 0.107692308 18 2.4 2.692308 -0.292307692 19 2.3 2.180000 0.120000000 20 2.7 2.692308 0.007692308 21 2.7 2.692308 0.007692308 22 2.9 2.692308 0.207692308 23 3.0 2.692308 0.307692308 24 2.2 2.180000 0.020000000 25 2.3 2.180000 0.120000000 26 2.8 2.692308 0.107692308 27 2.8 2.692308 0.107692308 28 2.8 2.692308 0.107692308 29 2.2 2.180000 0.020000000 30 2.6 2.180000 0.420000000 31 2.8 2.692308 0.107692308 32 2.5 2.180000 0.320000000 33 2.4 2.985714 -0.585714286 34 2.3 2.180000 0.120000000 35 1.9 1.555556 0.344444444 36 1.7 1.555556 0.144444444 37 2.0 1.555556 0.444444444 38 2.1 2.180000 -0.080000000 39 1.7 1.555556 0.144444444 40 1.8 1.555556 0.244444444 41 1.8 1.555556 0.244444444 42 1.8 1.555556 0.244444444 43 1.3 1.555556 -0.255555556 44 1.3 1.555556 -0.255555556 45 1.3 1.555556 -0.255555556 46 1.2 1.555556 -0.355555556 47 1.4 1.555556 -0.155555556 48 2.2 2.985714 -0.785714286 49 2.9 2.985714 -0.085714286 50 3.1 2.985714 0.114285714 51 3.5 2.985714 0.514285714 52 3.6 2.985714 0.614285714 53 4.4 5.125000 -0.725000000 54 4.1 5.125000 -1.025000000 55 5.1 5.125000 -0.025000000 56 5.8 5.125000 0.675000000 57 5.9 5.125000 0.775000000 58 5.4 5.125000 0.275000000 59 5.5 5.125000 0.375000000 60 4.8 5.125000 -0.325000000 61 3.2 2.985714 0.214285714 > 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/48lq81293188589.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/5c46w1293188589.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/6xmn11293188589.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/7j5371293188589.tab") + } > > try(system("convert tmp/2gu8n1293188589.ps tmp/2gu8n1293188589.png",intern=TRUE)) character(0) > try(system("convert tmp/3gu8n1293188589.ps tmp/3gu8n1293188589.png",intern=TRUE)) character(0) > try(system("convert tmp/48lq81293188589.ps tmp/48lq81293188589.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.300 0.640 2.917