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(12008.00 + ,4.00 + ,9169.00 + ,5.90 + ,8788.00 + ,7.10 + ,8417.00 + ,10.50 + ,8247.00 + ,15.10 + ,8197.00 + ,16.80 + ,8236.00 + ,15.30 + ,8253.00 + ,18.40 + ,7733.00 + ,16.10 + ,8366.00 + ,11.30 + ,8626.00 + ,7.90 + ,8863.00 + ,5.60 + ,10102.00 + ,3.40 + ,8463.00 + ,4.80 + ,9114.00 + ,6.50 + ,8563.00 + ,8.50 + ,8872.00 + ,15.10 + ,8301.00 + ,15.70 + ,8301.00 + ,18.70 + ,8278.00 + ,19.20 + ,7736.00 + ,12.90 + ,7973.00 + ,14.40 + ,8268.00 + ,6.20 + ,9476.00 + ,3.30 + ,11100.00 + ,4.60 + ,8962.00 + ,7.10 + ,9173.00 + ,7.80 + ,8738.00 + ,9.90 + ,8459.00 + ,13.60 + ,8078.00 + ,17.10 + ,8411.00 + ,17.80 + ,8291.00 + ,18.60 + ,7810.00 + ,14.70 + ,8616.00 + ,10.50 + ,8312.00 + ,8.60 + ,9692.00 + ,4.40 + ,9911.00 + ,2.30 + ,8915.00 + ,2.80 + ,9452.00 + ,8.80 + ,9112.00 + ,10.70 + ,8472.00 + ,13.90 + ,8230.00 + ,19.30 + ,8384.00 + ,19.50 + ,8625.00 + ,20.40 + ,8221.00 + ,15.30 + ,8649.00 + ,7.90 + ,8625.00 + ,8.30 + ,10443.00 + ,4.50 + ,10357.00 + ,3.20 + ,8586.00 + ,5.00 + ,8892.00 + ,6.60 + ,8329.00 + ,11.10 + ,8101.00 + ,12.80 + ,7922.00 + ,16.30 + ,8120.00 + ,17.40 + ,7838.00 + ,18.90 + ,7735.00 + ,15.80 + ,8406.00 + ,11.70 + ,8209.00 + ,6.40 + ,9451.00 + ,2.90 + ,10041.00 + ,4.70 + ,9411.00 + ,2.40 + ,10405.00 + ,7.20 + ,8467.00 + ,10.70 + ,8464.00 + ,13.40 + ,8102.00 + ,18.30 + ,7627.00 + ,18.40 + ,7513.00 + ,16.80 + ,7510.00 + ,16.60 + ,8291.00 + ,14.10 + ,8064.00 + ,6.10 + ,9383.00 + ,3.50 + ,9706.00 + ,1.70 + ,8579.00 + ,2.30 + ,9474.00 + ,4.50 + ,8318.00 + ,9.30 + ,8213.00 + ,14.20 + ,8059.00 + ,17.30 + ,9111.00 + ,23.00 + ,7708.00 + ,16.30 + ,7680.00 + ,18.40 + ,8014.00 + ,14.20 + ,8007.00 + ,9.10 + ,8718.00 + ,5.90 + ,9486.00 + ,7.20 + ,9113.00 + ,6.80 + ,9025.00 + ,8.00 + ,8476.00 + ,14.30 + ,7952.00 + ,14.60 + ,7759.00 + ,17.50 + ,7835.00 + ,17.20 + ,7600.00 + ,17.20 + ,7651.00 + ,14.10 + ,8319.00 + ,10.40 + ,8812.00 + ,6.80 + ,8630.00 + ,4.10) + ,dim=c(2 + ,96) + ,dimnames=list(c('Sterftecijfers' + ,'Temperatuur') + ,1:96)) > y <- array(NA,dim=c(2,96),dimnames=list(c('Sterftecijfers','Temperatuur'),1:96)) > 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 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] "Sterftecijfers" > x[,par1] [1] 12008 9169 8788 8417 8247 8197 8236 8253 7733 8366 8626 8863 [13] 10102 8463 9114 8563 8872 8301 8301 8278 7736 7973 8268 9476 [25] 11100 8962 9173 8738 8459 8078 8411 8291 7810 8616 8312 9692 [37] 9911 8915 9452 9112 8472 8230 8384 8625 8221 8649 8625 10443 [49] 10357 8586 8892 8329 8101 7922 8120 7838 7735 8406 8209 9451 [61] 10041 9411 10405 8467 8464 8102 7627 7513 7510 8291 8064 9383 [73] 9706 8579 9474 8318 8213 8059 9111 7708 7680 8014 8007 8718 [85] 9486 9113 9025 8476 7952 7759 7835 7600 7651 8319 8812 8630 > 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]) 7510 7513 7600 7627 7651 7680 7708 7733 7735 7736 7759 7810 7835 1 1 1 1 1 1 1 1 1 1 1 1 1 7838 7922 7952 7973 8007 8014 8059 8064 8078 8101 8102 8120 8197 1 1 1 1 1 1 1 1 1 1 1 1 1 8209 8213 8221 8230 8236 8247 8253 8268 8278 8291 8301 8312 8318 1 1 1 1 1 1 1 1 1 2 2 1 1 8319 8329 8366 8384 8406 8411 8417 8459 8463 8464 8467 8472 8476 1 1 1 1 1 1 1 1 1 1 1 1 1 8563 8579 8586 8616 8625 8626 8630 8649 8718 8738 8788 8812 8863 1 1 1 1 2 1 1 1 1 1 1 1 1 8872 8892 8915 8962 9025 9111 9112 9113 9114 9169 9173 9383 9411 1 1 1 1 1 1 1 1 1 1 1 1 1 9451 9452 9474 9476 9486 9692 9706 9911 10041 10102 10357 10405 10443 1 1 1 1 1 1 1 1 1 1 1 1 1 11100 12008 1 1 > colnames(x) [1] "Sterftecijfers" "Temperatuur" > colnames(x)[par1] [1] "Sterftecijfers" > x[,par1] [1] 12008 9169 8788 8417 8247 8197 8236 8253 7733 8366 8626 8863 [13] 10102 8463 9114 8563 8872 8301 8301 8278 7736 7973 8268 9476 [25] 11100 8962 9173 8738 8459 8078 8411 8291 7810 8616 8312 9692 [37] 9911 8915 9452 9112 8472 8230 8384 8625 8221 8649 8625 10443 [49] 10357 8586 8892 8329 8101 7922 8120 7838 7735 8406 8209 9451 [61] 10041 9411 10405 8467 8464 8102 7627 7513 7510 8291 8064 9383 [73] 9706 8579 9474 8318 8213 8059 9111 7708 7680 8014 8007 8718 [85] 9486 9113 9025 8476 7952 7759 7835 7600 7651 8319 8812 8630 > 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/1f3fk1323606748.tab") + } + } > m Conditional inference tree with 3 terminal nodes Response: Sterftecijfers Input: Temperatuur Number of observations: 96 1) Temperatuur <= 4.7; criterion = 1, statistic = 42.801 2)* weights = 17 1) Temperatuur > 4.7 3) Temperatuur <= 10.7; criterion = 1, statistic = 23.557 4)* weights = 32 3) Temperatuur > 10.7 5)* weights = 47 > postscript(file="/var/www/rcomp/tmp/2hcx21323606748.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/34g5k1323606748.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 12008 9804.647 2203.352941 2 9169 8760.344 408.656250 3 8788 8760.344 27.656250 4 8417 8760.344 -343.343750 5 8247 8116.809 130.191489 6 8197 8116.809 80.191489 7 8236 8116.809 119.191489 8 8253 8116.809 136.191489 9 7733 8116.809 -383.808511 10 8366 8116.809 249.191489 11 8626 8760.344 -134.343750 12 8863 8760.344 102.656250 13 10102 9804.647 297.352941 14 8463 8760.344 -297.343750 15 9114 8760.344 353.656250 16 8563 8760.344 -197.343750 17 8872 8116.809 755.191489 18 8301 8116.809 184.191489 19 8301 8116.809 184.191489 20 8278 8116.809 161.191489 21 7736 8116.809 -380.808511 22 7973 8116.809 -143.808511 23 8268 8760.344 -492.343750 24 9476 9804.647 -328.647059 25 11100 9804.647 1295.352941 26 8962 8760.344 201.656250 27 9173 8760.344 412.656250 28 8738 8760.344 -22.343750 29 8459 8116.809 342.191489 30 8078 8116.809 -38.808511 31 8411 8116.809 294.191489 32 8291 8116.809 174.191489 33 7810 8116.809 -306.808511 34 8616 8760.344 -144.343750 35 8312 8760.344 -448.343750 36 9692 9804.647 -112.647059 37 9911 9804.647 106.352941 38 8915 9804.647 -889.647059 39 9452 8760.344 691.656250 40 9112 8760.344 351.656250 41 8472 8116.809 355.191489 42 8230 8116.809 113.191489 43 8384 8116.809 267.191489 44 8625 8116.809 508.191489 45 8221 8116.809 104.191489 46 8649 8760.344 -111.343750 47 8625 8760.344 -135.343750 48 10443 9804.647 638.352941 49 10357 9804.647 552.352941 50 8586 8760.344 -174.343750 51 8892 8760.344 131.656250 52 8329 8116.809 212.191489 53 8101 8116.809 -15.808511 54 7922 8116.809 -194.808511 55 8120 8116.809 3.191489 56 7838 8116.809 -278.808511 57 7735 8116.809 -381.808511 58 8406 8116.809 289.191489 59 8209 8760.344 -551.343750 60 9451 9804.647 -353.647059 61 10041 9804.647 236.352941 62 9411 9804.647 -393.647059 63 10405 8760.344 1644.656250 64 8467 8760.344 -293.343750 65 8464 8116.809 347.191489 66 8102 8116.809 -14.808511 67 7627 8116.809 -489.808511 68 7513 8116.809 -603.808511 69 7510 8116.809 -606.808511 70 8291 8116.809 174.191489 71 8064 8760.344 -696.343750 72 9383 9804.647 -421.647059 73 9706 9804.647 -98.647059 74 8579 9804.647 -1225.647059 75 9474 9804.647 -330.647059 76 8318 8760.344 -442.343750 77 8213 8116.809 96.191489 78 8059 8116.809 -57.808511 79 9111 8116.809 994.191489 80 7708 8116.809 -408.808511 81 7680 8116.809 -436.808511 82 8014 8116.809 -102.808511 83 8007 8760.344 -753.343750 84 8718 8760.344 -42.343750 85 9486 8760.344 725.656250 86 9113 8760.344 352.656250 87 9025 8760.344 264.656250 88 8476 8116.809 359.191489 89 7952 8116.809 -164.808511 90 7759 8116.809 -357.808511 91 7835 8116.809 -281.808511 92 7600 8116.809 -516.808511 93 7651 8116.809 -465.808511 94 8319 8760.344 -441.343750 95 8812 8760.344 51.656250 96 8630 9804.647 -1174.647059 > 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/4iv4r1323606748.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/57im81323606748.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/6dd4y1323606748.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/70qf11323606748.tab") + } > > try(system("convert tmp/2hcx21323606748.ps tmp/2hcx21323606748.png",intern=TRUE)) character(0) > try(system("convert tmp/34g5k1323606748.ps tmp/34g5k1323606748.png",intern=TRUE)) character(0) > try(system("convert tmp/4iv4r1323606748.ps tmp/4iv4r1323606748.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.12 0.13 2.23