---
title: 'Statistical Thinking using Randomisation and Simulation'
subtitle:  'Ensemble models using bootstrap'
author: Di Cook (dicook@monash.edu, @visnut)
date: "W8.C1"
output:
  xaringan::moon_reader:
    css: ["default", "myremark.css"]
    self_contained: false
    nature:
      highlightStyle: github
      highlightLines: true
      countIncrementalSlides: false
---

```{r setup, include = FALSE}
library(knitr)
opts_chunk$set(
  message = FALSE,
  warning = FALSE,
  cache = FALSE,
  echo=FALSE,
  fig.height = 6,
  fig.width = 9,
  fig.align='center'
)
options(digits=3)
library(tidyverse)
library(gridExtra)
library(broom)
```

# Bootstrapping

- Taking bootstrap samples, then fitting a model to each sample, and combining model predictions improves predictive accuracy (reduces error)
- The famous example is a random forest. 
- A forest is constructed by taking a bootstrap sample, and a random sample of available explanatory variables, fitting a tree to each, then pooling the model predictions.

---
# Example - same as for trees

```{r}
set.seed(8000)
x=sort(runif(100)-0.5)
df <- data.frame(x, y=c(x[1:50]^2, x[51:75]*2, -x[76:100]^2)+rnorm(100)*0.1)
ggplot(df, aes(x=x, y=y)) + geom_point()
```

---
# Bootstrap samples

```{r fig.width=12, fig.height=5}
df1 <- df %>% sample_n(size=100, replace=TRUE)
df2 <- df %>% sample_n(size=100, replace=TRUE)
p1 <- ggplot(df1, aes(x=x, y=y)) + geom_point(alpha=0.3, size=3) +
  xlim(range(df$x)) + ylim(range(df$y)) + ggtitle("Sample 1")
p2 <- ggplot(df2, aes(x=x, y=y)) + geom_point(alpha=0.3, size=3) +
  xlim(range(df$x)) + ylim(range(df$y)) + ggtitle("Sample 2")
grid.arrange(p1, p2, ncol=2)
```

---

.pull-left[
Sample 1
```{r fig.width=5}
library(rpart)
library(rpart.plot)
df1_rp <- rpart(y~x, data=df1)
rpart.plot(df1_rp)
```
]

.pull-right[
Sample 2
```{r fig.width=5}
df2_rp <- rpart(y~x, data=df2)
rpart.plot(df2_rp)
```
]

---
# Combine predictions

```{r echo=TRUE}
df_rp_aug <- cbind(df, p1=predict(df1_rp, df), p2=predict(df2_rp, df))
df_rp_aug[91:100,]
```


---

.pull-left[
```{r echo=TRUE}
df_rp_aug <- df_rp_aug %>% mutate(p=(p1+p2)/2)
df_rp_aug[91:100,]
```
]

.pull-right[
```{r fig.width=4, fig.height=3}
ggplot(df_rp_aug, aes(x=x, y=y)) + geom_point() +
  geom_point(aes(y=p), colour="hotpink", alpha=0.7, size=3)
```
]
---
# Compare errors

```{r}
MSE <- function(obs, pred) {
  e <- obs - pred
  mse <- sum(e^2)/length(e)
  return(mse)
}
options(digits=5)
```

.pull-left[
Tree model on full data
```{r echo=TRUE}
df_rp <- rpart(y~x, df)
df_rp_aug <- cbind(df_rp_aug, full_p=predict(df_rp, df))
MSE(df_rp_aug$y, df_rp_aug$full_p)
```
]

.pull-right[
Average of two models
```{r echo=TRUE}
MSE(df_rp_aug$y, df_rp_aug$p)
```
]

Reduces the error (a little)!

---
# Add more trees

.pull-left[
MSE from averaging 5 trees

```{r}
df3 <- df %>% sample_n(size=100, replace=TRUE)
df3_rp <- rpart(y~x, df3)
df_rp_aug <- cbind(df_rp_aug, p3=predict(df3_rp, df))
df4 <- df %>% sample_n(size=100, replace=TRUE)
df4_rp <- rpart(y~x, df4)
df_rp_aug <- cbind(df_rp_aug, p4=predict(df4_rp, df))
df5 <- df %>% sample_n(size=100, replace=TRUE)
df5_rp <- rpart(y~x, df5)
df_rp_aug <- cbind(df_rp_aug, p5=predict(df5_rp, df))
df_rp_aug <- df_rp_aug %>% mutate(p=(p1+p2+p3+p4+p5)/5)
MSE(df_rp_aug$y, df_rp_aug$p)
```

reduces error by quite a bit more.
]

