\documentclass{article}
\usepackage[utf8]{inputenc}

\begin{document}
%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{Tree-based regression (Set 8)}

\title{8: Tree-based regression}
\author{John H Maindonald}
\maketitle

<<setup, cache=FALSE, echo=FALSE, include=FALSE>>=
library(knitr)
options(replace.assign=FALSE,width=72)
opts_chunk$set(fig.path='figs/tree-', cache.path='cache/tree-',
               fig.align='center', dev='pdf', fig.width=3.5,
               fig.height=3.5, out.width="0.8\\textwidth",
               fig.show='hold', par=TRUE,
               tidy=FALSE,  comment=NA)
knit_hooks$set(pars=function(before, options, envir){
if (before && options$fig.show!='none') par(mar=c(4,4,1.6,.1),
              font.main=1, cex.lab=.95,cex.axis=.9,mgp=c(2,.7,0),
              tcl=-.3)
              par(options$pars)
}, crop=hook_pdfcrop)
pdf.options(pointsize=12)
oldopt <- options(digits=4)
@ %

\subsection*{Ideas and issues illustrated by the graphs in this vignette}

The fitting of a tree proceeds by making a succession of splits on the
$x$-variable or variables. For tree-based regression, the splitting
criterion is named \texttt{Anova}.  (To explicitly request use of this
criterion, specify \texttt{method="anova"} in the call to
\texttt{rpart()}.)

<<figControl>>=
# To include the figures, change `showFigs <- FALSE`  
# to `showFigs <- TRUE` in the source `.Rnw` file,
# and regenerate the PDF.
#
showFigs <- FALSE
@ %

\section{Code for the Figures}

<<fig8_1A, eval=TRUE, echo=TRUE>>=
fig8.1A <- function(){
  if(!exists('car90.rpart'))
  car90.rpart <- rpart(Mileage ~ tonsWt, data=Car90)
  plot(car90.rpart)
text(car90.rpart, xpd=TRUE, digits=3)
mtext(side=3, line=1.25, "A: Regression tree", adj=0)
}
@ %

<<fig8_1B, eval=TRUE, echo=TRUE>>=
fig8.1B <- function(){
if(!exists('car90.rpart'))
  car90.rpart <- rpart(Mileage ~ tonsWt, data=Car90)
plot(Mileage ~ tonsWt, data=Car90)
wt <- with(Car90, tonsWt)
hat <- predict(car90.rpart)
addhlines(wt, hat, lwd=2, col="gray")
mtext(side=3, line=1.25, "B: Predicted values from tree", adj=0)
}
@ %

<<fig8_2, eval=TRUE, echo=TRUE>>=
fig8.2 <- function(){
BSS <- bssBYcut(tonsWt, Mileage, Car90)
with(BSS, plot(xOrd, bss, xlab="Cutpoint",
               ylab="Between groups sum of squares"))
abline(v=1.218, lty=2)
}
@ %

<<fig8_3A, eval=TRUE, echo=TRUE>>=
fig8.3A <- function(){
opar <- par(mar=c(4,4,2.6,1.6))
if(!exists('car90x.rpart'))
  car90x.rpart <- rpart(Mileage ~ tonsWt, data=Car90,
                        minbucket=5, minsplit=10,
                        cp=0.001)
plot(car90x.rpart, uniform=TRUE)
text(car90x.rpart, digits=3, xpd=TRUE)
mtext(side=3, line=0.75, "A: Decision tree", adj=0)
par(opar)
}
@ %

<<fig8_3B, eval=TRUE, echo=TRUE>>=
opar <- par(mar=c(4,4,2.6,1.6))
fig8.3B <- function(){
if(!exists('car90x.rpart'))
  car90x.rpart <- rpart(Mileage ~ tonsWt, data=Car90,
                        minbucket=5, minsplit=10,
                        cp=0.001)
plot(Mileage ~ tonsWt, data=Car90)
hat <- predict(car90x.rpart)
wt <- with(Car90, tonsWt)
addhlines(wt, hat, lwd=2, col="gray")
mtext(side=3, line=0.75, "B: Mileage vs tonsWt", adj=0)
par(opar)
}
@ %

