---
title: "Bagging, Randomforest, boosting"
author: ''
output:
  html_document:
    toc: yes
  pdf_document: null
  pdf_presentation: null
  slidy_presentation:
    css: slidy.css
    widescreen: yes
  word_document: null
---

# Objectives


 + CART (single tree)
 
 + Bagging
 
 + Random Forest
 
 + Boosting


 Main idea: Averaging. 
 
 We move to ensemble methods (pro: better prediction, con: worse interpretability)
 
 
# Principle of methods and DATASET 1

## Variability in tree estimates

   + There is considerability in the estimated tree and the associated predictions, even when cross validaion
is used to choose the size of the tree.
   + The following example illustrates this variability, using 5 random splits into training and test sets.

```{r}
heart=read.csv("http://faculty.marshall.usc.edu/gareth-james/ISL/Heart.csv")
#heart=read.csv("http://www-bcf.usc.edu/~gareth/ISL/Heart.csv")
library(tree)
nmiss=apply(is.na(heart),1,sum) #number of missing values, by row
heart=heart[nmiss==0,]  #remove cases with missing values
attach(heart)
nrec=nrow(heart);nrec
ntrain=ifelse(nrec%%2==0,nrec/2,(nrec+1)/2);ntrain

```
 
 
```{r echo=F}
Ntries=5
for (i in 1:Ntries){
  print(paste("random split number ",i," into training and test sets"))
train=sample(1:nrec,ntrain,replace=F)
AHD.test=AHD[-train]
AHD.train=AHD[train]
heart.train=heart[train,] #training set
heart.test=heart[-train,] #test set
heart.tree=tree(AHD~., data=heart.train)
heartcv=cv.tree(heart.tree, FUN=prune.misclass)
bestsize=heartcv$size[heartcv$dev==min(heartcv$dev)]
bestsize=ifelse(length(bestsize)==1,bestsize,bestsize[1])
heart.prune=prune.misclass(heart.tree,best=bestsize)
plot(heart.prune)
text(heart.prune)
pred.train=predict(heart.prune,data.frame(heart.train),type="class")
ctable=table(pred.train,AHD.train)
print("Confusion matrix for training set")
print(ctable)
trainrate= (sum(ctable)-sum(diag(ctable)))/sum(ctable)
pred.test=predict(heart.prune,data.frame(heart.test),type="class")
print("Confusion matrix for test set")
ctable=table(pred.test,AHD.test)
print(ctable) 
testrate= (sum(ctable)-sum(diag(ctable)))/sum(ctable)
print(paste("Training error rate = ",round(100*trainrate,2)," %"))
print(paste("Test error rate = ",round(100*testrate,2)," %"))
}
      
```

## Bagging

   + **Bagging**, or **bootstrap aggregation** is a procedure which was
developed to deal with the high degree of variability in the
tree prediction.

### Bagging for regression trees

   + Suppose that $\hat f^1(x)$, $\hat f^2(x), \ldots, \hat f^B(x)$, are
   predictions generated using $B$ independent training sets.
   + Average the classifications to form

$$\hat f_{avg}(x) = \frac{1}{B} \sum_{b=1}^B \hat f^b(x)$$


   + But how do we get independent training sets?
   + We bootstrap, by taking repeated samples from our **SINGLE** original
   training set.
   + This gives $B$ bootstrap training sets.

   + Train the method on the $b$'th bootstrapped training set to prediction
   $f^{*b}(x)$, for $b=1,2, \ldots , B$.

   + Average to get

$$\hat f_{bag}(x) = \frac{1}{B} \sum_{b=1}^B \hat f^{*b} (x)$$

   + This is **the bagged estimate* of $f(x)$.

   + The procedure is called **bagging**

### Bagging for classification trees

   + How do we average predictions for classification trees?
   
   + We use a **majority vote**.  For a given $x$, predict the outcome using the
   prediction generated by the **majority** of the $B$ predictions.
   
   + For example, for a fixed $x$, and predicting a binary outcome $AHD$ with $B=100$, 
   suppose that 73 of the trees predict $AHD=Y$, and the other 27 trees predict $AHD=N$.
   Then the majority prediction is $Y$.
   

```{r}
library(randomForest) #random forest library
```

### Out-of-bag (OOB) error rate

   + For each bootstrap sample, about 1/3 of the sample is left out. Why?
   + These left out observations, are referred to as *out of bag* (or OOB) observations.
   + A tree is estimated for each bootstrap sample.
   + For each data point, say case $i$,  not used in the tree construction,  the outcome is predicted.  A majority vote is taken using all trees which did not use the $i$'th observation in their construction, and this gives the overall classification for the outcome, say $Y_i$.  
   + This is repeated for each of the $n$ observations, leading to an overall mislassification rate, the **out of bag error rate**.  It is a valid measure
   of **test error** for the bagged model.
   + The OOB misclassification rate will typically be close, but not indentical, to the misclassificaton rate on an independent test set.

