R version 2.12.0 (2010-10-15) Copyright (C) 2010 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-pc-linux-gnu (32-bit) 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(235.1,9700,280.7,9081,264.6,9084,240.7,9743,201.4,8587,240.8,9731,241.1,9563,223.8,9998,206.1,9437,174.7,10038,203.3,9918,220.5,9252,299.5,9737,347.4,9035,338.3,9133,327.7,9487,351.6,8700,396.6,9627,438.8,8947,395.6,9283,363.5,8829,378.8,9947,357,9628,369,9318,464.8,9605,479.1,8640,431.3,9214,366.5,9567,326.3,8547,355.1,9185,331.6,9470,261.3,9123,249,9278,205.5,10170,235.6,9434,240.9,9655,264.9,9429,253.8,8739,232.3,9552,193.8,9687,177,9019,213.2,9672,207.2,9206,180.6,9069,188.6,9788,175.4,10312,199,10105,179.6,9863,225.8,9656,234,9295,200.2,9946,183.6,9701,178.2,9049,203.2,10190,208.5,9706,191.8,9765,172.8,9893,148,9994,159.4,10433,154.5,10073,213.2,10112,196.4,9266,182.8,9820,176.4,10097,153.6,9115,173.2,10411,171,9678,151.2,10408,161.9,10153,157.2,10368,201.7,10581,236.4,10597,356.1,10680,398.3,9738,403.7,9556),dim=c(2,75),dimnames=list(c('unemployment','birth'),1:75)) > y <- array(NA,dim=c(2,75),dimnames=list(c('unemployment','birth'),1:75)) > 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.26056479 0 1.0000000 1.0276889 0.1496079 2 0.05406311 1 0.7394352 0.7677842 0.1273217 3 0.02706758 3 0.6313090 0.8616137 0.1515033 4 0.02195487 4 0.6042414 0.8747003 0.1494317 5 0.01000000 5 0.5822865 0.8237284 0.1434389 > postscript(file="/var/www/rcomp/tmp/1i4yx1293205590.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/rcomp/tmp/2i4yx1293205590.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 235.1 237.6000 -2.5000000 2 280.7 226.0143 54.6857143 3 264.6 226.0143 38.5857143 4 240.7 237.6000 3.1000000 5 201.4 344.9286 -143.5285714 6 240.8 237.6000 3.2000000 7 241.1 346.8111 -105.7111111 8 223.8 198.7962 25.0038462 9 206.1 283.1643 -77.0642857 10 174.7 198.7962 -24.0961538 11 203.3 198.7962 4.5038462 12 220.5 283.1643 -62.6642857 13 299.5 237.6000 61.9000000 14 347.4 226.0143 121.3857143 15 338.3 283.1643 55.1357143 16 327.7 346.8111 -19.1111111 17 351.6 344.9286 6.6714286 18 396.6 346.8111 49.7888889 19 438.8 344.9286 93.8714286 20 395.6 283.1643 112.4357143 21 363.5 344.9286 18.5714286 22 378.8 198.7962 180.0038462 23 357.0 346.8111 10.1888889 24 369.0 283.1643 85.8357143 25 464.8 346.8111 117.9888889 26 479.1 344.9286 134.1714286 27 431.3 283.1643 148.1357143 28 366.5 346.8111 19.6888889 29 326.3 344.9286 -18.6285714 30 355.1 283.1643 71.9357143 31 331.6 346.8111 -15.2111111 32 261.3 283.1643 -21.8642857 33 249.0 283.1643 -34.1642857 34 205.5 198.7962 6.7038462 35 235.6 283.1643 -47.5642857 36 240.9 237.6000 3.3000000 37 264.9 283.1643 -18.2642857 38 253.8 344.9286 -91.1285714 39 232.3 346.8111 -114.5111111 40 193.8 237.6000 -43.8000000 41 177.0 226.0143 -49.0142857 42 213.2 237.6000 -24.4000000 43 207.2 283.1643 -75.9642857 44 180.6 226.0143 -45.4142857 45 188.6 198.7962 -10.1961538 46 175.4 198.7962 -23.3961538 47 199.0 198.7962 0.2038462 48 179.6 198.7962 -19.1961538 49 225.8 237.6000 -11.8000000 50 234.0 283.1643 -49.1642857 51 200.2 198.7962 1.4038462 52 183.6 237.6000 -54.0000000 53 178.2 226.0143 -47.8142857 54 203.2 198.7962 4.4038462 55 208.5 237.6000 -29.1000000 56 191.8 198.7962 -6.9961538 57 172.8 198.7962 -25.9961538 58 148.0 198.7962 -50.7961538 59 159.4 198.7962 -39.3961538 60 154.5 198.7962 -44.2961538 61 213.2 198.7962 14.4038462 62 196.4 283.1643 -86.7642857 63 182.8 198.7962 -15.9961538 64 176.4 198.7962 -22.3961538 65 153.6 226.0143 -72.4142857 66 173.2 198.7962 -25.5961538 67 171.0 237.6000 -66.6000000 68 151.2 198.7962 -47.5961538 69 161.9 198.7962 -36.8961538 70 157.2 198.7962 -41.5961538 71 201.7 198.7962 2.9038462 72 236.4 198.7962 37.6038462 73 356.1 198.7962 157.3038462 74 398.3 237.6000 160.7000000 75 403.7 346.8111 56.8888889 > myr <- residuals(m) > myp <- predict(m) > postscript(file="/var/www/rcomp/tmp/3i4yx1293205590.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/rcomp/tmp/4eew51293205590.tab") > > try(system("convert tmp/1i4yx1293205590.ps tmp/1i4yx1293205590.png",intern=TRUE)) character(0) > try(system("convert tmp/2i4yx1293205590.ps tmp/2i4yx1293205590.png",intern=TRUE)) character(0) > try(system("convert tmp/3i4yx1293205590.ps tmp/3i4yx1293205590.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.290 0.570 1.851