<<fig8_4, eval=TRUE, echo=TRUE>>=
fig8.4 <- function(){
if(!exists('car90x.rpart'))
  car90x.rpart <- rpart(Mileage ~ tonsWt, data=Car90,
                        minbucket=5, minsplit=10,
                        cp=0.001)
plotcp(car90x.rpart)
}
@ %

<<fig8_5, eval=TRUE, echo=TRUE>>=
fig8.5 <- function(){
if(!exists('car90.rf'))
  car90.rf <- randomForest(Mileage ~ tonsWt,
                          data=Car90)
plot(Mileage ~ tonsWt, data=Car90, type="n")
with(Car90, points(Mileage ~ tonsWt, cex=0.8))
hat <- predict(car90.rf)
with(Car90, points(hat ~ tonsWt, pch="-"))
}
@ %

<<fig8_6, eval=TRUE, echo=TRUE>>=
fig8.6 <- function(){
ran <- range(errsmat)
at <- round(ran+c(0.02,-0.02)*diff(ran),2)
lis <- list(limits=ran, at=at, labels=format(at, digits=2))
lims=list(lis,lis,lis,lis,lis,lis)
library(lattice)
splom(errsmat,
      pscales=lims,
      par.settings=simpleTheme(cex=0.75),
      col=adjustcolor("black", alpha=0.5),
      panel=function(x,y,...){lpoints(x,y,...)
      panel.abline(0,1,col="gray")}
)
}
@ %

\section{Show the Figures}

<<pkgs-figs8, eval=showFigs, message=FALSE, warning=FALSE>>=
pkgs <- c("rpart","mgcv","randomForest","gamclass")
z <- sapply(pkgs, require, character.only=TRUE, warn.conflicts=FALSE)
if(any(!z)){
  notAvail <- paste(names(z)[!z], collapse=", ")
  print(paste("The following packages should be installed:", notAvail))
}
@ %

<<Car90, eval=showFigs>>=
if(!exists('Car90'))
Car90 <- na.omit(car90[, c("Mileage","Weight")])
## Express weight in metric tonnes
Car90 <- within(Car90, tonsWt <- Weight/2240)
@ %

<<meuse>>=
getmeuse <- function(){
if(require('sp', quietly=TRUE)){
data("meuse", package="sp", envir = environment())
meuse <- within(meuse, {levels(soil) <- c("1","2","2")
                        ffreq <- as.numeric(ffreq)
                        loglead <- log(lead)
                        })
invisible(meuse)
} else if(!exists("meuse"))
  print("Dataset 'meuse' was not found, get from package 'sp'")
}
@ %

<<meuse-rf>>=
cfRF <- function(nrep=50){
form1 <- ~ dist + elev + soil + ffreq
form3 <- ~ s(dist, k=3) + s(elev,k=3) + soil +ffreq
form3x <- ~ s(dist, k=3) + s(elev,k=3) + s(x, k=3) + soil+ffreq
form8x <- ~ s(dist, k=8) + s(elev,k=8) + s(x, k=8) + soil+ffreq
formlist <- list("Hybrid1"=form1, "Hybrid3"=form3,
                 "Hybrid3x"=form3x, "Hybrid8x"=form8x)
## ----rfgam-setup----
rfVars <- c("dist", "elev", "soil", "ffreq", "x", "y")
errsmat <- matrix(0, nrep, length(formlist)+2)
dimnames(errsmat)[[2]] <- c(names(formlist), "rfTest", "rfOOB")
n <- 95
for(i in 1:nrep){
sub <- sample(1:nrow(meuse), n)
meuseOut <- meuse[-sub,]
meuseIn <- meuse[sub,]
errsmat[i, ] <- gamRF(formlist=formlist, yvar="loglead",
                      rfVars=rfVars,
                      data=meuseIn, newdata=meuseOut,
                      printit=FALSE)
}
invisible(errsmat)
}
@ %

