R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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.

  Natural language support but running in an English locale

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(1606
+ ,6
+ ,3.74
+ ,16
+ ,1391
+ ,1634
+ ,6.81
+ ,4.17
+ ,29
+ ,1621
+ ,2013
+ ,9.75
+ ,4.84
+ ,22
+ ,1837
+ ,1654
+ ,6.96
+ ,4.21
+ ,30
+ ,2132
+ ,1003
+ ,3.94
+ ,3.93
+ ,20
+ ,1489
+ ,1029
+ ,5
+ ,4.9
+ ,39
+ ,1817
+ ,1052
+ ,4.9
+ ,4.7
+ ,18
+ ,1586
+ ,1653
+ ,5.7
+ ,3.5
+ ,9.6
+ ,1565
+ ,1918
+ ,6.5
+ ,3.4
+ ,10.2
+ ,1787
+ ,1926
+ ,7.1
+ ,3.7
+ ,20.2
+ ,1804
+ ,1862
+ ,7.5
+ ,4
+ ,50
+ ,1763
+ ,1816
+ ,7.8
+ ,4.3
+ ,120
+ ,1675
+ ,1712
+ ,7
+ ,4.1
+ ,19.8
+ ,1575
+ ,1646
+ ,7.4
+ ,4.5
+ ,18
+ ,1524
+ ,1555
+ ,8.55
+ ,5.5
+ ,3
+ ,1686
+ ,1402
+ ,7.43
+ ,5.3
+ ,11
+ ,1800
+ ,1047
+ ,4.7
+ ,4.5
+ ,15
+ ,1442
+ ,891
+ ,4.7
+ ,5.3
+ ,27
+ ,1345
+ ,940
+ ,5.3
+ ,5.6
+ ,28
+ ,1500
+ ,1372
+ ,6.2
+ ,4.5
+ ,14
+ ,1556
+ ,2012
+ ,7.4
+ ,3.7
+ ,5.6
+ ,2012
+ ,1879
+ ,7.5
+ ,4
+ ,6.5
+ ,1618
+ ,1667
+ ,7.32
+ ,4.4
+ ,8.5
+ ,1487
+ ,1856
+ ,8.15
+ ,4.4
+ ,87.9
+ ,1607
+ ,1771
+ ,7.24
+ ,4.1
+ ,5.8
+ ,1308
+ ,1721
+ ,7.4
+ ,4.3
+ ,25.2
+ ,1429
+ ,1773
+ ,9.4
+ ,5.3
+ ,7.5
+ ,1596
+ ,1507
+ ,8.9
+ ,5.9
+ ,13.7
+ ,1884
+ ,1033
+ ,4.5
+ ,4.4
+ ,34
+ ,1262
+ ,1011
+ ,4.9
+ ,4.9
+ ,17
+ ,1283
+ ,1111
+ ,5.6
+ ,5.1
+ ,9
+ ,1346
+ ,1736
+ ,6.4
+ ,3.7
+ ,9.2
+ ,1505
+ ,1865
+ ,6
+ ,3.2
+ ,5
+ ,1151
+ ,2078
+ ,6.9
+ ,3.3
+ ,24
+ ,1600
+ ,1947
+ ,6.7
+ ,3.5
+ ,40
+ ,1420
+ ,1428
+ ,5.4
+ ,3.8
+ ,86.5
+ ,1073
+ ,1500
+ ,5.6
+ ,3.8
+ ,0.54
+ ,1076
+ ,1950
+ ,6.9
+ ,3.5
+ ,14
+ ,1510
+ ,1591
+ ,6.9
+ ,4.3
+ ,4.8
+ ,1345
+ ,1613
+ ,7
+ ,4.3
+ ,28
+ ,1631
+ ,1077
+ ,4
+ ,3.7
+ ,16
+ ,1135
+ ,880
+ ,3.7
+ ,4.2
+ ,5.8
+ ,1009
+ ,1128
+ ,4.9
+ ,4.3
+ ,16
+ ,1155
+ ,1320
+ ,5
+ ,3.8
+ ,9.1
+ ,1184
+ ,1692
+ ,5.7
+ ,3.4
+ ,6
+ ,1285
+ ,1575
+ ,6.1
+ ,3.9
+ ,17
+ ,1257
+ ,1478
+ ,5.3
+ ,3.6
+ ,26
+ ,1131
+ ,1500
+ ,5.5
+ ,3.6
+ ,99.6
+ ,1274
+ ,1368
+ ,5.7
+ ,4.2
+ ,41
+ ,235
+ ,1563
+ ,5.21
+ ,3.3
+ ,72
+ ,1299
+ ,1424
+ ,5.4
+ ,3.8
+ ,23
+ ,1460
+ ,1274
+ ,4.5
+ ,3.5
+ ,42
+ ,1455
+ ,1047
+ ,3.7
+ ,3.7
+ ,40
+ ,1113
+ ,1049
+ ,4.1
+ ,3.9
+ ,18
+ ,1263
+ ,1069
+ ,4.8
+ ,4.5
+ ,45
+ ,1401
+ ,981
+ ,4.1
+ ,4.2
+ ,18
+ ,1135
+ ,1540
+ ,5
+ ,3.2
+ ,2
+ ,1137
+ ,1559
+ ,5.2
+ ,3.3
+ ,10
+ ,1140
+ ,1459
+ ,5.5
+ ,3.8
+ ,13.6
+ ,1014
+ ,1559
+ ,5.9
+ ,3.8
+ ,160
+ ,1220)
+ ,dim=c(5
+ ,60)
+ ,dimnames=list(c('aanvoer'
+ ,'aanvoerwaarde'
+ ,'prijzen'
+ ,'interventie'
+ ,'visserijen')
+ ,1:60))
>  y <- array(NA,dim=c(5,60),dimnames=list(c('aanvoer','aanvoerwaarde','prijzen','interventie','visserijen'),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 = ''
> par2 = 'none'
> par1 = '4'
> #'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] "interventie"
> x[,par1]
 [1]  16.00  29.00  22.00  30.00  20.00  39.00  18.00   9.60  10.20  20.20
[11]  50.00 120.00  19.80  18.00   3.00  11.00  15.00  27.00  28.00  14.00
[21]   5.60   6.50   8.50  87.90   5.80  25.20   7.50  13.70  34.00  17.00
[31]   9.00   9.20   5.00  24.00  40.00  86.50   0.54  14.00   4.80  28.00
[41]  16.00   5.80  16.00   9.10   6.00  17.00  26.00  99.60  41.00  72.00
[51]  23.00  42.00  40.00  18.00  45.00  18.00   2.00  10.00  13.60 160.00
> 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])

