###################################################
### chunk number 1: 
###################################################
options(show.signif.stars=FALSE,digits=4,width=80)


###################################################
### chunk number 2: 
###################################################
loc <- "http://www.stat.umn.edu/~sandy/courses/8053/Data/diamonds.txt"
diamonds <- read.table(url(loc), header=TRUE)
dim(diamonds)


###################################################
### chunk number 3: 
###################################################
set.seed(10131985)
est.set <- sample(1:308, 200)
val.set <- (1:308)[-est.set]
diamonds1 <- diamonds[est.set, ]


###################################################
### chunk number 4: one
###################################################
library(rpart)
r1 <- rpart(Price ~., diamonds1, method = "anova")
par(mfrow=c(1, 2), xpd = NA)
plot(r1, margin=.1)
text(r1)
plot(r1, compress=TRUE, uniform=TRUE, branch=0.4, margin=.1)
text(r1)


###################################################
### chunk number 5: two
###################################################
m1 <- lm(Price~., diamonds1)
plot(predict(m1) ~ predict(r1), xlab="Tree based", 
     ylab="ols based", las=2)
abline(0, 1, col="gray", lwd=2)


###################################################
### chunk number 6: three
###################################################
r2 <- update(r1, cp=0.001, xval = 10)
printcp(r2)
plotcp(r2)


###################################################
### chunk number 7: threea
###################################################
par(mfrow=c(1, 2), xpd = NA)
r3 <- prune(r2, cp=r2$cptable[7, 1])
plot(r1, margin=.1)
text(r1)
plot(r3, margin=.1)
text(r3)


###################################################
### chunk number 8: 
###################################################
get.error <- function(m, which){
with(diamonds, sqrt(sum((Price[which] -
      predict(m, diamonds[which, ]))^2)/length(which)))}
summary(m2 <- step(m1, trace=0))$coef
(ans <- data.frame(rpart=c(get.error(r3, est.set), c(get.error(r3,val.set))),
                  ols = c(get.error(m1, est.set), c(get.error(m1, val.set))),
                  ols.step= c(get.error(m2, est.set), c(get.error(m2, val.set))),
                  row.names=c("Estimation", "Validation")))


###################################################
### chunk number 9: four
###################################################
library(alr3)
par(mfrow=c(1, 2), xpd=NA)
set.seed(123456)
const <- rep(FALSE, dim(blowdown)[1])
const[sample(1:dim(blowdown)[1], 2000)] <- TRUE
blowdown$y1 <- factor(blowdown$y)
r3 <- rpart(y1 ~ S + D + SPP, blowdown, subset=const, method="class")
plot(r3,compress=F, uniform=T, branch=0.5, margin=.2)
text(r3)
r4 <- update(r3, cp = 0.001)
plotcp(r4)
printcp(r4)


###################################################
### chunk number 10: 
###################################################
show <- function(tt){
  print(tt)
  cat(paste("Misclassification rate =",round(1-sum(diag(tt))/sum(tt), 2),"\n"))
  invisible()}
show(with(blowdown, table(actual=y1[const],
         predicted=predict(r3, type="class"))))
show(with(blowdown, table(actual=y1[!const],
         predicted=predict(r3, newdata=blowdown[!const,], type="class"))))


###################################################
### chunk number 11: 
###################################################
library(mgcv)
m3 <- gam(y ~ s(S) + log2(D) + SPP, data=blowdown, family=binomial, subset=const)
show(with(blowdown,table(actual=y[const],
          predicted=ifelse(predict(m3) <= 0, 0 ,1))))
show(with(blowdown,table(actual=y[!const],
          predicted=ifelse(predict(m3, newdata=blowdown[!const,]) <= 0, 0, 1))))
m4 <- glm(y ~ S + log2(D) + SPP, data=blowdown, family=binomial, subset=const)
show(with(blowdown, table(actual=y[const],
          predicted=ifelse(predict(m4)<= 0, 0, 1))))
show(with(blowdown,table(actual=y[!const],
          predicted=ifelse(predict(m4,newdata=blowdown[!const,]) <= 0, 0, 1))))


###################################################
### chunk number 12: 
###################################################
library(ipred)
b1 <- bagging(y1 ~ S + D + SPP, blowdown, subset=const,nbagg=100)
show(with(blowdown, table(actual=y[const],predicted=predict(b1))))
show(with(blowdown,table(actual=y[!const],
     predicted=predict(b1, newdata=blowdown[!const,]))))


###################################################
### chunk number 13: 
###################################################
library(randomForest)
b2 <- randomForest(y1 ~ S + D + SPP, blowdown, subset=const)
show(with(blowdown,table(actual=y1[const], predicted=predict(b2))))
show(with(blowdown,table(actual=y1[!const],
     predicted=predict(b2, newdata=blowdown
[!const,]))))


###################################################
### chunk number 14: fivea
###################################################
varImpPlot(b2,type=2)


###################################################
### chunk number 15:  eval=FALSE
###################################################
## source("Aaron/R/neteffectplots.R")
## source("Aaron/R/sliceplot.R")
## source("Aaron/R/sliceplot.R")
## source("Aaron/R/car.R")
## source("Aaron/ARMS/R/ARMS")
## source("Aaron/ARMS/R/arms2.R")
## showallpd<-function(model,var,ylim=NULL,...) {
##         pd1<-pd.plot(model,var,display=FALSE,ylim=ylim,...)
## #       pd2<-pd.plot(model,var, all.lines=TRUE,ylim=range(pd1$y),ylim=ylim,...)
##         pd2<-pd.plot(model,var, all.lines=TRUE,ylim=ylim)
##         points(pd1$x,pd1$y)
##         pd2
## }


###################################################
### chunk number 16: bh eval=FALSE
###################################################
## library(mgcv)
## islands <- read.csv("Aaron/data/islands.csv",header=TRUE)
## islands <- na.omit(islands)
## islands$Def <- ifelse(islands$Deforestation > 3.5, "High", "Low")
## gam1<-gam(Deforestation ~ s(Rainfall) + s(Latitude) + Age + Tephra + s(Dust) +
##     s(Elevation) + s(log(Area)) + s(Isolation25),data=islands)
## rf2<-mf<-randomForest(Deforestation ~ Rainfall + Latitude + Age + Tephra + Dust+
##     Elevation + Area + Isolation25, data=islands, na.action=na.omit)
## par(mfrow=c(1,2))
## plot(gam1,residuals=TRUE,cex=5,lwd=2,select=1,main="GAM Partial Residual",
##   ylim=c(-1.5,1.5))
## foo<-showallpd(rf2,"Rainfall",ylim=c(-1.5,1.5),main="Partial Dependence")
## mtext("Random Forest",line=1.4,cex=1.2)
## #partialPlot(rf2,islands,"Rainfall")


###################################################
### chunk number 17: five eval=FALSE
###################################################
## par(mfrow=c(1,2), xpd=NA)
## partialPlot(b2,blowdown,D)
## partialPlot(b2,blowdown,S)