R version 2.9.0 (2009-04-17) Copyright (C) 2009 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. 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(7361 + ,493 + ,797 + ,48 + ,1.5 + ,105.0 + ,508643 + ,7391 + ,514 + ,840 + ,49 + ,1.6 + ,104.0 + ,527568 + ,7420 + ,522 + ,988 + ,59 + ,1.8 + ,109.8 + ,520008 + ,7406 + ,490 + ,819 + ,56 + ,1.5 + ,98.6 + ,498484 + ,7439 + ,484 + ,831 + ,47 + ,1.3 + ,93.5 + ,523917 + ,7512 + ,506 + ,904 + ,56 + ,1.6 + ,98.2 + ,553522 + ,7579 + ,501 + ,814 + ,50 + ,1.6 + ,88.0 + ,558901 + ,7520 + ,462 + ,798 + ,54 + ,1.8 + ,85.3 + ,548933 + ,7453 + ,465 + ,828 + ,79 + ,1.8 + ,96.8 + ,567013 + ,7462 + ,454 + ,789 + ,50 + ,1.6 + ,98.8 + ,551085 + ,7472 + ,464 + ,930 + ,54 + ,1.8 + ,110.3 + ,588245 + ,7443 + ,427 + ,744 + ,56 + ,2 + ,111.6 + ,605010 + ,7439 + ,460 + ,832 + ,50 + ,1.3 + ,111.2 + ,631572 + ,7460 + ,473 + ,826 + ,46 + ,1.1 + ,106.9 + ,639180 + ,7482 + ,465 + ,907 + ,47 + ,1 + ,117.6 + ,653847 + ,7442 + ,422 + ,776 + ,43 + ,1.2 + ,97.0 + ,657073 + ,7454 + ,415 + ,835 + ,52 + ,1.2 + ,97.3 + ,626291 + ,7536 + ,413 + ,715 + ,48 + ,1.3 + ,98.4 + ,625616 + ,7616 + ,420 + ,729 + ,36 + ,1.3 + ,87.6 + ,633352 + ,7548 + ,363 + ,733 + ,41 + ,1.4 + ,87.4 + ,672820 + ,7507 + ,376 + ,736 + ,34 + ,1.1 + ,94.7 + ,691369 + ,7515 + ,380 + ,712 + ,37 + ,0.9 + ,101.5 + ,702595 + ,7549 + ,384 + ,711 + ,37 + ,1 + ,110.4 + ,692241 + ,7540 + ,346 + ,667 + ,34 + ,1.1 + ,108.4 + ,718722 + ,7525 + ,389 + ,799 + ,55 + ,1.4 + ,109.7 + ,732297 + ,7575 + ,407 + ,661 + ,37 + ,1.5 + ,105.2 + ,721798 + ,7621 + ,393 + ,692 + ,27 + ,1.8 + ,111.1 + ,766192 + ,7589 + ,346 + ,649 + ,38 + ,1.8 + ,96.2 + ,788456 + ,7606 + ,348 + ,729 + ,43 + ,1.8 + ,97.3 + ,806132 + ,7722 + ,353 + ,622 + ,26 + ,1.7 + ,98.9 + ,813944 + ,7788 + ,364 + ,671 + ,32 + ,1.5 + ,91.7 + ,788025 + ,7735 + ,305 + ,635 + ,29 + ,1.1 + ,90.9 + ,765985 + ,7654 + ,307 + ,648 + ,41 + ,1.3 + ,98.8 + ,702684 + ,7678 + ,312 + ,745 + ,55 + ,1.6 + ,111.5 + ,730159 + ,7688 + ,312 + ,624 + ,50 + ,1.9 + ,119.0 + ,678942 + ,7653 + ,286 + ,477 + ,30 + ,1.9 + ,115.3 + ,672527 + ,7688 + ,324 + ,710 + ,35 + ,2 + ,116.3 + ,594783 + ,7734 + ,336 + ,515 + ,29 + ,2.2 + ,113.6 + ,594575 + ,7754 + ,327 + ,461 + ,22 + ,2.2 + ,115.1 + ,576299 + ,7760 + ,302 + ,590 + ,39 + ,2 + ,109.7 + ,530770 + ,7770 + ,299 + ,415 + ,24 + ,2.3 + ,97.6 + ,524491 + ,7867 + ,311 + ,554 + ,38 + ,2.6 + ,100.8 + ,456590 + ,7938 + ,315 + ,585 + ,30 + ,3.2 + ,94.0 + ,428448 + ,7860 + ,264 + ,513 + ,31 + ,3.2 + ,87.2 + ,444937 + ,7793 + ,278 + ,591 + ,39 + ,3.1 + ,102.9 + ,372206 + ,7829 + ,278 + ,561 + ,33 + ,2.8 + ,111.3 + ,317272 + ,7828 + ,287 + ,684 + ,57 + ,2.3 + ,106.6 + ,297604 + ,7789 + ,279 + ,668 + ,49 + ,1.9 + ,108.9 + ,288561 + ,7820 + ,324 + ,795 + ,74 + ,1.9 + ,108.2 + ,289287 + ,7850 + ,354 + ,776 + ,74 + ,2 + ,100.2 + ,258923 + ,7860 + ,354 + ,1043 + ,115 + ,2 + ,104.0 + ,255493 + ,7836 + ,360 + ,964 + ,67 + ,1.8 + ,90.0 + ,277992 + ,7844 + ,363 + ,762 + ,51 + ,1.6 + ,87.4 + ,295474 + ,7915 + ,385 + ,1030 + ,114 + ,1.4 + ,91.9 + ,291680 + ,7971 + ,412 + ,939 + ,70 + ,0.2 + ,89.3 + ,318736 + ,7890 + ,370 + ,779 + ,73 + ,0.3 + ,81.3 + ,338463 + ,7807 + ,389 + ,918 + ,77 + ,0.4 + ,94.9 + ,351963 + ,7797 + ,395 + ,839 + ,67 + ,0.7 + ,102.6 + ,347240 + ,7788 + ,417 + ,874 + ,60 + ,1 + ,107.2 + ,347081 + ,7779 + ,404 + ,840 + ,73 + ,1.1 + ,114.0 + ,383486) + ,dim=c(7 + ,60) + ,dimnames=list(c('Beroepsbevolking' + ,'Werkloosheid' + ,'Faillissementen' + ,'Faillissementennijverheid' + ,'Inflatie' + ,'Nijverheid' + ,'Beurswaarde ') + ,1:60)) > y <- array(NA,dim=c(7,60),dimnames=list(c('Beroepsbevolking','Werkloosheid','Faillissementen','Faillissementennijverheid','Inflatie','Nijverheid','Beurswaarde '),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 = '3' > 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] "Beroepsbevolking" > x[,par1] [1] 7361 7391 7420 7406 7439 7512 7579 7520 7453 7462 7472 7443 7439 7460 7482 [16] 7442 7454 7536 7616 7548 7507 7515 7549 7540 7525 7575 7621 7589 7606 7722 [31] 7788 7735 7654 7678 7688 7653 7688 7734 7754 7760 7770 7867 7938 7860 7793 [46] 7829 7828 7789 7820 7850 7860 7836 7844 7915 7971 7890 7807 7797 7788 7779 > 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]) 7361 7391 7406 7420 7439 7442 7443 7453 7454 7460 7462 7472 7482 7507 7512 7515 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 7520 7525 7536 7540 7548 7549 7575 7579 7589 7606 7616 7621 7653 7654 7678 7688 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 7722 7734 7735 7754 7760 7770 7779 7788 7789 7793 7797 7807 7820 7828 7829 7836 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 7844 7850 7860 7867 7890 7915 7938 7971 1 1 2 1 1 1 1 1 > colnames(x) [1] "Beroepsbevolking" "Werkloosheid" [3] "Faillissementen" "Faillissementennijverheid" [5] "Inflatie" "Nijverheid" [7] "Beurswaarde." > colnames(x)[par1] [1] "Beroepsbevolking" > x[,par1] [1] 7361 7391 7420 7406 7439 7512 7579 7520 7453 7462 7472 7443 7439 7460 7482 [16] 7442 7454 7536 7616 7548 7507 7515 7549 7540 7525 7575 7621 7589 7606 7722 [31] 7788 7735 7654 7678 7688 7653 7688 7734 7754 7760 7770 7867 7938 7860 7793 [46] 7829 7828 7789 7820 7850 7860 7836 7844 7915 7971 7890 7807 7797 7788 7779 > if (par2 == 'none') { + m <- ctree(as.formula(paste(colnames(x)[par1],' ~ .',sep='')),data = x) + } > > #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/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/rcomp/tmp/19d261293378483.tab") + } + } > m Conditional inference tree with 4 terminal nodes Response: Beroepsbevolking Inputs: Werkloosheid, Faillissementen, Faillissementennijverheid, Inflatie, Nijverheid, Beurswaarde. Number of observations: 60 1) Werkloosheid <= 412; criterion = 1, statistic = 28.419 2) Beurswaarde. <= 594575; criterion = 1, statistic = 24.596 3) Nijverheid <= 94; criterion = 0.991, statistic = 9.988 4)* weights = 7 3) Nijverheid > 94 5)* weights = 15 2) Beurswaarde. > 594575 6)* weights = 18 1) Werkloosheid > 412 7)* weights = 20 > postscript(file="/var/www/html/rcomp/tmp/29d261293378483.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/rcomp/tmp/324k91293378483.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 7361 7483.750 -122.7500000 2 7391 7483.750 -92.7500000 3 7420 7483.750 -63.7500000 4 7406 7483.750 -77.7500000 5 7439 7483.750 -44.7500000 6 7512 7483.750 28.2500000 7 7579 7483.750 95.2500000 8 7520 7483.750 36.2500000 9 7453 7483.750 -30.7500000 10 7462 7483.750 -21.7500000 11 7472 7483.750 -11.7500000 12 7443 7483.750 -40.7500000 13 7439 7483.750 -44.7500000 14 7460 7483.750 -23.7500000 15 7482 7483.750 -1.7500000 16 7442 7483.750 -41.7500000 17 7454 7483.750 -29.7500000 18 7536 7483.750 52.2500000 19 7616 7483.750 132.2500000 20 7548 7621.167 -73.1666667 21 7507 7621.167 -114.1666667 22 7515 7621.167 -106.1666667 23 7549 7621.167 -72.1666667 24 7540 7621.167 -81.1666667 25 7525 7621.167 -96.1666667 26 7575 7621.167 -46.1666667 27 7621 7621.167 -0.1666667 28 7589 7621.167 -32.1666667 29 7606 7621.167 -15.1666667 30 7722 7621.167 100.8333333 31 7788 7621.167 166.8333333 32 7735 7621.167 113.8333333 33 7654 7621.167 32.8333333 34 7678 7621.167 56.8333333 35 7688 7621.167 66.8333333 36 7653 7621.167 31.8333333 37 7688 7621.167 66.8333333 38 7734 7802.467 -68.4666667 39 7754 7802.467 -48.4666667 40 7760 7802.467 -42.4666667 41 7770 7802.467 -32.4666667 42 7867 7802.467 64.5333333 43 7938 7893.429 44.5714286 44 7860 7893.429 -33.4285714 45 7793 7802.467 -9.4666667 46 7829 7802.467 26.5333333 47 7828 7802.467 25.5333333 48 7789 7802.467 -13.4666667 49 7820 7802.467 17.5333333 50 7850 7802.467 47.5333333 51 7860 7802.467 57.5333333 52 7836 7893.429 -57.4285714 53 7844 7893.429 -49.4285714 54 7915 7893.429 21.5714286 55 7971 7893.429 77.5714286 56 7890 7893.429 -3.4285714 57 7807 7802.467 4.5333333 58 7797 7802.467 -5.4666667 59 7788 7483.750 304.2500000 60 7779 7802.467 -23.4666667 > 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/rcomp/tmp/4ud1c1293378483.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/rcomp/tmp/5yez01293378483.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/rcomp/tmp/6jxy61293378483.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/rcomp/tmp/75fet1293378483.tab") + } > > try(system("convert tmp/29d261293378483.ps tmp/29d261293378483.png",intern=TRUE)) character(0) > try(system("convert tmp/324k91293378483.ps tmp/324k91293378483.png",intern=TRUE)) character(0) > try(system("convert tmp/4ud1c1293378483.ps tmp/4ud1c1293378483.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.258 0.580 5.127