R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-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(589
+ ,248.85
+ ,65453
+ ,559
+ ,249.68
+ ,65715
+ ,623
+ ,251.13
+ ,66261
+ ,617
+ ,251.24
+ ,66332
+ ,603
+ ,253.24
+ ,66229
+ ,558
+ ,254.66
+ ,66579
+ ,609
+ ,255.85
+ ,66817
+ ,583
+ ,256.93
+ ,67373
+ ,570
+ ,258.99
+ ,68078
+ ,543
+ ,258.30
+ ,69137
+ ,598
+ ,260.53
+ ,69816
+ ,569
+ ,260.65
+ ,70252
+ ,552
+ ,260.98
+ ,70389
+ ,514
+ ,262.09
+ ,70572
+ ,569
+ ,263.18
+ ,70780
+ ,529
+ ,262.62
+ ,70912
+ ,515
+ ,263.18
+ ,71594
+ ,481
+ ,264.91
+ ,72587
+ ,536
+ ,265.20
+ ,73677
+ ,498
+ ,266.14
+ ,74712
+ ,446
+ ,268.15
+ ,75208
+ ,503
+ ,270.62
+ ,75657
+ ,470
+ ,272.65
+ ,76011
+ ,458
+ ,274.50
+ ,76748
+ ,437
+ ,274.37
+ ,76537
+ ,502
+ ,277.85
+ ,76622
+ ,482
+ ,280.15
+ ,76404
+ ,474
+ ,280.67
+ ,76219
+ ,457
+ ,281.42
+ ,76875
+ ,522
+ ,283.23
+ ,77374
+ ,513
+ ,283.34
+ ,77743
+ ,515
+ ,284.09
+ ,78030
+ ,506
+ ,285.47
+ ,77805
+ ,576
+ ,287.27
+ ,77905
+ ,556
+ ,287.96
+ ,78158
+ ,559
+ ,289.05
+ ,78616
+ ,541
+ ,289.84
+ ,79740
+ ,606
+ ,292.68
+ ,80312
+ ,600
+ ,294.61
+ ,80921
+ ,588
+ ,296.22
+ ,81078
+ ,570
+ ,296.70
+ ,81394
+ ,626
+ ,300.82
+ ,81787
+ ,601
+ ,303.57
+ ,82252
+ ,588
+ ,304.32
+ ,82854
+ ,573
+ ,304.52
+ ,83498
+ ,622
+ ,306.69
+ ,83811
+ ,570
+ ,308.73
+ ,84531
+ ,547
+ ,308.30
+ ,85330
+ ,512
+ ,309.67
+ ,86247
+ ,554
+ ,311.68
+ ,86386
+ ,517
+ ,312.62
+ ,86918
+ ,506
+ ,315.18
+ ,87184
+ ,479
+ ,320.19
+ ,87843
+ ,527
+ ,325.96
+ ,88204
+ ,508
+ ,330.45
+ ,87675
+ ,532
+ ,329.16
+ ,85964
+ ,532
+ ,327.53
+ ,84387
+ ,588
+ ,326.87
+ ,84530
+ ,566
+ ,326.52
+ ,85497
+ ,573
+ ,326.65
+ ,85968
+ ,545
+ ,329.25
+ ,86030
+ ,597
+ ,333.11
+ ,86963
+ ,555
+ ,334.51
+ ,87324
+ ,548
+ ,336.21
+ ,87770
+ ,524
+ ,339.91
+ ,88534
+ ,572
+ ,344.53
+ ,88888)
+ ,dim=c(3
+ ,66)
+ ,dimnames=list(c('Werkloosheid'
+ ,'CPI'
+ ,'BBP')
+ ,1:66))
>  y <- array(NA,dim=c(3,66),dimnames=list(c('Werkloosheid','CPI','BBP'),1:66))
>  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 = '3'
> #'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
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] "BBP"
> x[,par1]
 [1] 65453 65715 66261 66332 66229 66579 66817 67373 68078 69137 69816 70252
[13] 70389 70572 70780 70912 71594 72587 73677 74712 75208 75657 76011 76748
[25] 76537 76622 76404 76219 76875 77374 77743 78030 77805 77905 78158 78616
[37] 79740 80312 80921 81078 81394 81787 82252 82854 83498 83811 84531 85330
[49] 86247 86386 86918 87184 87843 88204 87675 85964 84387 84530 85497 85968
[61] 86030 86963 87324 87770 88534 88888
> 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])

65453 65715 66229 66261 66332 66579 66817 67373 68078 69137 69816 70252 70389 
    1     1     1     1     1     1     1     1     1     1     1     1     1 
70572 70780 70912 71594 72587 73677 74712 75208 75657 76011 76219 76404 76537 
    1     1     1     1     1     1     1     1     1     1     1     1     1 
76622 76748 76875 77374 77743 77805 77905 78030 78158 78616 79740 80312 80921 
    1     1     1     1     1     1     1     1     1     1     1     1     1 
81078 81394 81787 82252 82854 83498 83811 84387 84530 84531 85330 85497 85964 
    1     1     1     1     1     1     1     1     1     1     1     1     1 
