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(296.95,17.20,296.84,17.20,287.54,17.20,287.81,17.20,283.99,20.63,275.79,20.63,269.52,20.63,278.35,20.63,283.43,19.32,289.46,19.32,282.30,19.32,293.55,19.32,304.78,12.99,300.99,12.99,315.29,12.99,316.21,12.99,331.79,18.13,329.38,18.13,317.27,18.13,317.98,18.13,340.28,28.37,339.21,28.37,336.71,28.37,340.11,28.37,347.72,24.35,328.68,24.35,303.05,24.35,299.83,24.35,320.04,24.99,317.94,24.99,303.31,24.99,308.85,24.99,319.19,28.84,314.52,28.84,312.39,28.84,315.77,28.84,320.23,37.88,309.45,37.88,296.54,37.88,297.28,37.88,301.39,54.04,306.68,54.04,305.91,54.04,314.76,54.04,323.34,64.93,341.58,64.93,330.12,64.93,318.16,64.93,317.84,71.81,325.39,71.81,327.56,71.81,329.77,71.81,333.29,99.75,346.10,99.75,358.00,99.75,344.82,99.75,313.30,61.25,301.26,61.25,306.38,61.25,319.31,61.25),dim=c(2,60),dimnames=list(c('GemPrijsTicket_$','GemOlieprijs_$'),1:60)) > y <- array(NA,dim=c(2,60),dimnames=list(c('GemPrijsTicket_$','GemOlieprijs_$'),1:60)) > 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] "GemPrijsTicket_." > x[,par1] [1] 296.95 296.84 287.54 287.81 283.99 275.79 269.52 278.35 283.43 289.46 [11] 282.30 293.55 304.78 300.99 315.29 316.21 331.79 329.38 317.27 317.98 [21] 340.28 339.21 336.71 340.11 347.72 328.68 303.05 299.83 320.04 317.94 [31] 303.31 308.85 319.19 314.52 312.39 315.77 320.23 309.45 296.54 297.28 [41] 301.39 306.68 305.91 314.76 323.34 341.58 330.12 318.16 317.84 325.39 [51] 327.56 329.77 333.29 346.10 358.00 344.82 313.30 301.26 306.38 319.31 > 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]) 269.52 275.79 278.35 282.3 283.43 283.99 287.54 287.81 289.46 293.55 296.54 1 1 1 1 1 1 1 1 1 1 1 296.84 296.95 297.28 299.83 300.99 301.26 301.39 303.05 303.31 304.78 305.91 1 1 1 1 1 1 1 1 1 1 1 306.38 306.68 308.85 309.45 312.39 313.3 314.52 314.76 315.29 315.77 316.21 1 1 1 1 1 1 1 1 1 1 1 317.27 317.84 317.94 317.98 318.16 319.19 319.31 320.04 320.23 323.34 325.39 1 1 1 1 1 1 1 1 1 1 1 327.56 328.68 329.38 329.77 330.12 331.79 333.29 336.71 339.21 340.11 340.28 1 1 1 1 1 1 1 1 1 1 1 341.58 344.82 346.1 347.72 358 1 1 1 1 1 > colnames(x) [1] "GemPrijsTicket_." "GemOlieprijs_." > colnames(x)[par1] [1] "GemPrijsTicket_." > x[,par1] [1] 296.95 296.84 287.54 287.81 283.99 275.79 269.52 278.35 283.43 289.46 [11] 282.30 293.55 304.78 300.99 315.29 316.21 331.79 329.38 317.27 317.98 [21] 340.28 339.21 336.71 340.11 347.72 328.68 303.05 299.83 320.04 317.94 [31] 303.31 308.85 319.19 314.52 312.39 315.77 320.23 309.45 296.54 297.28 [41] 301.39 306.68 305.91 314.76 323.34 341.58 330.12 318.16 317.84 325.39 [51] 327.56 329.77 333.29 346.10 358.00 344.82 313.30 301.26 306.38 319.31 > 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/1lxor1292944265.tab") + } + } > m Conditional inference tree with 4 terminal nodes Response: GemPrijsTicket_. Input: GemOlieprijs_. Number of observations: 60 1) GemOlieprijs_. <= 20.63; criterion = 1, statistic = 15.637 2) GemOlieprijs_. <= 18.13; criterion = 0.976, statistic = 5.064 3)* weights = 12 2) GemOlieprijs_. > 18.13 4)* weights = 8 1) GemOlieprijs_. > 20.63 5) GemOlieprijs_. <= 61.25; criterion = 0.981, statistic = 5.523 6)* weights = 28 5) GemOlieprijs_. > 61.25 7)* weights = 12 > postscript(file="/var/www/html/freestat/rcomp/tmp/2lxor1292944265.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/3lxor1292944265.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 296.95 308.5692 -11.619167 2 296.84 308.5692 -11.729167 3 287.54 308.5692 -21.029167 4 287.81 308.5692 -20.759167 5 283.99 282.0487 1.941250 6 275.79 282.0487 -6.258750 7 269.52 282.0487 -12.528750 8 278.35 282.0487 -3.698750 9 283.43 282.0487 1.381250 10 289.46 282.0487 7.411250 11 282.30 282.0487 0.251250 12 293.55 282.0487 11.501250 13 304.78 308.5692 -3.789167 14 300.99 308.5692 -7.579167 15 315.29 308.5692 6.720833 16 316.21 308.5692 7.640833 17 331.79 308.5692 23.220833 18 329.38 308.5692 20.810833 19 317.27 308.5692 8.700833 20 317.98 308.5692 9.410833 21 340.28 315.7175 24.562500 22 339.21 315.7175 23.492500 23 336.71 315.7175 20.992500 24 340.11 315.7175 24.392500 25 347.72 315.7175 32.002500 26 328.68 315.7175 12.962500 27 303.05 315.7175 -12.667500 28 299.83 315.7175 -15.887500 29 320.04 315.7175 4.322500 30 317.94 315.7175 2.222500 31 303.31 315.7175 -12.407500 32 308.85 315.7175 -6.867500 33 319.19 315.7175 3.472500 34 314.52 315.7175 -1.197500 35 312.39 315.7175 -3.327500 36 315.77 315.7175 0.052500 37 320.23 315.7175 4.512500 38 309.45 315.7175 -6.267500 39 296.54 315.7175 -19.177500 40 297.28 315.7175 -18.437500 41 301.39 315.7175 -14.327500 42 306.68 315.7175 -9.037500 43 305.91 315.7175 -9.807500 44 314.76 315.7175 -0.957500 45 323.34 332.9975 -9.657500 46 341.58 332.9975 8.582500 47 330.12 332.9975 -2.877500 48 318.16 332.9975 -14.837500 49 317.84 332.9975 -15.157500 50 325.39 332.9975 -7.607500 51 327.56 332.9975 -5.437500 52 329.77 332.9975 -3.227500 53 333.29 332.9975 0.292500 54 346.10 332.9975 13.102500 55 358.00 332.9975 25.002500 56 344.82 332.9975 11.822500 57 313.30 315.7175 -2.417500 58 301.26 315.7175 -14.457500 59 306.38 315.7175 -9.337500 60 319.31 315.7175 3.592500 > 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/4d6nu1292944265.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/5syk21292944265.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/6kpk51292944265.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/7oqib1292944265.tab") + } > > try(system("convert tmp/2lxor1292944265.ps tmp/2lxor1292944265.png",intern=TRUE)) character(0) > try(system("convert tmp/3lxor1292944265.ps tmp/3lxor1292944265.png",intern=TRUE)) character(0) > try(system("convert tmp/4d6nu1292944265.ps tmp/4d6nu1292944265.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.849 0.752 5.161