.pull-right[
```{r fig.width=4, fig.height=3}
ggplot(df_rp_aug, aes(x=x, y=y)) + geom_point() +
  geom_point(aes(y=p), colour="hotpink", alpha=0.7, size=3)
```
]

---
# Random forests

- Multiple trees, fit to bootstrap samples
- Sample cases, using bootstrapping (ones not chosen are called out-of-bag, used for testing purposes)
- Sample variables to choose from at each node of each tree
- Lots of diagnostics generated!
- Bagging stands for `bootstrap aggregation`. Combine the results from multiple models built on different bootstrap samples. Random forests are an example of bagging
- Bagging can be used with almost any classifier
- Bagging reduces variation in estimates

---
# Example: Tennis performance

Women's grand slam performance statistics for 2012. 

Response variable: `won` whether the player won the match or not

Explanatory variables: `Winners`, `Unforced.Errors`, `First.Serves.In`, `Total.Aces`, `Total.Double.Faults`, `First.Serve.Pts.Won`, `Receiving.Points.Won`, `Break.Point.Conversions`

Split the data into two pieces, training and test. Model will be fit on training set, error will be calculated on test set.

```{r echo=TRUE}
library(caret)
tennis <- read_csv("../data/tennis_2012_wta_fixed.csv")
indx <- createDataPartition(tennis$won, p=0.67, list=FALSE)
train <- tennis[indx,]
test <- tennis[-indx,]
#train %>% count(won)
#test %>% count(won)
```

---
# Random forest fit to training

```{r}
library(randomForest)
tennis_rf <- randomForest(won~Winners+Unforced.Errors+
                            RPW.TPW+W.minus.Aces+
                  SPW.TPW+First.Serves.In+Total.Aces+
                  Total.Double.Faults+First.Serve.Pts.Won+
                  Receiving.Points.Won+Break.Point.Conversions,
                  data=train, na.action=na.omit, importance=TRUE)
tennis_rf
```

Test error

```{r}
pred <- predict(tennis_rf, test)
MSE(test$won, pred)
```

---
# Diagnostics: variable importance

```{r}
as.data.frame(tennis_rf$importance) %>%
  mutate(variable=rownames(tennis_rf$importance)) %>%
  arrange(desc(`%IncMSE`))
```

With the sampling of variables, it is possible to assess the effect of each variable on the accuracy of predictions. For winning a tennis match `Receiving.Points.Won` is by far the best variable.  Keep this high, and your chances of winning are much improved.

---
# Effect of number of trees

.pull-left[
- The training and test error changes as number of trees averaged increases.
- Training error could go to zero. It is buffered here, because only the "out-of-bag" cases are predicted.
]

.pull-right[
```{r fig.width=6, fig.height=4}
ntree <- seq(1, 50, 1)
tr_err <- rep(100, length(ntree))
ts_err <- rep(100, length(ntree))
for (i in 1:length(ntree)) {
  tennis_rf <- randomForest(won~Winners+Unforced.Errors+
                            RPW.TPW+W.minus.Aces+
                  SPW.TPW+First.Serves.In+Total.Aces+
                  Total.Double.Faults+First.Serve.Pts.Won+
                  Receiving.Points.Won+Break.Point.Conversions,
                  data=train, na.action=na.omit, importance=TRUE,
                  ntree=ntree[i])
  tr_err[i] <- tennis_rf$mse[i]
  pred <- predict(tennis_rf, test)
  ts_err[i] <- MSE(test$won, pred)
}
df <- data.frame(ntree, tr_err, ts_err) %>% gather(type, error, -ntree)
ggplot(df, aes(x=ntree, y=error, colour=type)) + 
  geom_point() + geom_smooth(se=FALSE) +
  ylim(c(0, 0.35)) + scale_colour_brewer(palette="Dark2")
```
]

- Test error will decline, up to a point, and then plateau, or worse get higher.
- Where the curves cross is the best choice of number of trees, about 20 here.
- Balance complexity vs accuracy (bias vs variance trade-off)

---
# Resources

- [Tutorial on Tree Based Modeling ](https://www.analyticsvidhya.com/blog/2016/04/complete-tutorial-tree-based-modeling-scratch-in-python/)

---
class: inverse middle 
# Share and share alike

<a rel="license" href="http://creativecommons.org/licenses/by/4.0/"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by/4.0/88x31.png" /></a><br />This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">Creative Commons Attribution 4.0 International License</a>.