\begin{figure}[h]
<<car90-plot12, fig.width=3.35, fig.height=3, out.width="0.47\\textwidth", eval=showFigs, echo=FALSE>>=
fig8.1A()
fig8.1B()
@ %
\caption{Regression tree for predicting \texttt{Mileage} given
\texttt{Weight}. At each node, observations for which the criterion
is satisfied take the branch to the left.  Thus at the
first node, \texttt{tonsWt$\geq$1.218} chooses the branch to the
left, while \texttt{tonsWt$<$1.218} chooses the branch to the
right.  Panel B plots \texttt{Mileage} versus \texttt{tonsWt}, with
fitted values
from the \texttt{rpart} model shown as horizontal grey lines.}
\end{figure}

\begin{figure}
\begin{minipage}[c]{0.65\linewidth}
<<bss-plot, echo=FALSE, , eval=showFigs, fig.width=3.35, fig.height=3, out.width="0.8\\textwidth">>=
fig8.2()
@ %
\end{minipage}
\hspace*{0.05\linewidth}
\begin{minipage}[c]{0.28\linewidth}
\caption{Between group sum of squares for \texttt{Mileage}, as a
function of the value of \texttt{tonsWt} at which the split is
made. The choice $c = 1.218$ maximizes the between groups
sum of squares.\label{fig:splitbss}}
\end{minipage}
\vspace*{-18pt}
\end{figure}

\begin{figure}[h]
<<Car90-loosen, echo=FALSE, eval=showFigs, fig.width=3.5, fig.height=2.75, out.width="0.47\\textwidth">>=
fig8.3A()
fig8.3B()
@ %
\caption{For the decision tree to which these panels relate,
the minimum number at each terminal leaf (\texttt{minbucket})
has been reduced (from 10) to 5, the minimum number to allow
further splitting (\texttt{minsplit}) has been reduced (from 20)
to 10, and the complexity parameter has been reduced to
\texttt{cp = 0.001} .\label{fig:Car90-loosen}}
\end{figure}

\begin{figure}[h]
\begin{minipage}[c]{0.65\linewidth}
<<car90x-plotcp, eval=showFigs, fig.width=4.0, fig.height=3.25, out.width="0.92\\textwidth">>=
fig8.4()
@ %
\end{minipage}
\hspace*{0.05\linewidth}
\begin{minipage}[c]{0.28\linewidth}
\caption{Change in cross-validated error rate, relative
to baseline error, with successive splits.  Because of the
random element that arises from the cross-validation,
the tree that is fitted and the pattern of change of
cross-validated error will commonly change from one
run to the next.\label{fig:ofitx}}
\end{minipage}
\end{figure}

\begin{figure}[h]
\begin{minipage}[c]{0.65\linewidth}
<<Car90-rf, eval=showFigs, echo=FALSE, fig.width=4.0, fig.height=3.5>>=
fig8.5()
@ %
\end{minipage}
\hspace*{0.05\linewidth}
\begin{minipage}[c]{0.28\linewidth}
\caption{Plot of \texttt{Mileage} versus \texttt{tonsWt}, with
fitted values from a \texttt{randomForest} regression shown
as horizontal bars.\label{fig:car90-rf}
}
\end{minipage}
\end{figure}

<<errsmat, eval=showFigs>>=
if(exists('errsmat'))nrep <- nrow(errsmat) else {
  nrep <- 25
  errsmat <- cfRF(nrep=nrep)
}
@ %

\begin{figure}
<<cf-models-rfgam, echo=FALSE, eval=showFigs, fig.width=7, fig.height=7, out.width="0.8\\textwidth">>=
meuse <- getmeuse()
fig8.6()
@ %
\caption{Scatterplot matrix of accuracies
for the several models.
         Each panel shows the line $y=x$.
         The label rfOOB is out-of-bag accuracy
         for the 95 training set observations,
         while rfTest is test data accuracy,
         for a random forest model.
         Other results are test set accuracy from
         fitting a random forest model to residuals
         from a preliminary smooth. Labels are
         the name of the formula for the smooth.
         Random forest fits were from \Sexpr{if(exists('nrep'))nrep else NA}
         bootstrap samples.
         \label{fig:gamRF}
}
\end{figure}

\end{document}
