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(356,182,89,386,213,97,444,227,154,387,209,81,327,219,110,448,221,116,225,114,73,182,97,73,460,205,174,411,215,103,342,224,130,361,189,91,377,182,136,331,201,106,428,198,136,340,173,122,352,238,131,461,258,135,221,122,75,198,101,68,422,259,143,329,243,115,320,188,93,375,173,128,364,224,152,351,215,125,380,196,107,319,159,116,322,187,220,386,208,137,221,131,34,187,93,51,344,210,153,342,228,145,365,176,116,313,195,145,356,188,98,337,188,118,389,190,139,326,188,140,343,176,113,357,225,149,220,93,79,218,79,47,391,235,166,425,247,180,332,195,122,298,197,134,360,211,114,336,156,125,325,209,181,393,180,142,301,185,143,426,303,187,265,129,137,210,85,62,429,249,239,440,231,157,357,212,139,431,240,187,442,234,99,442,217,146,544,287,175,420,221,148,396,208,130,482,241,183,261,156,115,211,96,80,448,320,223,468,242,131,464,227,201,425,200,157),dim=c(3,72),dimnames=list(c('Vlaanderen','Wallonië','Brussel'),1:72)) > y <- array(NA,dim=c(3,72),dimnames=list(c('Vlaanderen','Wallonië','Brussel'),1:72)) > 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] "Vlaanderen" > x[,par1] [1] 356 386 444 387 327 448 225 182 460 411 342 361 377 331 428 340 352 461 221 [20] 198 422 329 320 375 364 351 380 319 322 386 221 187 344 342 365 313 356 337 [39] 389 326 343 357 220 218 391 425 332 298 360 336 325 393 301 426 265 210 429 [58] 440 357 431 442 442 544 420 396 482 261 211 448 468 464 425 > 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]) 182 187 198 210 211 218 220 221 225 261 265 298 301 313 319 320 322 325 326 327 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 329 331 332 336 337 340 342 343 344 351 352 356 357 360 361 364 365 375 377 380 1 1 1 1 1 1 2 1 1 1 1 2 2 1 1 1 1 1 1 1 386 387 389 391 393 396 411 420 422 425 426 428 429 431 440 442 444 448 460 461 2 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 2 1 1 464 468 482 544 1 1 1 1 > colnames(x) [1] "Vlaanderen" "Walloni.." "Brussel" > colnames(x)[par1] [1] "Vlaanderen" > x[,par1] [1] 356 386 444 387 327 448 225 182 460 411 342 361 377 331 428 340 352 461 221 [20] 198 422 329 320 375 364 351 380 319 322 386 221 187 344 342 365 313 356 337 [39] 389 326 343 357 220 218 391 425 332 298 360 336 325 393 301 426 265 210 429 [58] 440 357 431 442 442 544 420 396 482 261 211 448 468 464 425 > 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/1w0z81291987304.tab") + } + } > m Conditional inference tree with 3 terminal nodes Response: Vlaanderen Inputs: Walloni.., Brussel Number of observations: 72 1) Walloni.. <= 156; criterion = 1, statistic = 51.835 2)* weights = 13 1) Walloni.. > 156 3) Walloni.. <= 225; criterion = 1, statistic = 21.429 4)* weights = 41 3) Walloni.. > 225 5)* weights = 18 > postscript(file="/var/www/html/rcomp/tmp/2w0z81291987304.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/3w0z81291987304.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 356 364.6341 -8.6341463 2 386 364.6341 21.3658537 3 444 430.0000 14.0000000 4 387 364.6341 22.3658537 5 327 364.6341 -37.6341463 6 448 364.6341 83.3658537 7 225 227.3077 -2.3076923 8 182 227.3077 -45.3076923 9 460 364.6341 95.3658537 10 411 364.6341 46.3658537 11 342 364.6341 -22.6341463 12 361 364.6341 -3.6341463 13 377 364.6341 12.3658537 14 331 364.6341 -33.6341463 15 428 364.6341 63.3658537 16 340 364.6341 -24.6341463 17 352 430.0000 -78.0000000 18 461 430.0000 31.0000000 19 221 227.3077 -6.3076923 20 198 227.3077 -29.3076923 21 422 430.0000 -8.0000000 22 329 430.0000 -101.0000000 23 320 364.6341 -44.6341463 24 375 364.6341 10.3658537 25 364 364.6341 -0.6341463 26 351 364.6341 -13.6341463 27 380 364.6341 15.3658537 28 319 364.6341 -45.6341463 29 322 364.6341 -42.6341463 30 386 364.6341 21.3658537 31 221 227.3077 -6.3076923 32 187 227.3077 -40.3076923 33 344 364.6341 -20.6341463 34 342 430.0000 -88.0000000 35 365 364.6341 0.3658537 36 313 364.6341 -51.6341463 37 356 364.6341 -8.6341463 38 337 364.6341 -27.6341463 39 389 364.6341 24.3658537 40 326 364.6341 -38.6341463 41 343 364.6341 -21.6341463 42 357 364.6341 -7.6341463 43 220 227.3077 -7.3076923 44 218 227.3077 -9.3076923 45 391 430.0000 -39.0000000 46 425 430.0000 -5.0000000 47 332 364.6341 -32.6341463 48 298 364.6341 -66.6341463 49 360 364.6341 -4.6341463 50 336 227.3077 108.6923077 51 325 364.6341 -39.6341463 52 393 364.6341 28.3658537 53 301 364.6341 -63.6341463 54 426 430.0000 -4.0000000 55 265 227.3077 37.6923077 56 210 227.3077 -17.3076923 57 429 430.0000 -1.0000000 58 440 430.0000 10.0000000 59 357 364.6341 -7.6341463 60 431 430.0000 1.0000000 61 442 430.0000 12.0000000 62 442 364.6341 77.3658537 63 544 430.0000 114.0000000 64 420 364.6341 55.3658537 65 396 364.6341 31.3658537 66 482 430.0000 52.0000000 67 261 227.3077 33.6923077 68 211 227.3077 -16.3076923 69 448 430.0000 18.0000000 70 468 430.0000 38.0000000 71 464 430.0000 34.0000000 72 425 364.6341 60.3658537 > 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/4hjyw1291987304.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/521ek1291987304.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/6dsen1291987304.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/7gbct1291987304.tab") + } > > try(system("convert tmp/2w0z81291987304.ps tmp/2w0z81291987304.png",intern=TRUE)) character(0) > try(system("convert tmp/3w0z81291987304.ps tmp/3w0z81291987304.png",intern=TRUE)) character(0) > try(system("convert tmp/4hjyw1291987304.ps tmp/4hjyw1291987304.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.191 0.538 5.390