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(97.06 + ,21454 + ,631923 + ,130678 + ,97.73 + ,23899 + ,654294 + ,120877 + ,98 + ,24939 + ,671833 + ,137114 + ,97.76 + ,23580 + ,586840 + ,134406 + ,97.48 + ,24562 + ,600969 + ,120262 + ,97.77 + ,24696 + ,625568 + ,130846 + ,97.96 + ,23785 + ,558110 + ,120343 + ,98.22 + ,23812 + ,630577 + ,98881 + ,98.51 + ,21917 + ,628654 + ,115678 + ,98.19 + ,19713 + ,603184 + ,120796 + ,98.37 + ,19282 + ,656255 + ,94261 + ,98.31 + ,18788 + ,600730 + ,89151 + ,98.6 + ,21453 + ,670326 + ,119880 + ,98.96 + ,24482 + ,678423 + ,131468 + ,99.11 + ,27474 + ,641502 + ,155089 + ,99.64 + ,27264 + ,625311 + ,149581 + ,100.02 + ,27349 + ,628177 + ,122788 + ,99.98 + ,30632 + ,589767 + ,143900 + ,100.32 + ,29429 + ,582471 + ,112115 + ,100.44 + ,30084 + ,636248 + ,109600 + ,100.51 + ,26290 + ,599885 + ,117446 + ,101 + ,24379 + ,621694 + ,118456 + ,100.88 + ,23335 + ,637406 + ,101901 + ,100.55 + ,21346 + ,595994 + ,89940 + ,100.82 + ,21106 + ,696308 + ,129143 + ,101.5 + ,24514 + ,674201 + ,126102 + ,102.15 + ,28353 + ,648861 + ,143048 + ,102.39 + ,30805 + ,649605 + ,142258 + ,102.54 + ,31348 + ,672392 + ,131011 + ,102.85 + ,34556 + ,598396 + ,146471 + ,103.47 + ,33855 + ,613177 + ,114073 + ,103.56 + ,34787 + ,638104 + ,114642 + ,103.69 + ,32529 + ,615632 + ,118226 + ,103.49 + ,29998 + ,634465 + ,111338 + ,103.47 + ,29257 + ,638686 + ,108701 + ,103.45 + ,28155 + ,604243 + ,80512 + ,103.48 + ,30466 + ,706669 + ,146865 + ,103.93 + ,35704 + ,677185 + ,137179 + ,103.89 + ,39327 + ,644328 + ,166536 + ,104.4 + ,39351 + ,664825 + ,137070 + ,104.79 + ,42234 + ,605707 + ,127090 + ,104.77 + ,43630 + ,600136 + ,139966 + ,105.13 + ,43722 + ,612166 + ,122243 + ,105.26 + ,43121 + ,599659 + ,109097 + ,104.96 + ,37985 + ,634210 + ,116591 + ,104.75 + ,37135 + ,618234 + ,111964 + ,105.01 + ,34646 + ,613576 + ,109754 + ,105.15 + ,33026 + ,627200 + ,77609 + ,105.2 + ,35087 + ,668973 + ,138445 + ,105.77 + ,38846 + ,651479 + ,127901 + ,105.78 + ,42013 + ,619661 + ,156615 + ,106.26 + ,43908 + ,644260 + ,133264 + ,106.13 + ,42868 + ,579936 + ,143521 + ,106.12 + ,44423 + ,601752 + ,152139 + ,106.57 + ,44167 + ,595376 + ,131523 + ,106.44 + ,43636 + ,588902 + ,113925 + ,106.54 + ,44382 + ,634341 + ,86495 + ,107.1 + ,42142 + ,594305 + ,127877 + ,108.1 + ,43452 + ,606200 + ,107017 + ,108.4 + ,36912 + ,610926 + ,78716 + ,108.84 + ,42413 + ,633685 + ,138278 + ,109.62 + ,45344 + ,639696 + ,144238 + ,110.42 + ,44873 + ,659451 + ,143679 + ,110.67 + ,47510 + ,593248 + ,159932 + ,111.66 + ,49554 + ,606677 + ,136781 + ,112.28 + ,47369 + ,599434 + ,148173 + ,112.87 + ,45998 + ,569578 + ,125673 + ,112.18 + ,48140 + ,629873 + ,105573 + ,112.36 + ,48441 + ,613438 + ,122405 + ,112.16 + ,44928 + ,604172 + ,128045 + ,111.49 + ,40454 + ,658328 + ,94467 + ,111.25 + ,38661 + ,612633 + ,85573 + ,111.36 + ,37246 + ,707372 + ,121501 + ,111.74 + ,36843 + ,739770 + ,125074 + ,111.1 + ,36424 + ,777535 + ,144979 + ,111.33 + ,37594 + ,685030 + ,142120 + ,111.25 + ,38144 + ,730234 + ,124213 + ,111.04 + ,38737 + ,714154 + ,144407 + ,110.97 + ,34560 + ,630872 + ,125170 + ,111.31 + ,36080 + ,719492 + ,109267 + ,111.02 + ,33508 + ,677023 + ,122354 + ,111.07 + ,35462 + ,679272 + ,122589 + ,111.36 + ,33374 + ,718317 + ,104982 + ,111.54 + ,32110 + ,645672 + ,90542) + ,dim=c(4 + ,84) + ,dimnames=list(c('CPI' + ,'vacatures' + ,'werklozen' + ,'inschrijvingen') + ,1:84)) > y <- array(NA,dim=c(4,84),dimnames=list(c('CPI','vacatures','werklozen','inschrijvingen'),1:84)) > 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 = '4' > par2 = 'none' > par1 = '2' > #'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] "vacatures" > x[,par1] [1] 21454 23899 24939 23580 24562 24696 23785 23812 21917 19713 19282 18788 [13] 21453 24482 27474 27264 27349 30632 29429 30084 26290 24379 23335 21346 [25] 21106 24514 28353 30805 31348 34556 33855 34787 32529 29998 29257 28155 [37] 30466 35704 39327 39351 42234 43630 43722 43121 37985 37135 34646 33026 [49] 35087 38846 42013 43908 42868 44423 44167 43636 44382 42142 43452 36912 [61] 42413 45344 44873 47510 49554 47369 45998 48140 48441 44928 40454 38661 [73] 37246 36843 36424 37594 38144 38737 34560 36080 33508 35462 33374 32110 > 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]) 18788 19282 19713 21106 21346 21453 21454 21917 23335 23580 23785 23812 23899 1 1 1 1 1 1 1 1 1 1 1 1 1 24379 24482 24514 24562 24696 24939 26290 27264 27349 27474 28155 28353 29257 1 1 1 1 1 1 1 1 1 1 1 1 1 29429 29998 30084 30466 30632 30805 31348 32110 32529 33026 33374 33508 33855 1 1 1 1 1 1 1 1 1 1 1 1 1 34556 34560 34646 34787 35087 35462 35704 36080 36424 36843 36912 37135 37246 1 1 1 1 1 1 1 1 1 1 1 1 1 37594 37985 38144 38661 38737 38846 39327 39351 40454 42013 42142 42234 42413 1 1 1 1 1 1 1 1 1 1 1 1 1 42868 43121 43452 43630 43636 43722 43908 44167 44382 44423 44873 44928 45344 1 1 1 1 1 1 1 1 1 1 1 1 1 45998 47369 47510 48140 48441 49554 1 1 1 1 1 1 > colnames(x) [1] "CPI" "vacatures" "werklozen" "inschrijvingen" > colnames(x)[par1] [1] "vacatures" > x[,par1] [1] 21454 23899 24939 23580 24562 24696 23785 23812 21917 19713 19282 18788 [13] 21453 24482 27474 27264 27349 30632 29429 30084 26290 24379 23335 21346 [25] 21106 24514 28353 30805 31348 34556 33855 34787 32529 29998 29257 28155 [37] 30466 35704 39327 39351 42234 43630 43722 43121 37985 37135 34646 33026 [49] 35087 38846 42013 43908 42868 44423 44167 43636 44382 42142 43452 36912 [61] 42413 45344 44873 47510 49554 47369 45998 48140 48441 44928 40454 38661 [73] 37246 36843 36424 37594 38144 38737 34560 36080 33508 35462 33374 32110 > 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/13s7b1293006243.tab") + } + } > m Conditional inference tree with 4 terminal nodes Response: vacatures Inputs: CPI, werklozen, inschrijvingen Number of observations: 84 1) CPI <= 103.49; criterion = 1, statistic = 53.161 2) CPI <= 101.5; criterion = 1, statistic = 17.4 3)* weights = 26 2) CPI > 101.5 4)* weights = 9 1) CPI > 103.49 5) werklozen <= 613438; criterion = 0.999, statistic = 13.164 6)* weights = 18 5) werklozen > 613438 7)* weights = 31 > postscript(file="/var/www/html/rcomp/tmp/2w1oe1293006243.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/3w1oe1293006243.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 21454 24214.00 -2760.00000 2 23899 24214.00 -315.00000 3 24939 24214.00 725.00000 4 23580 24214.00 -634.00000 5 24562 24214.00 348.00000 6 24696 24214.00 482.00000 7 23785 24214.00 -429.00000 8 23812 24214.00 -402.00000 9 21917 24214.00 -2297.00000 10 19713 24214.00 -4501.00000 11 19282 24214.00 -4932.00000 12 18788 24214.00 -5426.00000 13 21453 24214.00 -2761.00000 14 24482 24214.00 268.00000 15 27474 24214.00 3260.00000 16 27264 24214.00 3050.00000 17 27349 24214.00 3135.00000 18 30632 24214.00 6418.00000 19 29429 24214.00 5215.00000 20 30084 24214.00 5870.00000 21 26290 24214.00 2076.00000 22 24379 24214.00 165.00000 23 23335 24214.00 -879.00000 24 21346 24214.00 -2868.00000 25 21106 24214.00 -3108.00000 26 24514 24214.00 300.00000 27 28353 30754.78 -2401.77778 28 30805 30754.78 50.22222 29 31348 30754.78 593.22222 30 34556 30754.78 3801.22222 31 33855 30754.78 3100.22222 32 34787 38065.55 -3278.54839 33 32529 38065.55 -5536.54839 34 29998 30754.78 -756.77778 35 29257 30754.78 -1497.77778 36 28155 30754.78 -2599.77778 37 30466 30754.78 -288.77778 38 35704 38065.55 -2361.54839 39 39327 38065.55 1261.45161 40 39351 38065.55 1285.45161 41 42234 44042.67 -1808.66667 42 43630 44042.67 -412.66667 43 43722 44042.67 -320.66667 44 43121 44042.67 -921.66667 45 37985 38065.55 -80.54839 46 37135 38065.55 -930.54839 47 34646 38065.55 -3419.54839 48 33026 38065.55 -5039.54839 49 35087 38065.55 -2978.54839 50 38846 38065.55 780.45161 51 42013 38065.55 3947.45161 52 43908 38065.55 5842.45161 53 42868 44042.67 -1174.66667 54 44423 44042.67 380.33333 55 44167 44042.67 124.33333 56 43636 44042.67 -406.66667 57 44382 38065.55 6316.45161 58 42142 44042.67 -1900.66667 59 43452 44042.67 -590.66667 60 36912 44042.67 -7130.66667 61 42413 38065.55 4347.45161 62 45344 38065.55 7278.45161 63 44873 38065.55 6807.45161 64 47510 44042.67 3467.33333 65 49554 44042.67 5511.33333 66 47369 44042.67 3326.33333 67 45998 44042.67 1955.33333 68 48140 38065.55 10074.45161 69 48441 44042.67 4398.33333 70 44928 44042.67 885.33333 71 40454 38065.55 2388.45161 72 38661 44042.67 -5381.66667 73 37246 38065.55 -819.54839 74 36843 38065.55 -1222.54839 75 36424 38065.55 -1641.54839 76 37594 38065.55 -471.54839 77 38144 38065.55 78.45161 78 38737 38065.55 671.45161 79 34560 38065.55 -3505.54839 80 36080 38065.55 -1985.54839 81 33508 38065.55 -4557.54839 82 35462 38065.55 -2603.54839 83 33374 38065.55 -4691.54839 84 32110 38065.55 -5955.54839 > 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/4osnz1293006243.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/5ab4n1293006243.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/6dtkb1293006243.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/7zc1h1293006243.tab") + } > > try(system("convert tmp/2w1oe1293006243.ps tmp/2w1oe1293006243.png",intern=TRUE)) character(0) > try(system("convert tmp/3w1oe1293006243.ps tmp/3w1oe1293006243.png",intern=TRUE)) character(0) > try(system("convert tmp/4osnz1293006243.ps tmp/4osnz1293006243.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.373 0.544 4.970