R version 2.15.2 (2012-10-26) -- "Trick or Treat" Copyright (C) 2012 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i686-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(1966 + ,1 + ,41 + ,1966 + ,2 + ,39 + ,1966 + ,3 + ,50 + ,1966 + ,4 + ,40 + ,1966 + ,5 + ,43 + ,1966 + ,6 + ,38 + ,1966 + ,7 + ,44 + ,1966 + ,8 + ,35 + ,1966 + ,9 + ,39 + ,1966 + ,10 + ,35 + ,1966 + ,11 + ,29 + ,1966 + ,12 + ,49 + ,1967 + ,1 + ,50 + ,1967 + ,2 + ,59 + ,1967 + ,3 + ,63 + ,1967 + ,4 + ,32 + ,1967 + ,5 + ,39 + ,1967 + ,6 + ,47 + ,1967 + ,7 + ,53 + ,1967 + ,8 + ,60 + ,1967 + ,9 + ,57 + ,1967 + ,10 + ,52 + ,1967 + ,11 + ,70 + ,1967 + ,12 + ,90 + ,1968 + ,1 + ,74 + ,1968 + ,2 + ,62 + ,1968 + ,3 + ,55 + ,1968 + ,4 + ,84 + ,1968 + ,5 + ,94 + ,1968 + ,6 + ,70 + ,1968 + ,7 + ,108 + ,1968 + ,8 + ,139 + ,1968 + ,9 + ,120 + ,1968 + ,10 + ,97 + ,1968 + ,11 + ,126 + ,1968 + ,12 + ,149 + ,1969 + ,1 + ,158 + ,1969 + ,2 + ,124 + ,1969 + ,3 + ,140 + ,1969 + ,4 + ,109 + ,1969 + ,5 + ,114 + ,1969 + ,6 + ,77 + ,1969 + ,7 + ,120 + ,1969 + ,8 + ,133 + ,1969 + ,9 + ,110 + ,1969 + ,10 + ,92 + ,1969 + ,11 + ,97 + ,1969 + ,12 + ,78 + ,1970 + ,1 + ,99 + ,1970 + ,2 + ,107 + ,1970 + ,3 + ,112 + ,1970 + ,4 + ,90 + ,1970 + ,5 + ,98 + ,1970 + ,6 + ,125 + ,1970 + ,7 + ,155 + ,1970 + ,8 + ,190 + ,1970 + ,9 + ,236 + ,1970 + ,10 + ,189 + ,1970 + ,11 + ,174 + ,1970 + ,12 + ,178 + ,1971 + ,1 + ,136 + ,1971 + ,2 + ,161 + ,1971 + ,3 + ,171 + ,1971 + ,4 + ,149 + ,1971 + ,5 + ,184 + ,1971 + ,6 + ,155 + ,1971 + ,7 + ,276 + ,1971 + ,8 + ,224 + ,1971 + ,9 + ,213 + ,1971 + ,10 + ,279 + ,1971 + ,11 + ,268 + ,1971 + ,12 + ,287 + ,1972 + ,1 + ,238 + ,1972 + ,2 + ,213 + ,1972 + ,3 + ,257 + ,1972 + ,4 + ,293 + ,1972 + ,5 + ,212 + ,1972 + ,6 + ,246 + ,1972 + ,7 + ,353 + ,1972 + ,8 + ,339 + ,1972 + ,9 + ,308 + ,1972 + ,10 + ,247 + ,1972 + ,11 + ,257 + ,1972 + ,12 + ,322 + ,1973 + ,1 + ,298 + ,1973 + ,2 + ,273 + ,1973 + ,3 + ,312 + ,1973 + ,4 + ,249 + ,1973 + ,5 + ,286 + ,1973 + ,6 + ,279 + ,1973 + ,7 + ,309 + ,1973 + ,8 + ,401 + ,1973 + ,9 + ,309 + ,1973 + ,10 + ,328 + ,1973 + ,11 + ,353 + ,1973 + ,12 + ,354 + ,1974 + ,1 + ,327 + ,1974 + ,2 + ,324 + ,1974 + ,3 + ,285 + ,1974 + ,4 + ,243 + ,1974 + ,5 + ,241 + ,1974 + ,6 + ,287 + ,1974 + ,7 + ,355 + ,1974 + ,8 + ,460 + ,1974 + ,9 + ,364 + ,1974 + ,10 + ,487 + ,1974 + ,11 + ,452 + ,1974 + ,12 + ,391 + ,1975 + ,1 + ,500 + ,1975 + ,2 + ,451 + ,1975 + ,3 + ,375 + ,1975 + ,4 + ,372 + ,1975 + ,5 + ,302 + ,1975 + ,6 + ,316 + ,1975 + ,7 + ,398 + ,1975 + ,8 + ,394 + ,1975 + ,9 + ,431 + ,1975 + ,10 + ,431) + ,dim=c(3 + ,118) + ,dimnames=list(c('Year' + ,'month' + ,'robberies') + ,1:118)) > y <- array(NA,dim=c(3,118),dimnames=list(c('Year','month','robberies'),1:118)) > 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 = 'none' > par2 = 'none' > par1 = '3' > 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, 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) Hmisc library by Frank E Harrell Jr Type library(help='Hmisc'), ?Overview, or ?Hmisc.Overview') to see overall documentation. NOTE:Hmisc no longer redefines [.factor to drop unused levels when subsetting. To get the old behavior of Hmisc type dropUnusedLevels(). 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) Warning message: NAs introduced by coercion > 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] "robberies" > x[,par1] [1] 41 39 50 40 43 38 44 35 39 35 29 49 50 59 63 32 39 47 [19] 53 60 57 52 70 90 74 62 55 84 94 70 108 139 120 97 126 149 [37] 158 124 140 109 114 77 120 133 110 92 97 78 99 107 112 90 98 125 [55] 155 190 236 189 174 178 136 161 171 149 184 155 276 224 213 279 268 287 [73] 238 213 257 293 212 246 353 339 308 247 257 322 298 273 312 249 286 279 [91] 309 401 309 328 353 354 327 324 285 243 241 287 355 460 364 487 452 391 [109] 500 451 375 372 302 316 398 394 431 431 > 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]) 29 32 35 38 39 40 41 43 44 47 49 50 52 53 55 57 59 60 62 63 1 1 2 1 3 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 70 74 77 78 84 90 92 94 97 98 99 107 108 109 110 112 114 120 124 125 2 1 1 1 1 2 1 1 2 1 1 1 1 1 1 1 1 2 1 1 126 133 136 139 140 149 155 158 161 171 174 178 184 189 190 212 213 224 236 238 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 2 1 1 1 241 243 246 247 249 257 268 273 276 279 285 286 287 293 298 302 308 309 312 316 1 1 1 1 1 2 1 1 1 2 1 1 2 1 1 1 1 2 1 1 322 324 327 328 339 353 354 355 364 372 375 391 394 398 401 431 451 452 460 487 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 1 500 1 > colnames(x) [1] "Year" "month" "robberies" > colnames(x)[par1] [1] "robberies" > x[,par1] [1] 41 39 50 40 43 38 44 35 39 35 29 49 50 59 63 32 39 47 [19] 53 60 57 52 70 90 74 62 55 84 94 70 108 139 120 97 126 149 [37] 158 124 140 109 114 77 120 133 110 92 97 78 99 107 112 90 98 125 [55] 155 190 236 189 174 178 136 161 171 149 184 155 276 224 213 279 268 287 [73] 238 213 257 293 212 246 353 339 308 247 257 322 298 273 312 249 286 279 [91] 309 401 309 328 353 354 327 324 285 243 241 287 355 460 364 487 452 391 [109] 500 451 375 372 302 316 398 394 431 431 > if (par2 == 'none') { + m <- ctree(as.formula(paste(colnames(x)[par1],' ~ .',sep='')),data = x) + } > > #Note: the /var/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/fisher/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/fisher/rcomp/tmp/1cjvq1354892553.tab") + } + } > m Conditional inference tree with 8 terminal nodes Response: robberies Inputs: Year, month Number of observations: 118 1) Year <= 1971; criterion = 1, statistic = 100.673 2) Year <= 1969; criterion = 1, statistic = 50.304 3) Year <= 1967; criterion = 1, statistic = 30.687 4) Year <= 1966; criterion = 0.991, statistic = 8.005 5)* weights = 12 4) Year > 1966 6)* weights = 12 3) Year > 1967 7)* weights = 24 2) Year > 1969 8) month <= 6; criterion = 0.999, statistic = 11.773 9)* weights = 12 8) month > 6 10)* weights = 12 1) Year > 1971 11) Year <= 1973; criterion = 1, statistic = 16.877 12) month <= 6; criterion = 0.969, statistic = 5.836 13)* weights = 12 12) month > 6 14)* weights = 12 11) Year > 1973 15)* weights = 22 > postscript(file="/var/fisher/rcomp/tmp/210ax1354892553.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/fisher/rcomp/tmp/3ab751354892553.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 41 40.16667 0.83333333 2 39 40.16667 -1.16666667 3 50 40.16667 9.83333333 4 40 40.16667 -0.16666667 5 43 40.16667 2.83333333 6 38 40.16667 -2.16666667 7 44 40.16667 3.83333333 8 35 40.16667 -5.16666667 9 39 40.16667 -1.16666667 10 35 40.16667 -5.16666667 11 29 40.16667 -11.16666667 12 49 40.16667 8.83333333 13 50 56.00000 -6.00000000 14 59 56.00000 3.00000000 15 63 56.00000 7.00000000 16 32 56.00000 -24.00000000 17 39 56.00000 -17.00000000 18 47 56.00000 -9.00000000 19 53 56.00000 -3.00000000 20 60 56.00000 4.00000000 21 57 56.00000 1.00000000 22 52 56.00000 -4.00000000 23 70 56.00000 14.00000000 24 90 56.00000 34.00000000 25 74 105.41667 -31.41666667 26 62 105.41667 -43.41666667 27 55 105.41667 -50.41666667 28 84 105.41667 -21.41666667 29 94 105.41667 -11.41666667 30 70 105.41667 -35.41666667 31 108 105.41667 2.58333333 32 139 105.41667 33.58333333 33 120 105.41667 14.58333333 34 97 105.41667 -8.41666667 35 126 105.41667 20.58333333 36 149 105.41667 43.58333333 37 158 105.41667 52.58333333 38 124 105.41667 18.58333333 39 140 105.41667 34.58333333 40 109 105.41667 3.58333333 41 114 105.41667 8.58333333 42 77 105.41667 -28.41666667 43 120 105.41667 14.58333333 44 133 105.41667 27.58333333 45 110 105.41667 4.58333333 46 92 105.41667 -13.41666667 47 97 105.41667 -8.41666667 48 78 105.41667 -27.41666667 49 99 132.25000 -33.25000000 50 107 132.25000 -25.25000000 51 112 132.25000 -20.25000000 52 90 132.25000 -42.25000000 53 98 132.25000 -34.25000000 54 125 132.25000 -7.25000000 55 155 222.41667 -67.41666667 56 190 222.41667 -32.41666667 57 236 222.41667 13.58333333 58 189 222.41667 -33.41666667 59 174 222.41667 -48.41666667 60 178 222.41667 -44.41666667 61 136 132.25000 3.75000000 62 161 132.25000 28.75000000 63 171 132.25000 38.75000000 64 149 132.25000 16.75000000 65 184 132.25000 51.75000000 66 155 132.25000 22.75000000 67 276 222.41667 53.58333333 68 224 222.41667 1.58333333 69 213 222.41667 -9.41666667 70 279 222.41667 56.58333333 71 268 222.41667 45.58333333 72 287 222.41667 64.58333333 73 238 263.00000 -25.00000000 74 213 263.00000 -50.00000000 75 257 263.00000 -6.00000000 76 293 263.00000 30.00000000 77 212 263.00000 -51.00000000 78 246 263.00000 -17.00000000 79 353 323.33333 29.66666667 80 339 323.33333 15.66666667 81 308 323.33333 -15.33333333 82 247 323.33333 -76.33333333 83 257 323.33333 -66.33333333 84 322 323.33333 -1.33333333 85 298 263.00000 35.00000000 86 273 263.00000 10.00000000 87 312 263.00000 49.00000000 88 249 263.00000 -14.00000000 89 286 263.00000 23.00000000 90 279 263.00000 16.00000000 91 309 323.33333 -14.33333333 92 401 323.33333 77.66666667 93 309 323.33333 -14.33333333 94 328 323.33333 4.66666667 95 353 323.33333 29.66666667 96 354 323.33333 30.66666667 97 327 372.09091 -45.09090909 98 324 372.09091 -48.09090909 99 285 372.09091 -87.09090909 100 243 372.09091 -129.09090909 101 241 372.09091 -131.09090909 102 287 372.09091 -85.09090909 103 355 372.09091 -17.09090909 104 460 372.09091 87.90909091 105 364 372.09091 -8.09090909 106 487 372.09091 114.90909091 107 452 372.09091 79.90909091 108 391 372.09091 18.90909091 109 500 372.09091 127.90909091 110 451 372.09091 78.90909091 111 375 372.09091 2.90909091 112 372 372.09091 -0.09090909 113 302 372.09091 -70.09090909 114 316 372.09091 -56.09090909 115 398 372.09091 25.90909091 116 394 372.09091 21.90909091 117 431 372.09091 58.90909091 118 431 372.09091 58.90909091 > if (par2 != 'none') { + print(cbind(as.factor(x[,par1]),predict(m))) + myt <- table(as.factor(x[,par1]),predict(m)) + print(myt) + } > postscript(file="/var/fisher/rcomp/tmp/456c61354892553.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/fisher/rcomp/tmp/5c6qi1354892553.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/fisher/rcomp/tmp/6v2081354892553.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/fisher/rcomp/tmp/7amkc1354892553.tab") + } > > try(system("convert tmp/210ax1354892553.ps tmp/210ax1354892553.png",intern=TRUE)) character(0) > try(system("convert tmp/3ab751354892553.ps tmp/3ab751354892553.png",intern=TRUE)) character(0) > try(system("convert tmp/456c61354892553.ps tmp/456c61354892553.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 4.597 0.571 5.155