---
title: "Tree based models"
author: ''
output:
  html_document:
    toc: yes
  pdf_document: null
  pdf_presentation: null
  slidy_presentation:
    css: slidy.css
    widescreen: yes
  word_document: null
---



# Introduction

## The regression (and later, classification) tree model

What should we do when we have no idea of the function that generates the data $y=f(x)$?

We can use polynomial regression (Monday's lecture)

Or we can use another idea which is based on the following 'decision' model.

Here is the magic formula. We will construct an estimate of God's function in the following manner:

$$ \hat{y} = \sum_i B_i(x) v_i $$
where each parameter $B_i$ is a boolean expression and only one boolean is true for any given value $x$ of the predictor variables.

We will construct the boolean value by combining <b>intervals for each predictor</b> of the form $x < c$ or $x > c$.

For example, if I have two predictors $x,z$, I may construct the boolean expression

(x < 1.5) and (z > 3.)

this is a rectangle in the space $(x,z)$ of predictors.

<img src="https://mathstat.dal.ca/~fullsack/stat2450/images/rtree.png">


# The sweet spot

As explained in the last lecture, we have to find the right model complexity

For regression trees, we can decimate (prune) trees to decrease complexity

We will use the idea of regularization to find the sweetest complexity

It will be given by an optimal value of a hyperparameter (penalizing complexity)

We will use cross-validation for finding the sweet spot



# Fitting a Regression Tree

Figure 8.1 of "Introduction to Statistical Learning" shows the initial part of
a regression tree which is used to model the logarithm of a baseball players
salary, using number of years played, and number of hits in the previous year
as predictors.

The data *Hitters* is part of the ISLR library.  The function **tree** used to
fit a regression tree is part of the library **tree**.  After loading the
data, records with missing salaries are removed, and only the variables
Salary, Years and Hits are retained.  Salary is then transformed to log(Salary).

```{r}
library(ISLR)  #contains the Hitters data set
library(tree)  #contains functions for tree construction
baseball=Hitters[!is.na(Hitters$Salary),] #remove records with missing salary
baseb=baseball[,c("Salary","Years","Hits")] #keep only 3 variables
baseb$Salary=log(baseb$Salary) #model log(salary)
```

The function tree is used to construct the regression tree, which is then
printed, plotted, and text is added to the plot.

```{r}
basebtree=tree(baseb) #construct the tree
basebtree #print the tree
plot(basebtree) #plot the tree
text(basebtree) #label the tree
```

Following are a number of functions which can be used to construct bits of the
tree.

```{r}
n=dim(baseb)[1] #number of cases
M=dim(baseb)[2] #variables

SSE=function(data,pred,cut){
#split data into two sets, data1 and data2
#returns sums of squares and means for each of data1,data2
 data1=data[pred<=cut]
 data2=data[pred>cut]
#calculate the error sum of squares for each set
 M1=mean(data1)
 M2=mean(data2)
 SSE1=sum((data1-M1)^2)
 SSE2=sum((data2-M2)^2)
#calculate the total error sum of squares
 return(c(SSE1,SSE2,M1,M2))
 }

gencuts=function(x){
 #calculates all possible cutpoints for x
xuo=sort(unique(x))  #sorted unique values of x
cuts=rep(0,(length(xuo)-1)) 
nc=length(cuts)
for (i in 1:nc)cuts[i]=mean(xuo[i:(i+1)]) #midpoints 
return(cuts)}

cut1=function(y,x){
 cuts=gencuts(x) #get cutpoints for variable x
 ncuts=length(cuts) #number of cut points
 SSEs=matrix(rep(0,5*ncuts),byrow=T,ncol=5)
 for (i in 1:length(cuts)){
   SSEs[i,2:5]=SSE(y,x,cuts[i])
   SSEs[i,1]=sum(SSEs[i,2:3])}
   minp=(1:length(cuts))[SSEs[,1]==min(SSEs[,1])]
#returns best cut point, sums of squares for 2 branches, means for 2 branches
 return(c(cuts[minp],SSEs[minp,]))      
}
```

Here's an example of what **gencuts** does.

```{r}
x <- baseb$Hits
xuo=sort(unique(x))  #sorted unique values of x
plot(xuo)

```

```{r}
table(baseb$Hits)
gencuts(baseb$Hits)
```

   + Try all possible cutpoints $C$, and partition the dataset according to whether $Hits<C$ or $Hits>C$.

   + For each subset,  calculate the mean and sum of squares of the dependent variable $Y$.  Here $Y$
is the $log(Salary)$.  

   + Calculate the  total sum of squares $SSE = SSE_1 + SSE_2$

where  SSE1 and SSE2 are the error sums of squares for the first and second subsets.
   
   + pick the cutopoint $C$ which minimizes $SSE$.  This
   gives us the minimum sum of squares when we partition
   the dataset using Hits.  
   
   + Do the same thing, but  using the variable Years
   
   + Choose the $best$ predictor variable at this stage
   as that variable which minimizes SSE.



```{r}
# find the best variable for the first partition
 for (j in 2:3)print(cut1(baseb[,1],baseb[,j]))
```

From the output, we see that the best choice is to 
partition using $Years<4.5$ and $Years>4.5$.

```{r}
  baseb1=baseb[baseb[,2]<4.5,]    
  baseb2=baseb[baseb[,2]>4.5,]
```



   + proceed recursively down each of the associated
branches of the tree.  Consider a further partition
of the branch identified by $Years>4.5$.

```{r}
for (j in 2:3)print(cut1(baseb2[,1],baseb2[,j]))  
```

  + the next split on this branch uses $Hits<117.5$ vs $Hits>117.5$.

  + continue recursively down all branches, until some
stopping criterion is satisfied.


## Example: use all predictors to grow the tree, work with Salary on original scale.

```{r}
baseball=Hitters[!is.na(Hitters$Salary),] 
basebtree=tree(Salary~.,data=baseball) #construct the tree
basebtree #print the tree
plot(basebtree) #plot the tree
text(basebtree) #label the tree
```





### pruning

Remember polynomial regression? Here the number of leaves is the complexity measure of the model!

Let us vary the complexity of our model!

```{r}
mysize = 2
bp2=prune.tree(basebtree,best=mysize)
plot(bp2)
text(bp2)
```

```{r}
mysize = 5
bp5=prune.tree(basebtree,best=mysize)
plot(bp5)
text(bp5)
```



```{r}
mysize = 8
bp8=prune.tree(basebtree,best=mysize)
plot(bp8)
text(bp8)
```


### Cross-validated choice of tree size

   + Algorithm 8.1 in ISLR

   + For each $\alpha$, evaluate the penalized sum of squares
   
$$PSSE = \sum_{m=1}^{|T|} \sum_{i: x_i \in R_m} (y_i - \hat y_{R_m})^2 + \alpha |T|$$

   + where $|T|$ is the number of terminal nodes in the tree

   + $R_m$ is the subset of the predictor space corresponding to the $m$'th
   terminal node.

   + for each $\alpha$ there is a tree $T$ such that PSSE is as small as possible.
   
   + use $k-fold$ cross validation to choose the value of $\alpha$ which
   minimizes PSSE.  Equivalently, choose the size of tree which minimizes PSSE. The default value of $k$ is 10.



<img src="https://mathstat.dal.ca/~fullsack/stat2450/images/algo81.png">


### Carry out the cross validation, and plot the penalized sum of squares vs tree size.

```{r}
basebcv=cv.tree(basebtree)
plot(basebcv$size,basebcv$dev,type='b')
bestsize=basebcv$size[basebcv$dev==min(basebcv$dev)]
```

   + The cross validated choice of tree size is a tree with `r bestsize` terminal nodes.

   +  Use **prune.tree** to find the best tree of this size.

```{r}
baseball.prune=prune.tree(basebtree,best=bestsize)
plot(baseball.prune)
text(baseball.prune)
```

   + Note:  there is sampling variability in the procedure due to the random partition into $k$ folds.  Each time you run the procedure, you are likely to get a different "best" tree.  Hopefully they are not too different.
   
   +** Random forests** and/or **bagging** tend to give gives good predictions, with reduced variability. 


## Fitting a random forest

### Example: fitting random forest using 3 predictors at each split.

```{r}
library(randomForest)
n=nrow(baseball)
attach(baseball)
n2=ifelse(n%%2==0,n/2,(n+1)/2)
index=sample(1:n,n2,replace=F)
bbtrain=baseball[index,]
bbtest=baseball[-index,]
saltrain=Salary[index]
saltest=Salary[-index]
salary.rf=randomForest(Salary~.,data=bbtrain,mtry=3,importance=T)
importance(salary.rf)
varImpPlot(salary.rf)
shattrain=predict(salary.rf,newdata=bbtrain)
shattest=predict(salary.rf,newdata=bbtest)
paste("training MSE = ", mean((saltrain-shattrain)^2))
paste("test MSE = ", mean((saltest-shattest)^2))
```

### Example: bagging, fits a random forest with mtry=19.

```{r}
salary.rf=randomForest(Salary~.,data=bbtrain,mtry=19,importance=T)
importance(salary.rf)
varImpPlot(salary.rf)
shattrain=predict(salary.rf,newdata=bbtrain)
shattest=predict(salary.rf,newdata=bbtest)
paste("training MSE = ", mean((saltrain-shattrain)^2))
paste("test MSE = ", mean((saltest-shattest)^2))
```

## Note on cross-validation






The general procedure is as follows:

<ol>
<li>Shuffle the dataset randomly.
<li>Split the dataset into k groups
<li>For each unique group:
<ol>
<li>Take the group as a hold out or test data set
<li>Take the remaining groups as a training data set
<li>Fit a model on the training set and evaluate it on the test set
<li>Retain the evaluation score and discard the model
</ol>
<li>Summarize the skill of the model using the sample of model evaluation scores
</ol>

Importantly, each observation in the data sample is assigned to an individual group and stays in that group for the duration of the procedure.

This means that each sample is given the opportunity to be used in the hold out set 1 time and used to train the model k-1 times.



```{r}
# example 3-fold cv
d <- c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6) # data

# partitiom data into 3 folds
d1 <- c(0.2, 0.5)
d2 <- c(0.1, 0.3)
d3 <- c(0.4, 0.6)

# we have 3 models
# model 1 is fitted on d-d1 = d2+d3 and we use d1 for predictions
# etc





```





<ul>
<li><a href=""></a>
<li><a href="https://machinelearningmastery.com/k-fold-cross-validation/">k-fold-cross-validation</a>
<li><a href="https://link.springer.com/content/pdf/10.1007%2F978-1-4899-7502-7_28-1.pdf">Bias-variance tradeoff</a>
</ul>

