R : Copyright 2003, The R Development Core Team
Version 1.6.2  (2003-01-10)

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.

Type `demo()' for some demos, `help()' for on-line help, or
`help.start()' for a HTML browser interface to help.
Type `q()' to quit R.

> invisible(options(echo = TRUE))
> library(rpart)
> set.seed(100542)
> 
> target.lrn <- read.table("../lrn/num/472.dat",header=T,colClasses="numeric")
> target.val <- read.table("../val/num/472.dat",header=T,colClasses="numeric")
> target.tst <- read.table("../tst/num/472.dat",header=T,colClasses="numeric")
> 
>  
> y.lrn <- target.lrn[,1]
> y.val <- target.val[,1]
> y.tst <- target.tst[,1]
> y     <- c(y.lrn,y.val,y.tst)
> 
> target.lrn <- read.table("../lrn/num/471.dat",header=T,colClasses="numeric")
> target.val <- read.table("../val/num/471.dat",header=T,colClasses="numeric")
> target.tst <- read.table("../tst/num/471.dat",header=T,colClasses="numeric")
>  
> b.lrn <- target.lrn[,1]
> b.val <- target.val[,1]
> b.tst <- target.tst[,1]
> b     <- c(b.lrn,b.val,b.tst)
> 
> n.lrn <- length(y.lrn)
> n.val <- length(y.val)
> n.tst <- length(y.tst)
> n     <- length(y)
> 
> rm(target.lrn,target.val,target.tst)
> 
> wts <- mat.or.vec(n,1) ; for (i in 1:n.lrn) wts[i]=1
> idx.lrn <- 1:n.lrn
> idx.val <- (n.lrn+1):(n.lrn+n.val)
> idx.tst <- (n.lrn+n.val+1):n
> 
> mod <- read.table("../cty_mod.txt",
+          header=F,colClasses="character",col.names=c("file","feature","type"))
> 
> n.mod <- length(mod$file)
> 
> first.time <- TRUE
> 
> for (i in 1:n.mod) {
+ 
+   fn.lrn <- paste("../lrn/",mod$type[i],"/",mod$file[i],".dat",sep="")
+   fn.val <- paste("../val/",mod$type[i],"/",mod$file[i],".dat",sep="")
+   fn.tst <- paste("../tst/",mod$type[i],"/",mod$file[i],".dat",sep="")
+   print(mod$feature[i])
+ 
+   if (mod$type[i]=="chr") {
+ 
+     f.lrn <- read.table(fn.lrn,
+                header=T,colClasses="character",blank.lines.skip=F)
+     f.val <- read.table(fn.val,
+                header=T,colClasses="character",blank.lines.skip=F)
+     f.tst <- read.table(fn.tst,
+                header=T,colClasses="character",blank.lines.skip=F)
+ 
+     f <- c(f.lrn[,1],f.val[,1],f.tst[,1])
+ 
+     if (mod$feature[i]=="STATE") {
+       f[f=="AS"|f=="DC"|f=="DE"|f=="MA"|f=="ME"|f=="NH"] <- "S1"
+       f[f=="OH"|f=="RI"|f=="VI"|f=="WV"]                 <- "S1"
+       f[f=="AA"|f=="AE"|f=="AP"|f=="CT"|f=="GU"|f=="MD"] <- "S2"
+       f[f=="NJ"|f=="NY"|f=="PA"|f=="VA"|f=="VT"|f=="WY"] <- "S2"
+       f[f=="AK"|f=="UT"|f=="MS"]                         <- "S3"
+       f[f=="NE"|f=="ND"]                                 <- "S4"
+       f[f=="SD"|f=="SC"]                                 <- "S5"
+     } 
+     
+     f <- as.factor(f)
+ 
+     n.lev <- nlevels(f)
+     f.name <- levels(f) 
+     if(n.lev==2) f.name <- c(mod$feature[i],mod$feature[i])
+     print(paste("  nlevels = ",n.lev))
+ 
+     f <- model.matrix(y ~ f - 1)      # Note: Intercept removed.
+     f <- f[,2:n.lev]                  # Note: First dummy deleted.
+     f.name <- f.name[2:n.lev]         # Note: First name deleted.
+ 
+   } else {
+     
+     f.lrn<-read.table(fn.lrn,
+              header=T,colClasses="numeric",blank.lines.skip=F)
+     f.val<-read.table(fn.val,
+              header=T,colClasses="numeric",blank.lines.skip=F)
+     f.tst<-read.table(fn.tst,
+              header=T,colClasses="numeric",blank.lines.skip=F)
+ 
+     f <- c(f.lrn[,1],f.val[,1],f.tst[,1])
+ 
+     f[is.na(f)] <- 0
+ 
+     f.name <- mod$feature[i]
+ 
+     if (mod$feature[i]=="DOB") { 
+       d <- f ; d[d>0] <- 1         # Note: Dummy for missing DOB 
+       f <- cbind(d,f,f^2)          # Note: Quadratic term added to DOB.
+       rm(d)
+       f.name <- c("DOB.0","DOB.1","DOB.2")
+     }
+   }
+ 
+   if (first.time) {
+     X <- f
+     X.names <- f.name
+     first.time <- FALSE
+   } else {
+     X <- cbind(prev.X,f)
+     X.names <- c(prev.X.names,f.name)
+   }
+ 
+   prev.X <- X
+   prev.X.names <- X.names
+ 
+ }
[1] "LASTGIFT"
[1] "PEPSTRFL"
[1] "  nlevels =  2"
[1] "STATE"
[1] "  nlevels =  33"
[1] "RECP3"
[1] "  nlevels =  2"
[1] "DOB"
[1] "MAILCODE"
[1] "  nlevels =  2"
[1] "MHUC2"
[1] "LASTDATE"
[1] "MINRAMNT"
> 
> rm(prev.X,prev.X.names)
> rm(f.lrn,f.val,f)
> 
> dimnames(X) <- list(NULL,X.names)
> 
> b.fit <- lm(b~X,weights=wts)       # Note: Validation sample not in fit.
> 
> bhat <- mat.or.vec(n,1) 
> bhat[b.fit$fitted.values>0.5] <- 1
> 
> ehat <- (b != bhat)
> 
> cer.nul <- mean(b[idx.lrn]) 
> cer.lrn <- mean(ehat[idx.lrn])
> cer.val <- mean(ehat[idx.val])
> cer.tst <- mean(ehat[idx.tst])
> 
> print(paste("  cer.nul = ",cer.nul))
[1] "  cer.nul =  0.0505530642750374"
> print(paste("  cer.lrn = ",cer.lrn))
[1] "  cer.lrn =  0.0505530642750374"
> print(paste("  cer.val = ",cer.val))
[1] "  cer.val =  0.0510037409768692"
> print(paste("  cer.tst = ",cer.tst))
[1] "  cer.tst =  0.0517150949333893"
> 
> c <- as.factor(b)
> 
> ctrl <- rpart.control(minsplit=20,maxsurrogate=0,maxdepth=25,cp=0.0001,xval=0)
> tree <- rpart(c~X,weights=wts,method="class",control=ctrl)
> 
> chat <- tree$frame[tree$where,5] - 1
> 
> c <- b
> ehat <- (c != chat)
> 
> cer.nul <- mean(b[idx.lrn])
> cer.lrn <- mean(ehat[idx.lrn])
> cer.val <- mean(ehat[idx.val])
> cer.tst <- mean(ehat[idx.tst])  
> 
> print(paste("  cer.nul = ",cer.nul))
[1] "  cer.nul =  0.0505530642750374"
> print(paste("  cer.lrn = ",cer.lrn))
[1] "  cer.lrn =  0.0503288490284006"
> print(paste("  cer.val = ",cer.val))
[1] "  cer.val =  0.0514779493123979"
> print(paste("  cer.tst = ",cer.tst))
[1] "  cer.tst =  0.0525542851148642"
> 
> y.fit <- lm(y~X,weights=wts)       # Note: Validation sample not in fit.
> 
> yhat <- y.fit$fitted.values
> 
> ehat <- y - mean(y)
> mse.nul <- mean(ehat[idx.lrn]^2) 
> 
> ehat <- y.fit$residuals
> mse.lrn <- mean(ehat[idx.lrn]^2)
> mse.val <- mean(ehat[idx.val]^2)
> mse.tst <- mean(ehat[idx.tst]^2)
> 
> print(paste("  mse.nul = ",mse.nul))
[1] "  mse.nul =  20.0992375595588"
> print(paste("  mse.lrn = ",mse.lrn))
[1] "  mse.lrn =  19.9608334474451"
> print(paste("  mse.val = ",mse.val))
[1] "  mse.val =  18.6770872778876"
> print(paste("  mse.tst = ",mse.tst))
[1] "  mse.tst =  17.8000310256753"
> 
> x0 <- c(0,n.lrn)/n.lrn
> y0 <- c(0,sum(y.lrn-0.68))
> 
> idx <- seq(1,n.lrn,length=500)
> 
> x1 <- (1:n.lrn)/n.lrn
> x1 <- x1[idx]
> 
> idx.yhat <- sort.list(yhat[idx.lrn],decreasing=TRUE)
> y1 <- y.lrn[idx.yhat]-0.68
> y1 <- cumsum(y1)
> y1 <- y1[idx]
> 
> idx.bhat <- sort.list(bhat[idx.lrn],decreasing=TRUE)  
> y2 <- y.lrn[idx.bhat]-0.68
> y2 <- cumsum(y2)
> y2 <- y2[idx]
> 
> idx.chat <- sort.list(chat[idx.lrn],decreasing=TRUE)  
> y3 <- y.lrn[idx.chat]-0.68
> y3 <- cumsum(y3)
> y3 <- y3[idx]
> 
> source("psopts.r");
> postscript(file="cty_class_lif_lrn.eps");
> 
> plot(x=c(x0,x1),y=c(y0,y1),ylab="dollars",xlab="proportion",type="n")
> lines(x=x0,y=y0,col="green")
> lines(x=x1,y=y1,col="blue")
> lines(x=x1,y=y2,col="red")
> lines(x=x1,y=y3,col="orange")
> 
> dev.off()
null device 
          1 
