R version 2.13.0 (2011-04-13) Copyright (C) 2011 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(210907 + ,79 + ,30 + ,2 + ,94 + ,112285 + ,179321 + ,108 + ,30 + ,4 + ,103 + ,101193 + ,149061 + ,43 + ,26 + ,0 + ,93 + ,116174 + ,237213 + ,78 + ,38 + ,0 + ,123 + ,66198 + ,173326 + ,86 + ,44 + ,-4 + ,148 + ,71701 + ,133131 + ,44 + ,30 + ,4 + ,90 + ,57793 + ,258873 + ,104 + ,40 + ,4 + ,124 + ,80444 + ,324799 + ,158 + ,47 + ,0 + ,168 + ,97668 + ,230964 + ,102 + ,30 + ,-1 + ,115 + ,133824 + ,236785 + ,77 + ,31 + ,0 + ,71 + ,101481 + ,344297 + ,80 + ,30 + ,1 + ,108 + ,67654 + ,174724 + ,123 + ,34 + ,0 + ,120 + ,69112 + ,174415 + ,73 + ,31 + ,3 + ,114 + ,82753 + ,223632 + ,105 + ,33 + ,-1 + ,120 + ,72654 + ,294424 + ,107 + ,33 + ,4 + ,124 + ,101494 + ,325107 + ,84 + ,36 + ,3 + ,126 + ,79215 + ,106408 + ,33 + ,14 + ,1 + ,37 + ,31081 + ,96560 + ,42 + ,17 + ,0 + ,38 + ,22996 + ,265769 + ,96 + ,32 + ,-2 + ,120 + ,83122 + ,269651 + ,106 + ,30 + ,-3 + ,93 + ,70106 + ,149112 + ,56 + ,35 + ,-4 + ,95 + ,60578 + ,152871 + ,59 + ,28 + ,2 + ,90 + ,79892 + ,362301 + ,76 + ,34 + ,2 + ,110 + ,100708 + ,183167 + ,91 + ,39 + ,-4 + ,138 + ,82875 + ,277965 + ,115 + ,39 + ,3 + ,133 + ,139077 + ,218946 + ,76 + ,29 + ,2 + ,96 + ,80670 + ,244052 + ,101 + ,44 + ,2 + ,164 + ,143558 + ,341570 + ,94 + ,21 + ,0 + ,78 + ,117105 + ,233328 + ,92 + ,28 + ,5 + ,102 + ,120733 + ,206161 + ,75 + ,28 + ,-2 + ,99 + ,73107 + ,311473 + ,128 + ,38 + ,0 + ,129 + ,132068 + ,207176 + ,56 + ,32 + ,-2 + ,114 + ,87011 + ,196553 + ,41 + ,29 + ,-3 + ,99 + ,95260 + ,143246 + ,67 + ,27 + ,2 + ,104 + ,106671 + ,182192 + ,77 + ,40 + ,2 + ,138 + ,70054 + ,194979 + ,66 + ,40 + ,2 + ,151 + ,74011 + ,167488 + ,69 + ,28 + ,0 + ,72 + ,83737 + ,143756 + ,105 + ,34 + ,4 + ,120 + ,69094 + ,275541 + ,116 + ,33 + ,4 + ,115 + ,93133 + ,152299 + ,62 + ,33 + ,2 + ,98 + ,61370 + ,193339 + ,100 + ,35 + ,2 + ,71 + ,84651 + ,130585 + ,67 + ,29 + ,-4 + ,107 + ,95364 + ,112611 + ,46 + ,20 + ,3 + ,73 + ,26706 + ,148446 + ,135 + ,37 + ,3 + ,129 + ,126846 + ,182079 + ,124 + ,33 + ,2 + ,118 + ,102860 + ,243060 + ,58 + ,29 + ,-1 + ,104 + ,111813 + ,162765 + ,68 + ,28 + ,-3 + ,107 + ,120293 + ,85574 + ,37 + ,21 + ,0 + ,36 + ,24266 + ,225060 + ,93 + ,41 + ,1 + ,139 + ,109825 + ,133328 + ,56 + ,20 + ,-3 + ,56 + ,40909 + ,100750 + ,83 + ,30 + ,3 + ,93 + ,140867 + ,101523 + ,59 + ,22 + ,0 + ,87 + ,61056 + ,243511 + ,133 + ,42 + ,0 + ,110 + ,101338 + ,152474 + ,106 + ,32 + ,0 + ,83 + ,65567 + ,132487 + ,71 + ,36 + ,3 + ,98 + ,40735 + ,317394 + ,116 + ,31 + ,-3 + ,82 + ,91413 + ,244749 + ,98 + ,33 + ,0 + ,115 + ,76643 + ,184510 + ,64 + ,40 + ,-4 + ,140 + ,110681 + ,128423 + ,32 + ,38 + ,2 + ,120 + ,92696 + ,97839 + ,25 + ,24 + ,-1 + ,66 + ,94785 + ,172494 + ,46 + ,43 + ,3 + ,139 + ,86687 + ,229242 + ,63 + ,31 + ,2 + ,119 + ,91721 + ,351619 + ,95 + ,40 + ,5 + ,141 + ,115168 + ,324598 + ,113 + ,37 + ,2 + ,133 + ,135777 + ,195838 + ,111 + ,31 + ,-2 + ,98 + ,102372 + ,254488 + ,120 + ,39 + ,0 + ,117 + ,103772 + ,199476 + ,87 + ,32 + ,3 + ,105 + ,135400 + ,92499 + ,25 + ,18 + ,-2 + ,55 + ,21399 + ,224330 + ,131 + ,39 + ,0 + ,132 + ,130115 + ,181633 + ,47 + ,30 + ,6 + ,73 + ,64466 + ,271856 + ,109 + ,37 + ,-3 + ,86 + ,54990 + ,95227 + ,37 + ,32 + ,3 + ,48 + ,34777 + ,98146 + ,15 + ,17 + ,0 + ,48 + ,27114 + ,118612 + ,54 + ,12 + ,-2 + ,43 + ,30080 + ,65475 + ,16 + ,13 + ,1 + ,46 + ,69008 + ,108446 + ,22 + ,17 + ,0 + ,65 + ,46300 + ,121848 + ,37 + ,17 + ,2 + ,52 + ,30594 + ,76302 + ,29 + ,20 + ,2 + ,68 + ,30976 + ,98104 + ,55 + ,17 + ,-3 + ,47 + ,25568 + ,30989 + ,5 + ,17 + ,-2 + ,41 + ,4154 + ,31774 + ,0 + ,17 + ,1 + ,47 + ,4143 + ,150580 + ,27 + ,22 + ,-4 + ,71 + ,45588 + ,54157 + ,37 + ,15 + ,0 + ,30 + ,18625 + ,59382 + ,29 + ,12 + ,1 + ,24 + ,26263 + ,84105 + ,17 + ,17 + ,0 + ,63 + ,20055) + ,dim=c(6 + ,85) + ,dimnames=list(c('time_in_rfc' + ,'blogged_computations' + ,'compendiums_reviewed' + ,'SCORE' + ,'feedback_messages_p120' + ,'totsize') + ,1:85)) > y <- array(NA,dim=c(6,85),dimnames=list(c('time_in_rfc','blogged_computations','compendiums_reviewed','SCORE','feedback_messages_p120','totsize'),1:85)) > 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 = '6' > par2 = 'none' > par1 = '2' > 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] "blogged_computations" > x[,par1] [1] 79 108 43 78 86 44 104 158 102 77 80 123 73 105 107 84 33 42 96 [20] 106 56 59 76 91 115 76 101 94 92 75 128 56 41 67 77 66 69 105 [39] 116 62 100 67 46 135 124 58 68 37 93 56 83 59 133 106 71 116 98 [58] 64 32 25 46 63 95 113 111 120 87 25 131 47 109 37 15 54 16 22 [77] 37 29 55 5 0 27 37 29 17 > 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]) 0 5 15 16 17 22 25 27 29 32 33 37 41 42 43 44 46 47 54 55 1 1 1 1 1 1 2 1 2 1 1 4 1 1 1 1 2 1 1 1 56 58 59 62 63 64 66 67 68 69 71 73 75 76 77 78 79 80 83 84 3 1 2 1 1 1 1 2 1 1 1 1 1 2 2 1 1 1 1 1 86 87 91 92 93 94 95 96 98 100 101 102 104 105 106 107 108 109 111 113 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 115 116 120 123 124 128 131 133 135 158 1 2 1 1 1 1 1 1 1 1 > colnames(x) [1] "time_in_rfc" "blogged_computations" "compendiums_reviewed" [4] "SCORE" "feedback_messages_p120" "totsize" > colnames(x)[par1] [1] "blogged_computations" > x[,par1] [1] 79 108 43 78 86 44 104 158 102 77 80 123 73 105 107 84 33 42 96 [20] 106 56 59 76 91 115 76 101 94 92 75 128 56 41 67 77 66 69 105 [39] 116 62 100 67 46 135 124 58 68 37 93 56 83 59 133 106 71 116 98 [58] 64 32 25 46 63 95 113 111 120 87 25 131 47 109 37 15 54 16 22 [77] 37 29 55 5 0 27 37 29 17 > if (par2 == 'none') { + m <- ctree(as.formula(paste(colnames(x)[par1],' ~ .',sep='')),data = x) + } > > #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/1g2rf1323702174.tab") + } + } > m Conditional inference tree with 4 terminal nodes Response: blogged_computations Inputs: time_in_rfc, compendiums_reviewed, SCORE, feedback_messages_p120, totsize Number of observations: 85 1) time_in_rfc <= 133328; criterion = 1, statistic = 43.934 2) time_in_rfc <= 98146; criterion = 0.994, statistic = 10.542 3)* weights = 14 2) time_in_rfc > 98146 4)* weights = 12 1) time_in_rfc > 133328 5) compendiums_reviewed <= 29; criterion = 0.995, statistic = 10.951 6)* weights = 12 5) compendiums_reviewed > 29 7)* weights = 47 > postscript(file="/var/wessaorg/rcomp/tmp/2lgke1323702174.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/wessaorg/rcomp/tmp/3jsxm1323702174.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 79 95.87234 -16.8723404 2 108 95.87234 12.1276596 3 43 64.08333 -21.0833333 4 78 95.87234 -17.8723404 5 86 95.87234 -9.8723404 6 44 50.33333 -6.3333333 7 104 95.87234 8.1276596 8 158 95.87234 62.1276596 9 102 95.87234 6.1276596 10 77 95.87234 -18.8723404 11 80 95.87234 -15.8723404 12 123 95.87234 27.1276596 13 73 95.87234 -22.8723404 14 105 95.87234 9.1276596 15 107 95.87234 11.1276596 16 84 95.87234 -11.8723404 17 33 50.33333 -17.3333333 18 42 26.35714 15.6428571 19 96 95.87234 0.1276596 20 106 95.87234 10.1276596 21 56 95.87234 -39.8723404 22 59 64.08333 -5.0833333 23 76 95.87234 -19.8723404 24 91 95.87234 -4.8723404 25 115 95.87234 19.1276596 26 76 64.08333 11.9166667 27 101 95.87234 5.1276596 28 94 64.08333 29.9166667 29 92 64.08333 27.9166667 30 75 64.08333 10.9166667 31 128 95.87234 32.1276596 32 56 95.87234 -39.8723404 33 41 64.08333 -23.0833333 34 67 64.08333 2.9166667 35 77 95.87234 -18.8723404 36 66 95.87234 -29.8723404 37 69 64.08333 4.9166667 38 105 95.87234 9.1276596 39 116 95.87234 20.1276596 40 62 95.87234 -33.8723404 41 100 95.87234 4.1276596 42 67 50.33333 16.6666667 43 46 50.33333 -4.3333333 44 135 95.87234 39.1276596 45 124 95.87234 28.1276596 46 58 64.08333 -6.0833333 47 68 64.08333 3.9166667 48 37 26.35714 10.6428571 49 93 95.87234 -2.8723404 50 56 50.33333 5.6666667 51 83 50.33333 32.6666667 52 59 50.33333 8.6666667 53 133 95.87234 37.1276596 54 106 95.87234 10.1276596 55 71 50.33333 20.6666667 56 116 95.87234 20.1276596 57 98 95.87234 2.1276596 58 64 95.87234 -31.8723404 59 32 50.33333 -18.3333333 60 25 26.35714 -1.3571429 61 46 95.87234 -49.8723404 62 63 95.87234 -32.8723404 63 95 95.87234 -0.8723404 64 113 95.87234 17.1276596 65 111 95.87234 15.1276596 66 120 95.87234 24.1276596 67 87 95.87234 -8.8723404 68 25 26.35714 -1.3571429 69 131 95.87234 35.1276596 70 47 95.87234 -48.8723404 71 109 95.87234 13.1276596 72 37 26.35714 10.6428571 73 15 26.35714 -11.3571429 74 54 50.33333 3.6666667 75 16 26.35714 -10.3571429 76 22 50.33333 -28.3333333 77 37 50.33333 -13.3333333 78 29 26.35714 2.6428571 79 55 26.35714 28.6428571 80 5 26.35714 -21.3571429 81 0 26.35714 -26.3571429 82 27 64.08333 -37.0833333 83 37 26.35714 10.6428571 84 29 26.35714 2.6428571 85 17 26.35714 -9.3571429 > if (par2 != 'none') { + print(cbind(as.factor(x[,par1]),predict(m))) + myt <- table(as.factor(x[,par1]),predict(m)) + print(myt) + } > postscript(file="/var/wessaorg/rcomp/tmp/4u6fj1323702174.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/wessaorg/rcomp/tmp/5ipsp1323702174.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/wessaorg/rcomp/tmp/6950z1323702174.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/wessaorg/rcomp/tmp/78a7i1323702174.tab") + } > > try(system("convert tmp/2lgke1323702174.ps tmp/2lgke1323702174.png",intern=TRUE)) character(0) > try(system("convert tmp/3jsxm1323702174.ps tmp/3jsxm1323702174.png",intern=TRUE)) character(0) > try(system("convert tmp/4u6fj1323702174.ps tmp/4u6fj1323702174.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.761 0.220 3.004