85968 86030 86247 86386 86918 86963 87184 87324 87675 87770 87843 88204 88534 
    1     1     1     1     1     1     1     1     1     1     1     1     1 
88888 
    1 
> colnames(x)
[1] "Werkloosheid" "CPI"          "BBP"         
> colnames(x)[par1]
[1] "BBP"
> x[,par1]
 [1] 65453 65715 66261 66332 66229 66579 66817 67373 68078 69137 69816 70252
[13] 70389 70572 70780 70912 71594 72587 73677 74712 75208 75657 76011 76748
[25] 76537 76622 76404 76219 76875 77374 77743 78030 77805 77905 78158 78616
[37] 79740 80312 80921 81078 81394 81787 82252 82854 83498 83811 84531 85330
[49] 86247 86386 86918 87184 87843 88204 87675 85964 84387 84530 85497 85968
[61] 86030 86963 87324 87770 88534 88888
> if (par2 == 'none') {
+ m <- ctree(as.formula(paste(colnames(x)[par1],' ~ .',sep='')),data = x)
+ }
> 
> #Note: the /var/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/1lsje1324134062.tab") 
+ }
+ }
> m

	 Conditional inference tree with 4 terminal nodes

Response:  BBP 
Inputs:  Werkloosheid, CPI 
Number of observations:  66 

1) CPI <= 289.05; criterion = 1, statistic = 60.605
  2) CPI <= 263.18; criterion = 1, statistic = 32.375
    3)*  weights = 17 
  2) CPI > 263.18
    4)*  weights = 19 
1) CPI > 289.05
  5) CPI <= 306.69; criterion = 1, statistic = 21.674
    6)*  weights = 10 
  5) CPI > 306.69
    7)*  weights = 20 
> postscript(file="/var/www/rcomp/tmp/2uifs1324134062.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/rcomp/tmp/38hzw1324134062.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    65453  68369.94 -2916.94118
2    65715  68369.94 -2654.94118
3    66261  68369.94 -2108.94118
4    66332  68369.94 -2037.94118
5    66229  68369.94 -2140.94118
6    66579  68369.94 -1790.94118
7    66817  68369.94 -1552.94118
8    67373  68369.94  -996.94118
9    68078  68369.94  -291.94118
10   69137  68369.94   767.05882
11   69816  68369.94  1446.05882
12   70252  68369.94  1882.05882
13   70389  68369.94  2019.05882
14   70572  68369.94  2202.05882
15   70780  68369.94  2410.05882
16   70912  68369.94  2542.05882
17   71594  68369.94  3224.05882
18   72587  76467.79 -3880.78947
19   73677  76467.79 -2790.78947
20   74712  76467.79 -1755.78947
21   75208  76467.79 -1259.78947
22   75657  76467.79  -810.78947
23   76011  76467.79  -456.78947
24   76748  76467.79   280.21053
25   76537  76467.79    69.21053
26   76622  76467.79   154.21053
27   76404  76467.79   -63.78947
28   76219  76467.79  -248.78947
29   76875  76467.79   407.21053
30   77374  76467.79   906.21053
31   77743  76467.79  1275.21053
32   78030  76467.79  1562.21053
33   77805  76467.79  1337.21053
34   77905  76467.79  1437.21053
35   78158  76467.79  1690.21053
36   78616  76467.79  2148.21053
37   79740  81764.70 -2024.70000
38   80312  81764.70 -1452.70000
39   80921  81764.70  -843.70000
40   81078  81764.70  -686.70000
41   81394  81764.70  -370.70000
42   81787  81764.70    22.30000
43   82252  81764.70   487.30000
44   82854  81764.70  1089.30000
45   83498  81764.70  1733.30000
46   83811  81764.70  2046.30000
47   84531  86608.65 -2077.65000
48   85330  86608.65 -1278.65000
49   86247  86608.65  -361.65000
50   86386  86608.65  -222.65000
51   86918  86608.65   309.35000
52   87184  86608.65   575.35000
53   87843  86608.65  1234.35000
54   88204  86608.65  1595.35000
55   87675  86608.65  1066.35000
56   85964  86608.65  -644.65000
57   84387  86608.65 -2221.65000
58   84530  86608.65 -2078.65000
59   85497  86608.65 -1111.65000
60   85968  86608.65  -640.65000
61   86030  86608.65  -578.65000
62   86963  86608.65   354.35000
63   87324  86608.65   715.35000
64   87770  86608.65  1161.35000
65   88534  86608.65  1925.35000
66   88888  86608.65  2279.35000
> 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/rcomp/tmp/4u3fl1324134062.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/rcomp/tmp/5o37s1324134063.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/rcomp/tmp/6snrp1324134063.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/rcomp/tmp/7z96u1324134063.tab") 
+ }
> 
> try(system("convert tmp/2uifs1324134062.ps tmp/2uifs1324134062.png",intern=TRUE))
character(0)
> try(system("convert tmp/38hzw1324134062.ps tmp/38hzw1324134062.png",intern=TRUE))
character(0)
> try(system("convert tmp/4u3fl1324134062.ps tmp/4u3fl1324134062.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  2.350   0.120   2.478