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(107.11
+ ,236.67
+ ,8.92
+ ,1
+ ,122.23
+ ,258.1
+ ,9.32
+ ,2
+ ,134.69
+ ,241.52
+ ,8.9
+ ,3
+ ,128.79
+ ,190.71
+ ,8.53
+ ,4
+ ,126.16
+ ,200.32
+ ,8.51
+ ,5
+ ,119.98
+ ,223.41
+ ,9.03
+ ,6
+ ,108.45
+ ,201.38
+ ,9.6
+ ,7
+ ,108.43
+ ,211.83
+ ,9.88
+ ,8
+ ,98.17
+ ,224.41
+ ,10.81
+ ,9
+ ,106.09
+ ,211.57
+ ,11.61
+ ,10
+ ,108.81
+ ,194.77
+ ,11.81
+ ,11
+ ,103.03
+ ,201.86
+ ,13.93
+ ,12
+ ,124.36
+ ,225
+ ,16.19
+ ,13
+ ,118.52
+ ,278.9
+ ,18.05
+ ,14
+ ,112.2
+ ,259.74
+ ,17.08
+ ,15
+ ,114.71
+ ,230.45
+ ,17.46
+ ,16
+ ,107.96
+ ,238.26
+ ,16.9
+ ,17
+ ,101.21
+ ,250.14
+ ,15.69
+ ,18
+ ,102.77
+ ,263.81
+ ,15.86
+ ,19
+ ,112.13
+ ,247.22
+ ,12.98
+ ,20
+ ,109.36
+ ,229.81
+ ,12.31
+ ,21
+ ,110.91
+ ,224.27
+ ,11.51
+ ,22
+ ,123.57
+ ,213.23
+ ,11.73
+ ,23
+ ,129.95
+ ,239.57
+ ,11.7
+ ,24
+ ,124.46
+ ,249.7
+ ,10.9
+ ,25
+ ,122.34
+ ,212.5
+ ,10.57
+ ,26
+ ,116.61
+ ,203.27
+ ,10.37
+ ,27
+ ,114.59
+ ,192.05
+ ,9.59
+ ,28
+ ,112.52
+ ,190.04
+ ,9.09
+ ,29
+ ,118.67
+ ,202.05
+ ,9.26
+ ,30
+ ,116.8
+ ,211.91
+ ,9.9
+ ,31
+ ,123.63
+ ,210.39
+ ,9.61
+ ,32
+ ,128.04
+ ,231.25
+ ,9.85
+ ,33
+ ,134.57
+ ,224.3
+ ,9.99
+ ,34
+ ,130.33
+ ,209.64
+ ,9.9
+ ,35
+ ,136.47
+ ,206.05
+ ,10.45
+ ,36
+ ,139.05
+ ,229.7
+ ,11.66
+ ,37
+ ,158.21
+ ,264.67
+ ,13.61
+ ,38
+ ,148.07
+ ,246.29
+ ,12.88
+ ,39
+ ,137.74
+ ,260.91
+ ,12.52
+ ,40
+ ,139.74
+ ,265.14
+ ,10.93
+ ,41
+ ,144.08
+ ,284.52
+ ,12.07
+ ,42
+ ,145.35
+ ,287.48
+ ,13.21
+ ,43
+ ,145.77
+ ,321.9
+ ,13.68
+ ,44
+ ,140.56
+ ,321.59
+ ,14.02
+ ,45
+ ,121.41
+ ,282.39
+ ,11.7
+ ,46
+ ,120.44
+ ,241
+ ,11.83
+ ,47
+ ,116.97
+ ,228.48
+ ,11.32
+ ,48
+ ,128.03
+ ,261.59
+ ,12.24
+ ,49
+ ,128.51
+ ,270
+ ,13.31
+ ,50
+ ,127.76
+ ,262.86
+ ,12.93
+ ,51
+ ,134.58
+ ,277.41
+ ,13.47
+ ,52
+ ,147.64
+ ,288
+ ,15.47
+ ,53
+ ,144.46
+ ,287.14
+ ,16.58
+ ,54
+ ,137.6
+ ,337.65
+ ,17.8
+ ,55
+ ,146.87
+ ,328.38
+ ,21.72
+ ,56
+ ,145.67
+ ,374.41
+ ,23.45
+ ,57
+ ,151.95
+ ,344.77
+ ,23.16
+ ,58
+ ,150.23
+ ,361.05
+ ,22.77
+ ,59
+ ,155.86
+ ,374.22
+ ,24.9
+ ,60)
+ ,dim=c(4
+ ,60)
+ ,dimnames=list(c('Coffee'
+ ,'Tea'
+ ,'Sugar'
+ ,'Month')
+ ,1:60))
>  y <- array(NA,dim=c(4,60),dimnames=list(c('Coffee','Tea','Sugar','Month'),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 = '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] "Coffee"
> x[,par1]
 [1] 107.11 122.23 134.69 128.79 126.16 119.98 108.45 108.43  98.17 106.09
[11] 108.81 103.03 124.36 118.52 112.20 114.71 107.96 101.21 102.77 112.13
[21] 109.36 110.91 123.57 129.95 124.46 122.34 116.61 114.59 112.52 118.67
[31] 116.80 123.63 128.04 134.57 130.33 136.47 139.05 158.21 148.07 137.74
[41] 139.74 144.08 145.35 145.77 140.56 121.41 120.44 116.97 128.03 128.51
[51] 127.76 134.58 147.64 144.46 137.60 146.87 145.67 151.95 150.23 155.86
> 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])

 98.17 101.21 102.77 103.03 106.09 107.11 107.96 108.43 108.45 108.81 109.36 
     1      1      1      1      1      1      1      1      1      1      1 