> 
> x0 <- c(0,n.val)/n.val
> y0 <- c(0,sum(y.val-0.68))
> 
> idx <- seq(1,n.val,length=500)
> 
> x1 <- (1:n.val)/n.val
> x1 <- x1[idx]
> 
> idx.yhat <- sort.list(yhat[idx.val],decreasing=TRUE)
> y1 <- y.val[idx.yhat]-0.68
> y1 <- cumsum(y1)
> y1 <- y1[idx]
> 
> idx.bhat <- sort.list(bhat[idx.val],decreasing=TRUE)  
> y2 <- y.val[idx.bhat]-0.68
> y2 <- cumsum(y2)
> y2 <- y2[idx]
> 
> idx.chat <- sort.list(chat[idx.val],decreasing=TRUE)  
> y3 <- y.val[idx.chat]-0.68
> y3 <- cumsum(y3)
> y3 <- y3[idx]
> 
> source("psopts.r");
> postscript(file="cty_class_lif_val.eps");
> 
> plot(x=c(x0,x1),y=c(y0,y1),ylab="dollars",xlab="proportion",type="n")
> lines(x=x0,y=y0,col="green")
> lines(x=x1,y=y1,col="blue")
> lines(x=x1,y=y2,col="red")
> lines(x=x1,y=y3,col="orange")
> 
> dev.off()
null device 
          1 
> 
> source("psopts.r"); 
> postscript(file="cty_class_tree.eps")
> 
> plot.rpart(tree,uniform=T,margin=0.01)
> #text.rpart(tree)
> 
> dev.off()
null device 
          1 
> 
> proc.time()
[1] 121.51   7.07 128.79   0.00   0.00
>