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(6.3,0,3,2.1,3.40602894496362,4,9.1,1.02325245963371,4,15.8,-1.69897000433602,1,5.2,2.20411998265592,4,10.9,0.51851393987789,1,8.3,1.71733758272386,1,11,-0.36653154442041,4,3.2,2.66745295288995,5,6.3,-1.09691001300806,1,6.6,-0.10237290870956,2,9.5,-0.69897000433602,2,3.3,1.44185217577329,5,11,-0.92081875395238,2,4.7,1.92941892571429,1,10.4,-1,3,7.4,0.01703333929878,4,2.1,2.71683772329952,5,17.9,-2,1,6.1,1.79239168949825,1,11.9,-1.69897000433602,3,13.8,0.23044892137827,1,14.3,0.54406804435028,1,15.2,-0.31875876262441,2,10,1,4,11.9,0.20951501454263,2,6.5,2.28330122870355,4,7.5,0.39794000867204,5,10.6,-0.55284196865778,3,7.4,0.62736585659273,1,8.4,0.83250891270624,2,5.7,-0.1249387366083,2,4.9,0.55630250076729,3,3.2,1.74429298312268,5,11,-0.045757490560675,2,4.9,0.30102999566398,3,13.2,-1,2,9.7,0.6222140229663,4,12.8,0.54406804435028,1),dim=c(3,39),dimnames=list(c('SWS','logWb','D'),1:39)) > y <- array(NA,dim=c(3,39),dimnames=list(c('SWS','logWb','D'),1:39)) > for (i in 1:dim(x)[1]) + { + for (j in 1:dim(x)[2]) + { + y[i,j] <- as.numeric(x[i,j]) + } + } > par2 = 'No' > par1 = '1' > par2 <- 'No' > 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(rpart) > library(partykit) Loading required package: grid Loading required package: mvtnorm > par1 <- as.numeric(par1) > autoprune <- function ( tree, method='Minimum CV'){ + xerr <- tree$cptable[,'xerror'] + cpmin.id <- which.min(xerr) + if (method == 'Minimum CV Error plus 1 SD'){ + xstd <- tree$cptable[,'xstd'] + errt <- xerr[cpmin.id] + xstd[cpmin.id] + cpSE1.min <- which.min( errt < xerr ) + mycp <- (tree$cptable[,'CP'])[cpSE1.min] + } + if (method == 'Minimum CV') { + mycp <- (tree$cptable[,'CP'])[cpmin.id] + } + return (mycp) + } > conf.multi.mat <- function(true, new) + { + if ( all( is.na(match( levels(true),levels(new) ) )) ) + stop ( 'conflict of vector levels') + multi.t <- list() + for (mylev in levels(true) ) { + true.tmp <- true + new.tmp <- new + left.lev <- levels (true.tmp)[- match(mylev,levels(true) ) ] + levels(true.tmp) <- list ( mylev = mylev, all = left.lev ) + levels(new.tmp) <- list ( mylev = mylev, all = left.lev ) + curr.t <- conf.mat ( true.tmp , new.tmp ) + multi.t[[mylev]] <- curr.t + multi.t[[mylev]]$precision <- + round( curr.t$conf[1,1] / sum( curr.t$conf[1,] ), 2 ) + } + return (multi.t) + } > x <- t(y) > k <- length(x[1,]) > n <- length(x[,1]) > x1 <- cbind(x[,par1], x[,1:k!=par1]) > mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1]) > colnames(x1) <- mycolnames #colnames(x)[par1] > m <- rpart(as.data.frame(x1)) > par2 [1] "No" > if (par2 != 'No') { + mincp <- autoprune(m,method=par2) + print(mincp) + m <- prune(m,cp=mincp) + } > m$cptable CP nsplit rel error xerror xstd 1 0.4060272 0 1.0000000 1.0802417 0.1987379 2 0.1054015 1 0.5939728 0.7730237 0.1465899 3 0.0100000 2 0.4885713 0.7272346 0.1270163 > postscript(file="/var/www/html/rcomp/tmp/1658r1272715458.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(as.party(m),tp_args=list(id=FALSE)) > dev.off() null device 1 > postscript(file="/var/www/html/rcomp/tmp/2658r1272715458.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plotcp(m) > dev.off() null device 1 > cbind(y=m$y,pred=predict(m),res=residuals(m)) y pred res 1 6.3 9.033333 -2.73333333 2 2.1 4.470000 -2.37000000 3 9.1 9.033333 0.06666667 4 15.8 12.072727 3.72727273 5 5.2 4.470000 0.73000000 6 10.9 9.033333 1.86666667 7 8.3 4.470000 3.83000000 8 11.0 12.072727 -1.07272727 9 3.2 4.470000 -1.27000000 10 6.3 12.072727 -5.77272727 11 6.6 9.033333 -2.43333333 12 9.5 12.072727 -2.57272727 13 3.3 4.470000 -1.17000000 14 11.0 12.072727 -1.07272727 15 4.7 4.470000 0.23000000 16 10.4 12.072727 -1.67272727 17 7.4 9.033333 -1.63333333 18 2.1 4.470000 -2.37000000 19 17.9 12.072727 5.82727273 20 6.1 4.470000 1.63000000 21 11.9 12.072727 -0.17272727 22 13.8 9.033333 4.76666667 23 14.3 9.033333 5.26666667 24 15.2 12.072727 3.12727273 25 10.0 9.033333 0.96666667 26 11.9 9.033333 2.86666667 27 6.5 4.470000 2.03000000 28 7.5 9.033333 -1.53333333 29 10.6 12.072727 -1.47272727 30 7.4 9.033333 -1.63333333 31 8.4 9.033333 -0.63333333 32 5.7 9.033333 -3.33333333 33 4.9 9.033333 -4.13333333 34 3.2 4.470000 -1.27000000 35 11.0 9.033333 1.96666667 36 4.9 9.033333 -4.13333333 37 13.2 12.072727 1.12727273 38 9.7 9.033333 0.66666667 39 12.8 9.033333 3.76666667 > myr <- residuals(m) > myp <- predict(m) > postscript(file="/var/www/html/rcomp/tmp/3658r1272715458.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > op <- par(mfrow=c(2,2)) > plot(myr,ylab='residuals') > plot(density(myr),main='Residual Kernel Density') > plot(myp,myr,xlab='predicted',ylab='residuals',main='Predicted vs Residuals') > plot(density(myp),main='Prediction Kernel Density') > par(op) > dev.off() null device 1 > > #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/rcomp/createtable") > > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Model Performance',6,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'#',header=TRUE) > a<-table.element(a,'Complexity',header=TRUE) > a<-table.element(a,'split',header=TRUE) > a<-table.element(a,'relative error',header=TRUE) > a<-table.element(a,'CV error',header=TRUE) > a<-table.element(a,'CV S.D.',header=TRUE) > a<-table.row.end(a) > for (i in 1:length(m$cptable[,1])) { + a<-table.row.start(a) + a<-table.element(a,i,header=TRUE) + a<-table.element(a,round(m$cptable[i,'CP'],3)) + a<-table.element(a,m$cptable[i,'nsplit']) + a<-table.element(a,round(m$cptable[i,'rel error'],3)) + a<-table.element(a,round(m$cptable[i,'xerror'],3)) + a<-table.element(a,round(m$cptable[i,'xstd'],3)) + a<-table.row.end(a) + } > a<-table.end(a) > table.save(a,file="/var/www/html/rcomp/tmp/42xni1272715458.tab") > > try(system("convert tmp/1658r1272715458.ps tmp/1658r1272715458.png",intern=TRUE)) character(0) > try(system("convert tmp/2658r1272715458.ps tmp/2658r1272715458.png",intern=TRUE)) character(0) > try(system("convert tmp/3658r1272715458.ps tmp/3658r1272715458.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 0.958 0.451 1.154