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.65321251377534 + ,0 + ,0.81954393554187 + ,1.6232492903979 + ,3 + ,1 + ,3 + ,2.1 + ,1.83884909073726 + ,3.40602894496362 + ,3.66304097489397 + ,2.79518458968242 + ,3 + ,5 + ,4 + ,9.1 + ,1.43136376415899 + ,1.02325245963371 + ,2.25406445291434 + ,2.25527250510331 + ,4 + ,4 + ,4 + ,15.8 + ,1.27875360095283 + ,-1.69897000433602 + ,-0.52287874528034 + ,1.54406804435028 + ,1 + ,1 + ,1 + ,5.2 + ,1.48287358360875 + ,2.20411998265592 + ,2.22788670461367 + ,2.59328606702046 + ,4 + ,5 + ,4 + ,10.9 + ,1.44715803134222 + ,0.51851393987789 + ,1.40823996531185 + ,1.79934054945358 + ,1 + ,2 + ,1 + ,8.3 + ,1.69897000433602 + ,1.71733758272386 + ,2.64345267648619 + ,2.36172783601759 + ,1 + ,1 + ,1 + ,11 + ,0.84509804001426 + ,-0.36653154442041 + ,0.80617997398389 + ,2.04921802267018 + ,5 + ,4 + ,4 + ,3.2 + ,1.47712125471966 + ,2.66745295288995 + ,2.62634036737504 + ,2.44870631990508 + ,5 + ,5 + ,5 + ,6.3 + ,0.54406804435028 + ,-1.09691001300806 + ,0.079181246047625 + ,1.6232492903979 + ,1 + ,1 + ,1 + ,6.6 + ,0.77815125038364 + ,-0.10237290870956 + ,0.54406804435028 + ,1.6232492903979 + ,2 + ,2 + ,2 + ,9.5 + ,1.01703333929878 + ,-0.69897000433602 + ,0.69897000433602 + ,2.07918124604762 + ,2 + ,2 + ,2 + ,3.3 + ,1.30102999566398 + ,1.44185217577329 + ,2.06069784035361 + ,2.17026171539496 + ,5 + ,5 + ,5 + ,11 + ,0.5910646070265 + ,-0.92081875395238 + ,0 + ,1.20411998265592 + ,3 + ,1 + ,2 + ,4.7 + ,1.61278385671974 + ,1.92941892571429 + ,2.51188336097887 + ,2.49136169383427 + ,1 + ,3 + ,1 + ,10.4 + ,0.95424250943932 + ,-1 + ,0.60205999132796 + ,1.44715803134222 + ,5 + ,1 + ,3 + ,7.4 + ,0.88081359228079 + ,0.01703333929878 + ,0.74036268949424 + ,1.83250891270624 + ,5 + ,3 + ,4 + ,2.1 + ,1.66275783168157 + ,2.71683772329952 + ,2.81624129999178 + ,2.52633927738984 + ,5 + ,5 + ,5 + ,17.9 + ,1.38021124171161 + ,-2 + ,-0.60205999132796 + ,1.69897000433602 + ,1 + ,1 + ,1 + ,6.1 + ,2 + ,1.79239168949825 + ,3.12057393120585 + ,2.42651126136458 + ,1 + ,1 + ,1 + ,11.9 + ,0.50514997831991 + ,-1.69897000433602 + ,-0.39794000867204 + ,1.27875360095283 + ,4 + ,1 + ,3 + ,13.8 + ,0.69897000433602 + ,0.23044892137827 + ,0.79934054945358 + ,1.07918124604762 + ,2 + ,1 + ,1 + ,14.3 + ,0.81291335664286 + ,0.54406804435028 + ,1.03342375548695 + ,2.07918124604762 + ,2 + ,1 + ,1 + ,15.2 + ,1.07918124604762 + ,-0.31875876262441 + ,1.19033169817029 + ,2.14612803567824 + ,2 + ,2 + ,2 + ,10 + ,1.30535136944662 + ,1 + ,2.06069784035361 + ,2.23044892137827 + ,4 + ,4 + ,4 + ,11.9 + ,1.11394335230684 + ,0.20951501454263 + ,1.05690485133647 + ,1.23044892137827 + ,2 + ,1 + ,2 + ,6.5 + ,1.43136376415899 + ,2.28330122870355 + ,2.25527250510331 + ,2.06069784035361 + ,4 + ,4 + ,4 + ,7.5 + ,1.25527250510331 + ,0.39794000867204 + ,1.08278537031645 + ,1.49136169383427 + ,5 + ,5 + ,5 + ,10.6 + ,0.67209785793572 + ,-0.55284196865778 + ,0.27875360095283 + ,1.32221929473392 + ,3 + ,1 + ,3 + ,7.4 + ,0.99122607569249 + ,0.62736585659273 + ,1.70243053644553 + ,1.7160033436348 + ,1 + ,1 + ,1 + ,8.4 + ,1.46239799789896 + ,0.83250891270624 + ,2.25285303097989 + ,2.2148438480477 + ,2 + ,3 + ,2 + ,5.7 + ,0.84509804001426 + ,-0.1249387366083 + ,1.0899051114394 + ,2.35218251811136 + ,2 + ,2 + ,2 + ,4.9 + ,0.77815125038364 + ,0.55630250076729 + ,1.32221929473392 + ,2.35218251811136 + ,3 + ,2 + ,3 + ,3.2 + ,1.30102999566398 + ,1.74429298312268 + ,2.24303804868629 + ,2.17897694729317 + ,5 + ,5 + ,5 + ,11 + ,0.65321251377534 + ,-0.045757490560675 + ,0.41497334797082 + ,1.77815125038364 + ,2 + ,1 + ,2 + ,4.9 + ,0.8750612633917 + ,0.30102999566398 + ,1.0899051114394 + ,2.30102999566398 + ,3 + ,1 + ,3 + ,13.2 + ,0.36172783601759 + ,-1 + ,0.39794000867204 + ,1.66275783168157 + ,3 + ,2 + ,2 + ,9.7 + ,1.38021124171161 + ,0.6222140229663 + ,1.76342799356294 + ,2.32221929473392 + ,4 + ,3 + ,4 + ,12.8 + ,0.47712125471966 + ,0.54406804435028 + ,0.5910646070265 + ,1.14612803567824 + ,2 + ,1 + ,1) + ,dim=c(8 + ,39) + ,dimnames=list(c('SWS' + ,'logL' + ,'logWb' + ,'logWbr' + ,'logtg' + ,'P' + ,'S' + ,'D') + ,1:39)) > y <- array(NA,dim=c(8,39),dimnames=list(c('SWS','logL','logWb','logWbr','logtg','P','S','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' > #'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.41879268 0 1.0000000 1.0758879 0.2087285 2 0.07282071 1 0.5812073 0.8674751 0.1803420 3 0.01000000 2 0.5083866 0.8290468 0.1642184 > postscript(file="/var/www/html/rcomp/tmp/1e5t51292594318.ps",horizontal=F,onefile=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/2oes81292594318.ps",horizontal=F,onefile=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 8.95000 -2.65000 2 2.1 5.68125 -3.58125 3 9.1 5.68125 3.41875 4 15.8 11.84000 3.96000 5 5.2 5.68125 -0.48125 6 10.9 11.84000 -0.94000 7 8.3 5.68125 2.61875 8 11.0 8.95000 2.05000 9 3.2 5.68125 -2.48125 10 6.3 11.84000 -5.54000 11 6.6 11.84000 -5.24000 12 9.5 11.84000 -2.34000 13 3.3 5.68125 -2.38125 14 11.0 11.84000 -0.84000 15 4.7 5.68125 -0.98125 16 10.4 8.95000 1.45000 17 7.4 8.95000 -1.55000 18 2.1 5.68125 -3.58125 19 17.9 11.84000 6.06000 20 6.1 5.68125 0.41875 21 11.9 8.95000 2.95000 22 13.8 11.84000 1.96000 23 14.3 11.84000 2.46000 24 15.2 11.84000 3.36000 25 10.0 5.68125 4.31875 26 11.9 11.84000 0.06000 27 6.5 8.95000 -2.45000 28 7.5 8.95000 -1.45000 29 10.6 8.95000 1.65000 30 7.4 11.84000 -4.44000 31 8.4 5.68125 2.71875 32 5.7 5.68125 0.01875 33 4.9 5.68125 -0.78125 34 3.2 5.68125 -2.48125 35 11.0 11.84000 -0.84000 36 4.9 5.68125 -0.78125 37 13.2 11.84000 1.36000 38 9.7 5.68125 4.01875 39 12.8 11.84000 0.96000 > myr <- residuals(m) > myp <- predict(m) > postscript(file="/var/www/html/rcomp/tmp/3oes81292594318.ps",horizontal=F,onefile=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/42o8z1292594318.tab") > > try(system("convert tmp/1e5t51292594318.ps tmp/1e5t51292594318.png",intern=TRUE)) character(0) > try(system("convert tmp/2oes81292594318.ps tmp/2oes81292594318.png",intern=TRUE)) character(0) > try(system("convert tmp/3oes81292594318.ps tmp/3oes81292594318.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.027 0.491 2.244