```{r}
train=sample(1:nrec,ntrain,replace=F)
AHD.train=AHD[train]
AHD.test=AHD[-train]
heart.bag=randomForest(AHD~.,data=heart,subset=train,mtry=14,importance=T)
heart.bag
pred.test=predict(heart.bag,newdata=heart[-train,])
print("Confusion matrix for test set")
table(pred.test,AHD.test)
```
### Variable importance measures

* for each observation in training or test sets, bagging will give a prediction, and the predictions tend to be better than for individual trees
* one downside of the bagging procedure is that it doesn't give a single tree, so has less interprability than the single tree estimate.  **Bagging improves predictability at the expense of interpretation**
* We can record the total amount that the Gini index is decreased by splits over a particular
   predictor, averaged overall all $B$ trees. This gives an overall estimate of importance of
   each predictor variable.
*  We can do the same using prediction accuracy, measuring the average 
   decrease in prediction error resulting from splits over each predictor.
   
```{r}   
importance(heart.bag)
varImpPlot(heart.bag)
```

### Random forests

* A **random** forest is like a bagged estimate, except that the
predictors chosen at each split are a random sample of the complete
set of predictors. 

* The number of variables $m$ used at each split is a tuning parameter.
   The default is $\sqrt p$ (rounded), where $p$ is the total number of predictors.
   
* Bagging is a particular case of the random forest algorithm, with 
   $m=p$.
   
* Rationale:  
    + suppose there is one very strong predictor, and several moderaly strong predictors.  Then each tree in the bagged collection will typically have 
   the same variable at the first split. 
    + so all bagged trees will look fairly similar
    + so predictions for the different trees will be quite similar
    + so variance of the overall prediction will not be dramatically reduced. Why?
* by randomly choosing variables at each split, the correlation among trees is reduced.
* With a large number of correlated predictors, a small value of $m$ typically works well.   


```{r}
heart.rf=randomForest(AHD~.,data=heart, subset=train,importance=T)
heart.rf
importance(heart.rf)
varImpPlot(heart.rf)
pred.test=predict(heart.rf,newdata=heart[-train,])
print("Confusion matrix for test set")
table(pred.test,AHD.test)
```


### Boosting (you won't be asked any questions on boosting)

* Fit a tree, and obtain residuals
* Fit a second tree on the residuals from the first tree
* Update the residuals by combining predictions from the second tree with
predictions from the first tree.
* Iterate


```{r}
library(gbm)          #boosting library
Y=ifelse(AHD.train=="Yes",1,0)  #need a 0/1 outcome for boosting
heart.boost=gbm(Y~.-AHD,data=heart[train,],distribution="bernoulli",n.trees=5000)
summary(heart.boost)
pred.train=predict(heart.boost,newdata=heart[train,],n.trees=5000,type="response")
pred.train=ifelse(pred.train<=.5,0,pred.train)
pred.train=ifelse(pred.train>.5,1,pred.train)
table(pred.train,AHD.train)
pred.test=predict(heart.boost,newdata=heart[-train,],n.trees=5000)
 pred.test=ifelse(pred.test<=.5,0,pred.test)
pred.test=ifelse(pred.test>.5,1,pred.test)
table(pred.test,AHD.test)
```



# DATASET 2: Boston : regression of median value response


## read the data

```{r}
library(MASS)
boston <- Boston #I generally prefer lower case
#setDT(boston)
str(boston)
```



Let's split the data and build one regression tree

```{r}
#library(magrittr)
# for sample_frac
library(dplyr)
library(tree)
set.seed(1)
# infix notation (pipe)
#boston_train = Boston %>% sample_frac(.5)
#boston_test = Boston %>%
#  setdiff(boston_train)
# magrittr for this syntax
# https://stackoverflow.com/questions/24536154/what-does-mean-in-r
# iris %>% head() is equivalent to head(iris)

# I (pf) prefer this syntax:
boston_train = sample_frac(Boston,.5)
boston_test = setdiff(Boston, boston_train)

tree_boston=tree(medv~., boston_train)

summary(tree_boston)
```

## Let's bag it!

Let's bag it!

Use $m=p=13$

mtry = 13 means that we use all 13 predictors in each split of the tree.

This implies that we are using RF for BAGGING ($S x B$)

```{r}
library(randomForest)
set.seed(1)
bag_boston = randomForest(medv~., 
                          data = boston_train, 
                          mtry = 13, 
                          importance = TRUE)
bag_boston
```



Now evaluate the performance of the BAGGED model on the test set:

```{r}
library(ggplot2)
bagged_estimate = predict(bag_boston, 
                          newdata = boston_test)

ggplot() + 
    geom_point(aes(x = boston_test$medv, y = bagged_estimate)) +
    geom_abline()
```