0.54    2    3  4.8    5  5.6  5.8    6  6.5  7.5  8.5    9  9.1  9.2  9.6   10 
   1    1    1    1    1    1    2    1    1    1    1    1    1    1    1    1 
10.2   11 13.6 13.7   14   15   16   17   18 19.8   20 20.2   22   23   24 25.2 
   1    1    1    1    2    1    3    2    4    1    1    1    1    1    1    1 
  26   27   28   29   30   34   39   40   41   42   45   50   72 86.5 87.9 99.6 
   1    1    2    1    1    1    1    2    1    1    1    1    1    1    1    1 
 120  160 
   1    1 
> colnames(x)
[1] "aanvoer"       "aanvoerwaarde" "prijzen"       "interventie"  
[5] "visserijen"   
> colnames(x)[par1]
[1] "interventie"
> x[,par1]
 [1]  16.00  29.00  22.00  30.00  20.00  39.00  18.00   9.60  10.20  20.20
[11]  50.00 120.00  19.80  18.00   3.00  11.00  15.00  27.00  28.00  14.00
[21]   5.60   6.50   8.50  87.90   5.80  25.20   7.50  13.70  34.00  17.00
[31]   9.00   9.20   5.00  24.00  40.00  86.50   0.54  14.00   4.80  28.00
[41]  16.00   5.80  16.00   9.10   6.00  17.00  26.00  99.60  41.00  72.00
[51]  23.00  42.00  40.00  18.00  45.00  18.00   2.00  10.00  13.60 160.00
> if (par2 == 'none') {
+ m <- ctree(as.formula(paste(colnames(x)[par1],' ~ .',sep='')),data = x)
+ }
> 
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/1ei6v1292348126.tab") 
+ }
+ }
> m

	 Conditional inference tree with 1 terminal nodes

