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(0.30102999566398 + ,0.65321251377534 + ,0 + ,0.81954393554187 + ,1.6232492903979 + ,3 + ,1 + ,3 + ,0.25527250510331 + ,1.83884909073726 + ,3.40602894496362 + ,3.66304097489397 + ,2.79518458968242 + ,3 + ,5 + ,4 + ,-0.15490195998574 + ,1.43136376415899 + ,1.02325245963371 + ,2.25406445291434 + ,2.25527250510331 + ,4 + ,4 + ,4 + ,0.5910646070265 + ,1.27875360095283 + ,-1.69897000433602 + ,-0.52287874528034 + ,1.54406804435028 + ,1 + ,1 + ,1 + ,0 + ,1.48287358360875 + ,2.20411998265592 + ,2.22788670461367 + ,2.59328606702046 + ,4 + ,5 + ,4 + ,0.55630250076729 + ,1.44715803134222 + ,0.51851393987789 + ,1.40823996531185 + ,1.79934054945358 + ,1 + ,2 + ,1 + ,0.14612803567824 + ,1.69897000433602 + ,1.71733758272386 + ,2.64345267648619 + ,2.36172783601759 + ,1 + ,1 + ,1 + ,0.17609125905568 + ,0.84509804001426 + ,-0.36653154442041 + ,0.80617997398389 + ,2.04921802267018 + ,5 + ,4 + ,4 + ,-0.15490195998574 + ,1.47712125471966 + ,2.66745295288995 + ,2.62634036737504 + ,2.44870631990508 + ,5 + ,5 + ,5 + ,0.32221929473392 + ,0.54406804435028 + ,-1.09691001300806 + ,0.079181246047625 + ,1.6232492903979 + ,1 + ,1 + ,1 + ,0.61278385671974 + ,0.77815125038364 + ,-0.10237290870956 + ,0.54406804435028 + ,1.6232492903979 + ,2 + ,2 + ,2 + ,0.079181246047625 + ,1.01703333929878 + ,-0.69897000433602 + ,0.69897000433602 + ,2.07918124604762 + ,2 + ,2 + ,2 + ,-0.30102999566398 + ,1.30102999566398 + ,1.44185217577329 + ,2.06069784035361 + ,2.17026171539496 + ,5 + ,5 + ,5 + ,0.53147891704226 + ,0.5910646070265 + ,-0.92081875395238 + ,0 + ,1.20411998265592 + ,3 + ,1 + ,2 + ,0.17609125905568 + ,1.61278385671974 + ,1.92941892571429 + ,2.51188336097887 + ,2.49136169383427 + ,1 + ,3 + ,1 + ,0.53147891704226 + ,0.95424250943932 + ,-1 + ,0.60205999132796 + ,1.44715803134222 + ,5 + ,1 + ,3 + ,-0.096910013008056 + ,0.88081359228079 + ,0.01703333929878 + ,0.74036268949424 + ,1.83250891270624 + ,5 + ,3 + ,4 + ,-0.096910013008056 + ,1.66275783168157 + ,2.71683772329952 + ,2.81624129999178 + ,2.52633927738984 + ,5 + ,5 + ,5 + ,0.30102999566398 + ,1.38021124171161 + ,-2 + ,-0.60205999132796 + ,1.69897000433602 + ,1 + ,1 + ,1 + ,0.27875360095283 + ,2 + ,1.79239168949825 + ,3.12057393120585 + ,2.42651126136458 + ,1 + ,1 + ,1 + ,0.11394335230684 + ,0.50514997831991 + ,-1.69897000433602 + ,-0.39794000867204 + ,1.27875360095283 + ,4 + ,1 + ,3 + ,0.7481880270062 + ,0.69897000433602 + ,0.23044892137827 + ,0.79934054945358 + ,1.07918124604762 + ,2 + ,1 + ,1 + ,0.49136169383427 + ,0.81291335664286 + ,0.54406804435028 + ,1.03342375548695 + ,2.07918124604762 + ,2 + ,1 + ,1 + ,0.25527250510331 + ,1.07918124604762 + ,-0.31875876262441 + ,1.19033169817029 + ,2.14612803567824 + ,2 + ,2 + ,2 + ,-0.045757490560675 + ,1.30535136944662 + ,1 + ,2.06069784035361 + ,2.23044892137827 + ,4 + ,4 + ,4 + ,0.25527250510331 + ,1.11394335230684 + ,0.20951501454263 + ,1.05690485133647 + ,1.23044892137827 + ,2 + ,1 + ,2 + ,0.27875360095283 + ,1.43136376415899 + ,2.28330122870355 + ,2.25527250510331 + ,2.06069784035361 + ,4 + ,4 + ,4 + ,-0.045757490560675 + ,1.25527250510331 + ,0.39794000867204 + ,1.08278537031645 + ,1.49136169383427 + ,5 + ,5 + ,5 + ,0.41497334797082 + ,0.67209785793572 + ,-0.55284196865778 + ,0.27875360095283 + ,1.32221929473392 + ,3 + ,1 + ,3 + ,0.38021124171161 + ,0.99122607569249 + ,0.62736585659273 + ,1.70243053644553 + ,1.7160033436348 + ,1 + ,1 + ,1 + ,0.079181246047625 + ,1.46239799789896 + ,0.83250891270624 + ,2.25285303097989 + ,2.2148438480477 + ,2 + ,3 + ,2 + ,-0.045757490560675 + ,0.84509804001426 + ,-0.1249387366083 + ,1.0899051114394 + ,2.35218251811136 + ,2 + ,2 + ,2 + ,-0.30102999566398 + ,0.77815125038364 + ,0.55630250076729 + ,1.32221929473392 + ,2.35218251811136 + ,3 + ,2 + ,3 + ,-0.22184874961636 + ,1.30102999566398 + ,1.74429298312268 + ,2.24303804868629 + ,2.17897694729317 + ,5 + ,5 + ,5 + ,0.36172783601759 + ,0.65321251377534 + ,-0.045757490560675 + ,0.41497334797082 + ,1.77815125038364 + ,2 + ,1 + ,2 + ,-0.30102999566398 + ,0.8750612633917 + ,0.30102999566398 + ,1.0899051114394 + ,2.30102999566398 + ,3 + ,1 + ,3 + ,0.41497334797082 + ,0.36172783601759 + ,-1 + ,0.39794000867204 + ,1.66275783168157 + ,3 + ,2 + ,2 + ,-0.22184874961636 + ,1.38021124171161 + ,0.6222140229663 + ,1.76342799356294 + ,2.32221929473392 + ,4 + ,3 + ,4 + ,0.81954393554187 + ,0.47712125471966 + ,0.54406804435028 + ,0.5910646070265 + ,1.14612803567824 + ,2 + ,1 + ,1) + ,dim=c(8 + ,39) + ,dimnames=list(c('logPS' + ,'logL' + ,'logWb' + ,'logWbr' + ,'logtg' + ,'P' + ,'S' + ,'D') + ,1:39)) > y <- array(NA,dim=c(8,39),dimnames=list(c('logPS','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.48752266 0 1.0000000 1.0531531 0.1802172 2 0.09264545 1 0.5124773 0.6774661 0.1240291 3 0.01000000 2 0.4198319 0.6895767 0.1059018 > postscript(file="/var/www/html/rcomp/tmp/1ji2u1292594548.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/2ji2u1292594548.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 0.30103000 0.27895288 0.0220771168 2 0.25527251 -0.05684936 0.3121218647 3 -0.15490196 -0.05684936 -0.0980526004 4 0.59106461 0.52624016 0.0648244450 5 0.00000000 -0.05684936 0.0568493596 6 0.55630250 0.52624016 0.0300623387 7 0.14612804 -0.05684936 0.2029773953 8 0.17609126 0.27895288 -0.1028616198 9 -0.15490196 -0.05684936 -0.0980526004 10 0.32221929 0.52624016 -0.2040208673 11 0.61278386 0.27895288 0.3338309778 12 0.07918125 0.27895288 -0.1997716328 13 -0.30103000 -0.05684936 -0.2441806361 14 0.53147892 0.27895288 0.2525260381 15 0.17609126 -0.05684936 0.2329406186 16 0.53147892 0.27895288 0.2525260381 17 -0.09691001 0.27895288 -0.3758628919 18 -0.09691001 -0.05684936 -0.0400606534 19 0.30103000 0.52624016 -0.2252101664 20 0.27875360 -0.05684936 0.3356029605 21 0.11394335 0.27895288 -0.1650095266 22 0.74818803 0.52624016 0.2219478650 23 0.49136169 0.52624016 -0.0348784682 24 0.25527251 0.27895288 -0.0236803738 25 -0.04575749 -0.05684936 0.0110918690 26 0.25527251 0.27895288 -0.0236803738 27 0.27875360 0.27895288 -0.0001992779 28 -0.04575749 0.27895288 -0.3247103695 29 0.41497335 0.27895288 0.1360204691 30 0.38021124 0.52624016 -0.1460289203 31 0.07918125 -0.05684936 0.1360306056 32 -0.04575749 -0.05684936 0.0110918690 33 -0.30103000 -0.05684936 -0.2441806361 34 -0.22184875 -0.05684936 -0.1649993900 35 0.36172784 0.27895288 0.0827749571 36 -0.30103000 -0.05684936 -0.2441806361 37 0.41497335 0.27895288 0.1360204691 38 -0.22184875 -0.05684936 -0.1649993900 39 0.81954394 0.52624016 0.2933037735 > myr <- residuals(m) > myp <- predict(m) > postscript(file="/var/www/html/rcomp/tmp/3t9kx1292594548.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/4pjh51292594548.tab") > > try(system("convert tmp/1ji2u1292594548.ps tmp/1ji2u1292594548.png",intern=TRUE)) character(0) > try(system("convert tmp/2ji2u1292594548.ps tmp/2ji2u1292594548.png",intern=TRUE)) character(0) > try(system("convert tmp/3t9kx1292594548.ps tmp/3t9kx1292594548.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.029 0.478 2.233