110.91 112.13  112.2 112.52 114.59 114.71 116.61  116.8 116.97 118.52 118.67 
     1      1      1      1      1      1      1      1      1      1      1 
119.98 120.44 121.41 122.23 122.34 123.57 123.63 124.36 124.46 126.16 127.76 
     1      1      1      1      1      1      1      1      1      1      1 
128.03 128.04 128.51 128.79 129.95 130.33 134.57 134.58 134.69 136.47  137.6 
     1      1      1      1      1      1      1      1      1      1      1 
137.74 139.05 139.74 140.56 144.08 144.46 145.35 145.67 145.77 146.87 147.64 
     1      1      1      1      1      1      1      1      1      1      1 
148.07 150.23 151.95 155.86 158.21 
     1      1      1      1      1 
> colnames(x)
[1] "Coffee" "Tea"    "Sugar"  "Month" 
> colnames(x)[par1]
[1] "Coffee"
> x[,par1]
 [1] 107.11 122.23 134.69 128.79 126.16 119.98 108.45 108.43  98.17 106.09
[11] 108.81 103.03 124.36 118.52 112.20 114.71 107.96 101.21 102.77 112.13
[21] 109.36 110.91 123.57 129.95 124.46 122.34 116.61 114.59 112.52 118.67
[31] 116.80 123.63 128.04 134.57 130.33 136.47 139.05 158.21 148.07 137.74
[41] 139.74 144.08 145.35 145.77 140.56 121.41 120.44 116.97 128.03 128.51
[51] 127.76 134.58 147.64 144.46 137.60 146.87 145.67 151.95 150.23 155.86
> 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/1znog1292276640.tab") 
+ }
+ }
> m

	 Conditional inference tree with 3 terminal nodes

Response:  Coffee 
Inputs:  Tea, Sugar, Month 
Number of observations:  60 

1) Month <= 32; criterion = 1, statistic = 30.038
  2)*  weights = 32 
1) Month > 32
  3) Sugar <= 13.47; criterion = 0.996, statistic = 10.204
    4)*  weights = 17 
  3) Sugar > 13.47
    5)*  weights = 11 
> postscript(file="/var/www/html/freestat/rcomp/tmp/2znog1292276640.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/3znog1292276640.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   107.11  115.2878  -8.17781250
2   122.23  115.2878   6.94218750
3   134.69  115.2878  19.40218750
4   128.79  115.2878  13.50218750
5   126.16  115.2878  10.87218750
6   119.98  115.2878   4.69218750
7   108.45  115.2878  -6.83781250
8   108.43  115.2878  -6.85781250
9    98.17  115.2878 -17.11781250
10  106.09  115.2878  -9.19781250
11  108.81  115.2878  -6.47781250
12  103.03  115.2878 -12.25781250
13  124.36  115.2878   9.07218750
14  118.52  115.2878   3.23218750
15  112.20  115.2878  -3.08781250
16  114.71  115.2878  -0.57781250
17  107.96  115.2878  -7.32781250
18  101.21  115.2878 -14.07781250
19  102.77  115.2878 -12.51781250
20  112.13  115.2878  -3.15781250
21  109.36  115.2878  -5.92781250
22  110.91  115.2878  -4.37781250
23  123.57  115.2878   8.28218750
24  129.95  115.2878  14.66218750
25  124.46  115.2878   9.17218750
26  122.34  115.2878   7.05218750
27  116.61  115.2878   1.32218750
28  114.59  115.2878  -0.69781250
29  112.52  115.2878  -2.76781250
30  118.67  115.2878   3.38218750
31  116.80  115.2878   1.51218750
32  123.63  115.2878   8.34218750
33  128.04  133.0082  -4.96823529
34  134.57  133.0082   1.56176471
35  130.33  133.0082  -2.67823529
36  136.47  133.0082   3.46176471
37  139.05  133.0082   6.04176471
38  158.21  147.7109  10.49909091
39  148.07  133.0082  15.06176471
40  137.74  133.0082   4.73176471
41  139.74  133.0082   6.73176471
42  144.08  133.0082  11.07176471
43  145.35  133.0082  12.34176471
44  145.77  147.7109  -1.94090909
45  140.56  147.7109  -7.15090909
46  121.41  133.0082 -11.59823529
47  120.44  133.0082 -12.56823529
48  116.97  133.0082 -16.03823529
49  128.03  133.0082  -4.97823529
50  128.51  133.0082  -4.49823529
51  127.76  133.0082  -5.24823529
52  134.58  133.0082   1.57176471
53  147.64  147.7109  -0.07090909
54  144.46  147.7109  -3.25090909
55  137.60  147.7109 -10.11090909
56  146.87  147.7109  -0.84090909
57  145.67  147.7109  -2.04090909
58  151.95  147.7109   4.23909091
59  150.23  147.7109   2.51909091
60  155.86  147.7109   8.14909091
> 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/4ren01292276640.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/55ol91292276640.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/6962x1292276640.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/freestat/rcomp/tmp/7ky101292276640.tab") 
+ }
> 
> try(system("convert tmp/2znog1292276640.ps tmp/2znog1292276640.png",intern=TRUE))
character(0)
> try(system("convert tmp/3znog1292276640.ps tmp/3znog1292276640.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ren01292276640.ps tmp/4ren01292276640.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  3.605   0.733   3.752