Response:  interventie 
Inputs:  aanvoer, aanvoerwaarde, prijzen, visserijen 
Number of observations:  60 

1)*  weights = 60 
> postscript(file="/var/www/html/freestat/rcomp/tmp/2ei6v1292348126.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/freestat/rcomp/tmp/3ei6v1292348126.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    16.00  27.21067 -11.2106667
2    29.00  27.21067   1.7893333
3    22.00  27.21067  -5.2106667
4    30.00  27.21067   2.7893333
5    20.00  27.21067  -7.2106667
6    39.00  27.21067  11.7893333
7    18.00  27.21067  -9.2106667
8     9.60  27.21067 -17.6106667
9    10.20  27.21067 -17.0106667
10   20.20  27.21067  -7.0106667
11   50.00  27.21067  22.7893333
12  120.00  27.21067  92.7893333
13   19.80  27.21067  -7.4106667
14   18.00  27.21067  -9.2106667
15    3.00  27.21067 -24.2106667
16   11.00  27.21067 -16.2106667
17   15.00  27.21067 -12.2106667
18   27.00  27.21067  -0.2106667
19   28.00  27.21067   0.7893333
20   14.00  27.21067 -13.2106667
21    5.60  27.21067 -21.6106667
22    6.50  27.21067 -20.7106667
23    8.50  27.21067 -18.7106667
24   87.90  27.21067  60.6893333
25    5.80  27.21067 -21.4106667
26   25.20  27.21067  -2.0106667
27    7.50  27.21067 -19.7106667
28   13.70  27.21067 -13.5106667
29   34.00  27.21067   6.7893333
30   17.00  27.21067 -10.2106667
31    9.00  27.21067 -18.2106667
32    9.20  27.21067 -18.0106667
33    5.00  27.21067 -22.2106667
34   24.00  27.21067  -3.2106667
35   40.00  27.21067  12.7893333
36   86.50  27.21067  59.2893333
37    0.54  27.21067 -26.6706667
38   14.00  27.21067 -13.2106667
39    4.80  27.21067 -22.4106667
40   28.00  27.21067   0.7893333
41   16.00  27.21067 -11.2106667
42    5.80  27.21067 -21.4106667
43   16.00  27.21067 -11.2106667
44    9.10  27.21067 -18.1106667
45    6.00  27.21067 -21.2106667
46   17.00  27.21067 -10.2106667
47   26.00  27.21067  -1.2106667
48   99.60  27.21067  72.3893333
49   41.00  27.21067  13.7893333
50   72.00  27.21067  44.7893333
51   23.00  27.21067  -4.2106667
52   42.00  27.21067  14.7893333
53   40.00  27.21067  12.7893333
54   18.00  27.21067  -9.2106667
55   45.00  27.21067  17.7893333
56   18.00  27.21067  -9.2106667
57    2.00  27.21067 -25.2106667
58   10.00  27.21067 -17.2106667
59   13.60  27.21067 -13.6106667
60  160.00  27.21067 132.7893333
> 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/freestat/rcomp/tmp/4h14j1292348126.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/freestat/rcomp/tmp/52j3p1292348126.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/freestat/rcomp/tmp/6daka1292348126.tab") 
+ }
Warning message:
In cor(result$Forecasts, result$Actuals) : the standard deviation is zero
> 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/freestat/rcomp/tmp/7hbig1292348126.tab") 
+ }
> try(system("convert tmp/2ei6v1292348126.ps tmp/2ei6v1292348126.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ei6v1292348126.ps tmp/3ei6v1292348126.png",intern=TRUE))
character(0)
> try(system("convert tmp/4h14j1292348126.ps tmp/4h14j1292348126.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  3.464   0.707   3.655