x <- array(list(0 ,3289.5 ,0.66 ,0.814 ,0.526 ,1 ,6 ,25299.2 ,0.384 ,0.743 ,0.587 ,2 ,0 ,52.8 ,NA ,0.891 ,0.799 ,NA ,0 ,10335.1 ,NA ,0.334 ,0.464 ,2 ,0 ,62.2 ,NA ,0.771 ,0.676 ,1 ,1 ,32642.4 ,0.681 ,0.813 ,0.612 ,0 ,64 ,17096.2 ,0.952 ,0.895 ,0.779 ,0 ,4 ,7670.5 ,0.709 ,0.875 ,0.793 ,0 ,1 ,256.1 ,NA ,0.778 ,0.778 ,0 ,0 ,492.9 ,0.602 ,0.825 ,0.754 ,1 ,0 ,105256 ,0.253 ,0.623 ,0.277 ,0 ,0 ,259.5 ,NA ,0.862 ,0.72 ,0 ,4 ,9949 ,0.764 ,0.884 ,0.791 ,0 ,0 ,190.2 ,0.615 ,0.828 ,0.54 ,0 ,0 ,4773.1 ,0.203 ,0.453 ,0.342 ,0 ,0 ,558.5 ,NA ,0.513 ,0.415 ,2 ,0 ,6658.5 ,0.604 ,0.612 ,0.476 ,0 ,0 ,4308.2 ,NA ,0.737 ,NA ,2 ,0 ,1382.4 ,0.497 ,0.696 ,0.607 ,0 ,12 ,149650.2 ,NA ,0.73 ,0.608 ,0 ,0 ,9324 ,NA ,0.45 ,0.283 ,1 ,0 ,12180.8 ,0.346 ,0.525 ,0.428 ,2 ,47 ,27700.9 ,0.875 ,0.904 ,0.796 ,0 ,0 ,2934.8 ,0.219 ,0.455 ,0.3 ,1 ,0 ,6011.2 ,NA ,0.484 ,0.3 ,2 ,0 ,13187.8 ,0.677 ,0.847 ,0.592 ,0 ,1 ,33203.3 ,0.39 ,0.762 ,0.582 ,1 ,0 ,36406.2 ,0.235 ,0.424 ,NA ,2 ,0 ,3070.2 ,0.548 ,0.879 ,0.587 ,0 ,4 ,4517.2 ,0.624 ,0.819 ,NA ,1 ,93 ,10570.3 ,0.69 ,0.859 ,0.523 ,2 ,0 ,766.7 ,0.626 ,0.892 ,0.747 ,0 ,25 ,10302.7 ,0.762 ,0.819 ,NA ,NA ,25 ,20143.2 ,NA ,0.801 ,NA ,2 ,11 ,5141 ,0.774 ,0.866 ,0.79 ,0 ,0 ,562.3 ,NA ,0.495 ,NA ,2 ,0 ,7194.7 ,0.499 ,0.751 ,0.513 ,0 ,0 ,10260.6 ,0.486 ,0.771 ,0.557 ,0 ,0 ,56843.3 ,0.376 ,0.663 ,0.493 ,1 ,0 ,5332.8 ,0.386 ,0.726 ,0.513 ,1 ,0 ,373.9 ,NA ,0.422 ,0.419 ,2 ,6 ,1567.6 ,0.716 ,0.78 ,0.661 ,1 ,7 ,48333.3 ,NA ,0.427 ,0.245 ,1 ,0 ,728.4 ,0.674 ,0.718 ,0.503 ,1 ,11 ,4986.4 ,0.74 ,0.871 ,0.776 ,0 ,66 ,56708.3 ,0.666 ,0.895 ,0.787 ,0 ,0 ,929 ,0.484 ,0.652 ,0.702 ,1 ,235 ,79098.1 ,0.721 ,0.874 ,0.796 ,0 ,1 ,14793.4 ,0.406 ,0.581 ,0.31 ,1 ,43 ,57214.5 ,0.688 ,0.877 ,0.781 ,0 ,10 ,10160.5 ,0.671 ,0.9 ,0.742 ,0 ,0 ,96.2 ,NA ,0.774 ,0.554 ,0 ,0 ,8923.1 ,0.296 ,0.666 ,0.499 ,1 ,0 ,5759.4 ,NA ,0.373 ,0.292 ,1 ,0 ,724.9 ,0.541 ,0.646 ,0.335 ,1 ,0 ,7124.9 ,0.294 ,0.554 ,0.384 ,2 ,0 ,4889.3 ,0.405 ,0.731 ,0.457 ,0 ,0 ,5793.9 ,0.687 ,0.907 ,0.779 ,NA ,86 ,10376.3 ,0.661 ,0.778 ,0.683 ,0 ,0 ,254.8 ,0.727 ,0.917 ,0.79 ,0 ,0 ,873785.4 ,0.318 ,0.605 ,0.359 ,1 ,15 ,184345.9 ,0.39 ,0.664 ,0.43 ,1 ,0 ,17373.8 ,0.345 ,0.749 ,NA ,2 ,7 ,3531.2 ,0.76 ,0.863 ,0.729 ,0 ,4 ,54870.6 ,NA ,0.659 ,0.592 ,2 ,3 ,4499.9 ,0.784 ,0.891 ,0.739 ,0 ,48 ,56832.3 ,0.636 ,0.897 ,0.781 ,0 ,7 ,2364.9 ,0.592 ,0.8 ,0.571 ,0 ,42 ,122251.2 ,0.761 ,0.931 ,0.798 ,0 ,0 ,3415.6 ,0.526 ,0.795 ,0.493 ,1 ,20 ,23447.2 ,0.408 ,0.621 ,0.375 ,1 ,82 ,42980.4 ,0.738 ,0.816 ,0.678 ,0 ,0 ,2087.7 ,0.514 ,0.828 ,0.848 ,1 ,5 ,2663.9 ,0.652 ,0.772 ,0.661 ,1 ,0 ,2948.4 ,NA ,0.768 ,0.607 ,1 ,0 ,1639.2 ,0.431 ,0.623 ,0.388 ,1 ,0 ,29 ,NA ,0.885 ,0.886 ,0 ,6 ,3695.9 ,0.688 ,0.801 ,NA ,0 ,0 ,381.2 ,0.656 ,0.868 ,0.86 ,0 ,0 ,11280.6 ,NA ,0.485 ,0.329 ,1 ,0 ,9380.9 ,0.239 ,0.428 ,0.242 ,2 ,1 ,18208.6 ,0.534 ,0.789 ,0.595 ,1 ,0 ,219.5 ,0.447 ,0.645 ,NA ,2 ,0 ,8672.9 ,0.082 ,0.382 ,0.27 ,0 ,0 ,367.5 ,0.686 ,0.872 ,0.713 ,0 ,0 ,1995.5 ,0.194 ,0.568 ,0.4 ,2 ,0 ,1059.5 ,NA ,0.78 ,0.588 ,0 ,2 ,84306.6 ,0.518 ,0.802 ,0.657 ,1 ,0 ,30.9 ,NA ,0.902 ,NA ,NA ,2 ,2192.6 ,0.58 ,0.639 ,0.426 ,0 ,8 ,24781.1 ,0.254 ,0.696 ,0.466 ,1 ,0 ,13547.1 ,0.116 ,0.367 ,0.189 ,1 ,0 ,39268.3 ,0.258 ,0.588 ,0.175 ,NA ,4 ,1414.8 ,0.524 ,0.644 ,0.53 ,0 ,0 ,19081.1 ,0.257 ,0.536 ,0.287 ,0 ,29 ,14891.7 ,0.814 ,0.899 ,0.797 ,0 ,18 ,3398 ,0.873 ,0.871 ,0.745 ,0 ,0 ,7788.2 ,0.079 ,0.338 ,0.27 ,1 ,7 ,97552.1 ,NA ,0.404 ,0.364 ,1 ,19 ,4241.5 ,0.82 ,0.891 ,0.823 ,0 ,0 ,1868.1 ,0.71 ,0.798 ,0.718 ,1 ,1 ,111844.7 ,0.241 ,0.642 ,0.411 ,1 ,0 ,2415.9 ,0.599 ,0.825 ,0.581 ,1 ,0 ,4157.7 ,0.222 ,0.561 ,0.4 ,0 ,0 ,4243.9 ,0.466 ,0.757 ,0.532 ,1 ,140 ,1145195.2 ,0.437 ,0.78 ,0.345 ,2 ,2 ,21685.5 ,NA ,0.719 ,0.54 ,1 ,1 ,61628.7 ,0.578 ,0.712 ,0.452 ,1 ,37 ,38056.2 ,0.676 ,0.803 ,NA ,0 ,0 ,9925.5 ,0.568 ,0.856 ,0.729 ,0 ,40 ,23206.7 ,0.705 ,0.779 ,0.624 ,1 ,0 ,7109.5 ,0.217 ,0.202 ,0.286 ,2 ,0 ,161.3 ,NA ,0.708 ,NA ,0 ,0 ,161.3 ,NA ,0.708 ,NA ,0 ,0 ,24.1 ,NA ,0.925 ,NA ,0 ,0 ,16139.1 ,0.57 ,0.767 ,0.762 ,2 ,0 ,7241.6 ,0.247 ,0.524 ,0.375 ,1 ,0 ,71 ,NA ,0.794 ,0.698 ,1 ,0 ,3981.6 ,0.182 ,0.296 ,0.258 ,2 ,0 ,3016.6 ,NA ,0.877 ,0.788 ,1 ,2 ,1926.7 ,0.67 ,0.837 ,NA ,0 ,0 ,309.5 ,NA ,0.578 ,NA ,0 ,4 ,36793.9 ,0.572 ,0.655 ,0.622 ,1 ,81 ,38889.2 ,0.619 ,0.9 ,0.756 ,0 ,0 ,17337 ,0.59 ,0.781 ,0.43 ,1 ,0 ,107.4 ,0.665 ,0.777 ,0.536 ,0 ,0 ,26494.2 ,0.163 ,0.513 ,0.315 ,2 ,1 ,406.9 ,NA ,0.748 ,0.567 ,1 ,0 ,862.9 ,0.451 ,0.622 ,0.517 ,1 ,23 ,8558.8 ,0.759 ,0.909 ,0.786 ,0 ,5 ,6673.7 ,0.762 ,0.909 ,0.836 ,0 ,0 ,12324.1 ,0.425 ,0.806 ,0.48 ,2 ,1 ,57072.1 ,0.417 ,0.828 ,0.526 ,1 ,0 ,966.2 ,0.179 ,0.522 ,0.34 ,0 ,0 ,3665.5 ,0.312 ,0.521 ,0.306 ,2 ,0 ,95.2 ,NA ,0.782 ,0.491 ,1 ,0 ,1215.5 ,0.609 ,0.771 ,0.659 ,0 ,0 ,8215.1 ,0.396 ,0.767 ,0.524 ,1 ,16 ,54130.3 ,0.41 ,0.679 ,0.622 ,1 ,0 ,17699.7 ,0.262 ,0.431 ,0.235 ,2 ,0 ,1808.6 ,0.462 ,0.817 ,0.871 ,1 ,290 ,253339.1 ,0.917 ,0.871 ,0.825 ,0 ,0 ,3109.1 ,0.64 ,0.828 ,0.609 ,0 ,0 ,146.6 ,NA ,0.679 ,0.512 ,0 ,0 ,19685.2 ,0.476 ,0.805 ,0.651 ,1 ,0 ,67101.6 ,0.371 ,0.719 ,0.307 ,NA ,0 ,11948.2 ,0.096 ,0.571 ,NA ,NA ,0 ,7860.1 ,0.406 ,0.434 ,0.346 ,0 ,0 ,10469.2 ,0.451 ,0.64 ,0.266 ,1) ,dim=c(6 ,149) ,dimnames=list(c('Totalpoints' ,'POP' ,'Education' ,'Health' ,'GNI/Cap' ,'Democracy') ,1:149)) y <- array(NA,dim=c(6,149),dimnames=list(c('Totalpoints','POP','Education','Health','GNI/Cap','Democracy'),1:149)) 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' par4 <- 'no' par3 <- '3' par2 <- 'none' par1 <- '1' #'GNU S' R Code compiled by R2WASP v. 1.2.291 () #Author: root #To cite this work: Wessa P., 2012, Recursive Partitioning (Regression Trees) (v1.0.3) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_regression_trees.wasp/ #Source of accompanying publication: # library(party) library(Hmisc) par1 <- as.numeric(par1) par3 <- as.numeric(par3) x <- data.frame(t(y)) is.data.frame(x) x <- x[!is.na(x[,par1]),] k <- length(x[1,]) n <- length(x[,1]) colnames(x)[par1] x[,par1] 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]) colnames(x) colnames(x)[par1] x[,par1] 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/1yj6w1356614191.tab") } } m postscript(file="/var/wessaorg/rcomp/tmp/2r7861356614191.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) plot(m) dev.off() postscript(file="/var/wessaorg/rcomp/tmp/3bce01356614191.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() 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) } 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/4q2ks1356614191.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() 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/5dxda1356614191.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/6nsh51356614191.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/7oxdj1356614191.tab") } try(system("convert tmp/2r7861356614191.ps tmp/2r7861356614191.png",intern=TRUE)) try(system("convert tmp/3bce01356614191.ps tmp/3bce01356614191.png",intern=TRUE)) try(system("convert tmp/4q2ks1356614191.ps tmp/4q2ks1356614191.png",intern=TRUE))