## MSE of the bagged model

Now estimate the 'MSE' of the BAGGED model

```{r}
bag_boston_25_trees = randomForest(medv~., data =  boston_train, mtry = 13, ntree = 25)
bagged_estimate_25_trees = predict(bag_boston_25_trees, newdata = boston_test)
mse_baggedmodel <- mean((bagged_estimate_25_trees - boston_test$medv)^2)
mse_baggedmodel
```

What if we take less trees?



```{r}
mse <- function(n,test,train){
bagm = randomForest(medv~., data =  train, mtry = 13, ntree = n)
bagpred = predict(bagm, newdata = test)
return(mean((bagpred - test$medv)^2))
}

mse(5,boston_test,boston_train)
```


What if we take more trees?



```{r}
mse(100,boston_test,boston_train)
```


## Was bagging useful?

Now, was bagging useful?

To test this, we are going to compare (compare what? your answer here) the BAGGED model to a SINGLE tree (CART model) that is optimized by cross-validation.


Let us build the CONCURRENT MODEL (CART-CV)

It will be a cross-validated CART model


Let us run thi

```{r}
tree_boston=tree(medv~., boston_train)
cv_boston = cv.tree(tree_boston)
plot(cv_boston$size, cv_boston$dev, type = 'b')
```

hmm, looks like cv is selecting 6 or 7 as best size (minimum deviance). But 7 more often than 6.
OK I take 7 as optimal tree size.


And I prune the full tree to this size (russian dolls):

```{r}
prune_boston = prune.tree(tree_boston, 
                          best = 7)
plot(prune_boston)
text(prune_boston, pretty = 0)
```


Now that I have built my alternative model, I can compare it to the bagged model.


But how should I compare it?


Simply by calculating the MSE of the CART model on the test set.


```{r}
single_tree_estimate = predict(prune_boston, 
                               newdata = boston_train)
ggplot() + 
    geom_point(aes(x = boston_test$medv, y = single_tree_estimate)) +
    geom_abline()

mse_cartmodel <- mean((single_tree_estimate - boston_test$medv)^2)
mse_cartmodel
```



So now we compare the two MSEs. Which one is better, what is your conclusion?


```{r}
c(mse_cartmodel,mse_baggedmodel)
```

Yes BAGGING WAS USEFUL!


# OK but can we do better ? (1)


How?


RF!

By default, when no value is passed to mtry, the function randomForest uses m=p/3 variables (features)
for regression, and sqrt(m) features for classification.

Let us try $m=6$


```{r}
mserf <- function(m,test,train){
rfb = randomForest(medv~., 
                         data = train, 
                         mtry = m, 
                         importance = TRUE)
rfest = predict(rfb,newdata = test)
return(mean((rfest - test$medv)^2))
}


m=6
mb=13
h=NULL
hb = NULL
for(i in 1:50){
h=c(h,mserf(m,boston_test,boston_train))
hb = c(hb,mserf(mb,boston_test,boston_train))
}
hist(h)
hist(hb)


```


Did RF improve over bagging?


```{r}
c(median(h),median(hb))
```


what do you think?


# OK but can we do better ? (2)

How?


Boosting!


I leave this as an exercise.

# DATASET 3: Carseats sales

Can you replicate the same methology for Carseats sales?

Response : either choose the original Sales response (continuous: regression exercise), or transform
the response into a two levels factor (low, high: classification exercise)


<a href="https://daviddalpiaz.github.io/r4sl/trees.html">(src:daviddalpiaz.github.io/r4sl/trees.html)</a> 

## Full tree


```{r eval=F}
library(tree)
library(ISLR)

data(Carseats)
#?Carseats
str(Carseats)

Carseats$Sales = as.factor(ifelse(Carseats$Sales <= 8, "Low", "High"))
str(Carseats)

# Fit an unpruned classification tree using all of the predictors.

seat_tree = tree(Sales ~ ., data = Carseats)
# seat_tree = tree(Sales ~ ., data = Carseats, 
#                  control = tree.control(nobs = nrow(Carseats), minsize = 10))
summary(seat_tree)
plot(seat_tree)
text(seat_tree, pretty = 0)
title(main = "Sales of carseats: full classification tree: no pruning!")
```



```{r eval=F}

#  details of the splits.
seat_tree

```


## Variability in tree estimates

   + There is considerability in the estimated tree and the associated predictions, even when cross validaion
is used to choose the size of the tree.
   + The following example illustrates this variability, using 5 random splits into training and test sets.

```{r eval=F}
# Carseats
attach(Carseats)
nrec=nrow(heart);nrec
ntrain=ifelse(nrec%%2==0,nrec/2,(nrec+1)/2);ntrain

```