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(16198.9 + ,16896.2 + ,0 + ,16554.2 + ,16698 + ,0 + ,19554.2 + ,19691.6 + ,0 + ,15903.8 + ,15930.7 + ,0 + ,18003.8 + ,17444.6 + ,0 + ,18329.6 + ,17699.4 + ,0 + ,16260.7 + ,15189.8 + ,0 + ,14851.9 + ,15672.7 + ,0 + ,18174.1 + ,17180.8 + ,0 + ,18406.6 + ,17664.9 + ,0 + ,18466.5 + ,17862.9 + ,0 + ,16016.5 + ,16162.3 + ,0 + ,17428.5 + ,17463.6 + ,0 + ,17167.2 + ,16772.1 + ,0 + ,19630 + ,19106.9 + ,0 + ,17183.6 + ,16721.3 + ,0 + ,18344.7 + ,18161.3 + ,0 + ,19301.4 + ,18509.9 + ,0 + ,18147.5 + ,17802.7 + ,0 + ,16192.9 + ,16409.9 + ,0 + ,18374.4 + ,17967.7 + ,0 + ,20515.2 + ,20286.6 + ,0 + ,18957.2 + ,19537.3 + ,0 + ,16471.5 + ,18021.9 + ,0 + ,18746.8 + ,20194.3 + ,0 + ,19009.5 + ,19049.6 + ,0 + ,19211.2 + ,20244.7 + ,0 + ,20547.7 + ,21473.3 + ,0 + ,19325.8 + ,19673.6 + ,0 + ,20605.5 + ,21053.2 + ,0 + ,20056.9 + ,20159.5 + ,0 + ,16141.4 + ,18203.6 + ,0 + ,20359.8 + ,21289.5 + ,0 + ,19711.6 + ,20432.3 + ,1 + ,15638.6 + ,17180.4 + ,1 + ,14384.5 + ,15816.8 + ,1 + ,13855.6 + ,15071.8 + ,1 + ,14308.3 + ,14521.1 + ,1 + ,15290.6 + ,15668.8 + ,1 + ,14423.8 + ,14346.9 + ,1 + ,13779.7 + ,13881 + ,1 + ,15686.3 + ,15465.9 + ,1 + ,14733.8 + ,14238.2 + ,1 + ,12522.5 + ,13557.7 + ,1 + ,16189.4 + ,16127.6 + ,1 + ,16059.1 + ,16793.9 + ,1 + ,16007.1 + ,16014 + ,1 + ,15806.8 + ,16867.9 + ,1 + ,15160 + ,16014.6 + ,0 + ,15692.1 + ,15878.6 + ,0 + ,18908.9 + ,18664.9 + ,0 + ,16969.9 + ,17962.5 + ,0 + ,16997.5 + ,17332.7 + ,0 + ,19858.9 + ,19542.1 + ,0 + ,17681.2 + ,17203.6 + ,0) + ,dim=c(3 + ,55) + ,dimnames=list(c('uitvoer' + ,'invoer' + ,'crisis') + ,1:55)) > y <- array(NA,dim=c(3,55),dimnames=list(c('uitvoer','invoer','crisis'),1:55)) > 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 = '0' > 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] "uitvoer" > x[,par1] [1] 16198.9 16554.2 19554.2 15903.8 18003.8 18329.6 16260.7 14851.9 18174.1 [10] 18406.6 18466.5 16016.5 17428.5 17167.2 19630.0 17183.6 18344.7 19301.4 [19] 18147.5 16192.9 18374.4 20515.2 18957.2 16471.5 18746.8 19009.5 19211.2 [28] 20547.7 19325.8 20605.5 20056.9 16141.4 20359.8 19711.6 15638.6 14384.5 [37] 13855.6 14308.3 15290.6 14423.8 13779.7 15686.3 14733.8 12522.5 16189.4 [46] 16059.1 16007.1 15806.8 15160.0 15692.1 18908.9 16969.9 16997.5 19858.9 [55] 17681.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]) 12522.5 13779.7 13855.6 14308.3 14384.5 14423.8 14733.8 14851.9 15160 15290.6 1 1 1 1 1 1 1 1 1 1 15638.6 15686.3 15692.1 15806.8 15903.8 16007.1 16016.5 16059.1 16141.4 16189.4 1 1 1 1 1 1 1 1 1 1 16192.9 16198.9 16260.7 16471.5 16554.2 16969.9 16997.5 17167.2 17183.6 17428.5 1 1 1 1 1 1 1 1 1 1 17681.2 18003.8 18147.5 18174.1 18329.6 18344.7 18374.4 18406.6 18466.5 18746.8 1 1 1 1 1 1 1 1 1 1 18908.9 18957.2 19009.5 19211.2 19301.4 19325.8 19554.2 19630 19711.6 19858.9 1 1 1 1 1 1 1 1 1 1 20056.9 20359.8 20515.2 20547.7 20605.5 1 1 1 1 1 > colnames(x) [1] "uitvoer" "invoer" "crisis" > colnames(x)[par1] [1] "uitvoer" > x[,par1] [1] 16198.9 16554.2 19554.2 15903.8 18003.8 18329.6 16260.7 14851.9 18174.1 [10] 18406.6 18466.5 16016.5 17428.5 17167.2 19630.0 17183.6 18344.7 19301.4 [19] 18147.5 16192.9 18374.4 20515.2 18957.2 16471.5 18746.8 19009.5 19211.2 [28] 20547.7 19325.8 20605.5 20056.9 16141.4 20359.8 19711.6 15638.6 14384.5 [37] 13855.6 14308.3 15290.6 14423.8 13779.7 15686.3 14733.8 12522.5 16189.4 [46] 16059.1 16007.1 15806.8 15160.0 15692.1 18908.9 16969.9 16997.5 19858.9 [55] 17681.2 > 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/10v261292345039.tab") + } + } > m Conditional inference tree with 4 terminal nodes Response: uitvoer Inputs: invoer, crisis Number of observations: 55 1) invoer <= 17180.4; criterion = 1, statistic = 47.242 2) invoer <= 15816.8; criterion = 1, statistic = 15.967 3)* weights = 11 2) invoer > 15816.8 4)* weights = 14 1) invoer > 17180.4 5) invoer <= 18203.6; criterion = 1, statistic = 19.183 6)* weights = 14 5) invoer > 18203.6 7)* weights = 16 > postscript(file="/var/www/html/freestat/rcomp/tmp/2bmkr1292345039.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/3bmkr1292345039.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 16198.9 16126.44 72.45714 2 16554.2 16126.44 427.75714 3 19554.2 19643.79 -89.58750 4 15903.8 16126.44 -222.64286 5 18003.8 17709.80 294.00000 6 18329.6 17709.80 619.80000 7 16260.7 14554.34 1706.36364 8 14851.9 14554.34 297.56364 9 18174.1 17709.80 464.30000 10 18406.6 17709.80 696.80000 11 18466.5 17709.80 756.70000 12 16016.5 16126.44 -109.94286 13 17428.5 17709.80 -281.30000 14 17167.2 16126.44 1040.75714 15 19630.0 19643.79 -13.78750 16 17183.6 16126.44 1057.15714 17 18344.7 17709.80 634.90000 18 19301.4 19643.79 -342.38750 19 18147.5 17709.80 437.70000 20 16192.9 16126.44 66.45714 21 18374.4 17709.80 664.60000 22 20515.2 19643.79 871.41250 23 18957.2 19643.79 -686.58750 24 16471.5 17709.80 -1238.30000 25 18746.8 19643.79 -896.98750 26 19009.5 19643.79 -634.28750 27 19211.2 19643.79 -432.58750 28 20547.7 19643.79 903.91250 29 19325.8 19643.79 -317.98750 30 20605.5 19643.79 961.71250 31 20056.9 19643.79 413.11250 32 16141.4 17709.80 -1568.40000 33 20359.8 19643.79 716.01250 34 19711.6 19643.79 67.81250 35 15638.6 16126.44 -487.84286 36 14384.5 14554.34 -169.83636 37 13855.6 14554.34 -698.73636 38 14308.3 14554.34 -246.03636 39 15290.6 14554.34 736.26364 40 14423.8 14554.34 -130.53636 41 13779.7 14554.34 -774.63636 42 15686.3 14554.34 1131.96364 43 14733.8 14554.34 179.46364 44 12522.5 14554.34 -2031.83636 45 16189.4 16126.44 62.95714 46 16059.1 16126.44 -67.34286 47 16007.1 16126.44 -119.34286 48 15806.8 16126.44 -319.64286 49 15160.0 16126.44 -966.44286 50 15692.1 16126.44 -434.34286 51 18908.9 19643.79 -734.88750 52 16969.9 17709.80 -739.90000 53 16997.5 17709.80 -712.30000 54 19858.9 19643.79 215.11250 55 17681.2 17709.80 -28.60000 > 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/4lv1c1292345039.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/57ez01292345039.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/6sfgn1292345039.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/7efwt1292345039.tab") + } > > try(system("convert tmp/2bmkr1292345039.ps tmp/2bmkr1292345039.png",intern=TRUE)) character(0) > try(system("convert tmp/3bmkr1292345039.ps tmp/3bmkr1292345039.png",intern=TRUE)) character(0) > try(system("convert tmp/4lv1c1292345039.ps tmp/4lv1c1292345039.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.696 0.739 3.846