index.Rmd

---
title: "Chapter 8: <br>Tree-based Models"
author: "Jeremy Selva"
date: ""
output:
  xaringan::moon_reader:
    chakra: libs/remark-latest.min.js #Enable "offline" use, xaringan::summon_remark()
    css:
      - default
      - libs/remark.css
    lib_dir: libs
    nature:
      ratio: "16:9"
      highlightStyle: googlecode #Code colour style
      highlightLines: true
      highlightLanguage: ["r"]
      countIncrementalSlides: false
      slideNumberFormat: |
        <div class="progress-bar-container">
          <div class="progress-bar" style="width: calc(%current% / %total% * 100%);">
          </div>
        </div>
---

```{r initial, echo = FALSE, cache = FALSE, results = 'hide'}
library(knitr)
library(glue)
library(metathis)
library(webshot2)
library(xaringanExtra)

options(htmltools.dir.version = FALSE, tibble.width = 60, tibble.print_min = 6)
opts_chunk$set(
  echo = FALSE, 
  warning = FALSE, 
  message = FALSE, 
  comment = "#>",
  #fig.path = 'figure/', 
  cache.path = 'cache/',
  dpi = 300, 
  fig.align = 'center', 
  #fig.width = 6,
  #ig.asp = 0.618,
  #fig.show = 'hold',
  cache = TRUE
  #dev = 'svglite'
)
```

```{r screenshot function, include=FALSE}
#' Screenshot Your Title Slide for Share Image
#'
#' Takes a screenshot of your title slide for sharing on Twitter
#' (and other social media sites).
#'
#' @param slides_rmd Your slides file
#' @param path Path to new share image
screenshot_share_image <- function(
  slides_rmd,
  path_image = "share-card.png"
) {
  if (!requireNamespace("webshot2", quietly = TRUE)) {
    stop(
      "`webshot2` is required: ", 
      'remotes::install_github("rstudio/webshot2")'
    )
  }
  
  webshot2::rmdshot(
    doc = slides_rmd,
    file = path_image,
    vheight = 600,
    vwidth = 600 * 191 / 100,
    rmd_args = list(
      output_options = list(
        nature = list(ratio = "191:100"),
        self_contained = TRUE
      )
    )
  )
  
  path_image
}
#screenshot_share_image("index.Rmd")
```

```{r metathis, echo=FALSE}
meta() %>%
  meta_name("github-repo" = "JauntyJJS/islr2-bookclub-cohort3-chapter8") %>% 
  meta_social(
    title = "Chapter 8: Tree-based Models",
    description = glue::glue("
    Presented at R for Data Science book club cohort 3 on \\
    An Introduction to Statistical Learning: With \\
    Applications in R Second Edition.
    "),
    url = "https://jauntyjjs.github.io/islr2-bookclub-cohort3-chapter8",
    image = "https://raw.githubusercontent.com/JauntyJJS/islr2-bookclub-cohort3-chapter8/main/share-card.png",
    image_alt = glue::glue("
      Title slide of Chapter 8: Tree-based Models \\
      of An Introduction to Statistical Learning: With \\
      Applications in R Second Edition for R for Data Science book club cohort 3.
      "),
    og_type = "website",
    og_author = "Jeremy Selva",
    twitter_card_type = "summary_large_image",
    twitter_creator = "@JauntyJJS"
  )
```

```{r xaringanExtra-share-again, echo=FALSE}
# xaringanExtra::embed_xaringan(
#  url = "https://jauntyjjs.github.io/islr2-bookclub-cohort3-chapter8",
#  ratio = "16:9"
# )
xaringanExtra::use_share_again()
```

# Introduction

This chapter introduces tree-based methods for regression and classification

- Decision Trees
- Bagging
- Random Forest
- Boosting
- Bayesian Additive Regression Trees 

---

# 8.1.1 Regression Trees

1. Divide the predictor space — that is, the set of possible values
for $X_1,X_2, . . . ,X_p$ — into $J$ distinct and .blue[**non-overlapping**] regions,
$R_1,R_2, . . . ,R_J$ .

2. For every observation that falls into the region $R_j$, where $1 \le j\le J$, take the .blue[**mean**] of the response values for the training observations in $R_j$.

```{r echo = FALSE, out.width="60%"}
knitr::include_graphics("img/figure_8_2.jpg")
```

---

# 8.1.1 Regression Trees

The goal is to find boxes $R_1, . . . ,R_J$ that .blue[**minimize**] the $RSS$, given by


$$\mathrm{RSS}=\sum_{j=1}^{J}\sum_{i{\in}R_j}^{}(y_i - \hat{y}_{R_j})^2$$

where $\hat{y}_{R_j}$ is the .blue[**mean**] response for the training observations within the $j$th box.

Unfortunately, it is .red[**computationally infeasible**] to consider every
possible partition of the feature space into $J$ boxes.

---

# 8.1.1 Regression Trees

One way is given all predictors $X_1,X_2, . . . ,X_p$, find the .blue[**best**] cut off point $s$ for a given predictor $X_k$, where $1 \le k\le p$ such that the resulting two regions

.center[
$R_1(k,s)=\left\{X|X_k<s\right\}$ and $R_2(k,s)=\left\{X|X_k{\ge}s\right\}$
]

have the .blue[**least**] $RSS$ for all predictors

$$RSS=\sum_{i:x_i{\in}R_1(k,s)}(y_i-\hat{y}_{R_1})^2 +  \sum_{i:x_i{\in}R_2(k,s)}(y_i-\hat{y}_{R_2})^2$$
where $\hat{y}_{R_1}$ is the .blue[**mean**] response for the training observations in $R_1(k,s)$, and $\hat{y}_{R_2}$ is the .blue[**mean**] response for the training observations in $R_2(k,s)$.

.blue[**Repeat**] the process until no region contains more than five observations.

---

# 8.1.1 Regression Trees

See https://bookdown.org/tpinto_home/Beyond-Additivity/regression-and-classification-trees.html for animation.

.center[![DecisionTree](DecisionTree.gif)]

---

# 8.1.1 Regression Trees

See [StatQuest with Josh Starmer](https://www.youtube.com/watch?v=g9c66TUylZ4) 22 min animation.

.center[
<iframe width="853" height="480" src="https://www.youtube.com/embed/g9c66TUylZ4" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>]

---

# 8.1.1 Regression Trees

However, the previous method may result in a tree that .red[**overfits**] the data.

A better strategy is to 

1. grow a very large tree $T_0$ as before 
2. tune $\alpha$ in the cost complexity pruning formula to prune/penalise $T_0$ to give .blue[**a sequence**] of sub trees $T$s. 
3. apply cross validation or validation set on those sub trees $T$s and pick $\alpha$ that gives the .blue[**best sub tree**] (smallest $RSS$).

The cost complexity pruning formula is 

$$\sum_{m=1}^{|T|}\sum_{x_i{\in}R_m}(y_i-\hat{y}_{R_m})^2 + \alpha|T|$$
where $\alpha\ge0$, $|T|$ is the number of .blue[**terminal nodes**] the sub tree $T$ holds. $R_m$ is the .blue[**box/region**] corresponding to the $m$th terminal node and $y_i-\hat{y}_{R_m}$ is the .blue[**mean**] response for the training observation in $R_m$

---

# 8.1.1 Regression Trees

$$\sum_{m=1}^{|T|}\sum_{x_i{\in}R_m}(y_i-\hat{y}_{R_m})^2 + \alpha|T|$$

When $\alpha = 0$, the sub tree $=T_0$ as the cost complexity pruning formula equals to the $RSS$.

The cost complexity .red[**penalises**] complex trees with a lot of terminal nodes.

Similar to lasso, as $\alpha$ increase, some terminal nodes will not be used.

---

# 8.1.1 Regression Trees

See [StatQuest with Josh Starmer](https://www.youtube.com/watch?v=D0efHEJsfHo) 16 min animation.

.center[
<iframe width="853" height="480" src="https://www.youtube.com/embed/D0efHEJsfHo" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>]

---

# 8.1.2 Classification Trees

In the case of classification, instead of taking the mean of a given region, we take the .blue[**most commonly occurring class**] (or mode) of training observations in the region to which it belongs.

For the measurement of error, instead of $RSS$, one way is to use the fraction of the training observations in that region that do not belong to the most common class.

$$E = 1 - \max_k(\hat{p}_{mk})$$
where $\hat{p}_{mk}$ the .blue[**proportion of training observations**] in the $m$th region that are from the $k$th class. Hence, $0\le\hat{p}_{mk}\le1$

However, this error rate is unsuited for tree-based classification because $E$ does not change much as the tree grows (.red[**lack sensitivity**]).

Instead, .green[**more sensitive**] classification error measurements like the Gini index or the entropy are used.

---

# 8.1.2 Classification Trees

The Gini index and the entropy are very similar measures, where .green[**small values**] indicate that the node contains .blue[**large number of observations from a single class**].

The Gini index is defined by

$$G = \sum_{k=1}^{K}\hat{p}_{mk}(1-\hat{p}_{mk})$$
The value varies between $0$ and $1$,
where, $0$ denotes that .blue[**all**] observations belong to a certain class, $1$ denotes that the observations are .blue[**randomly**] distributed across various classes and
$0.5$ denotes .blue[**equally**] distributed observations into some classes.

The entropy is defined by

$$D = - \sum_{k=1}^{K}\hat{p}_{mk}\log(\hat{p}_{mk})$$

Both Gini index and Entropy are .blue[**small**] when all of the $\hat{p}_{mk}$'s are close to .blue[**zero**] or .blue[**one**].

---

# 8.1.2 Classification Trees

See https://bookdown.org/tpinto_home/Beyond-Additivity/regression-and-classification-trees.html

```{r echo = FALSE}
knitr::include_graphics("img/gini_vs_entropy.jpg")
```

---

# 8.1.3 Trees Versus Linear Models

In terms of prediction accuracy, ... it depends

```{r echo = FALSE, out.width="50%"}
knitr::include_graphics("img/figure_8_7.jpg")
```

---

# 8.1.3 Trees Versus Linear Models

Compared to linear models, decision trees are .green[**easier to explain**] to others because they mirror human decision-making. They can be displayed graphically and easily interpreted by non-experts.

However, it is .red[**less robust**] than the linear models. A small change in the data set can cause a large change to the decision tree. This may hurt the decision trees performance in prediction.

To improve its performance, methods like bagging, random forests, and boosting are used (sometimes at the cost of interpretability). 

---

# 8.2.1 Bagging

Bagging is an approach that tries to “stabilise”(make it more robust to change of data) the predictions from trees

The idea is to generate $B$ different .blue[**bootstrapped training datasets**] to create $B$ different regression trees $\hat{f}^{*b}$ respectively for $1\le b \le B$. Given a new observation to predict called $x$, we use the $B$ regression trees to create $B$ predictions. The final prediction value is calculated by .blue[**averaging**] all the predictions.

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

In the case of classification, we record the class predicted by each of the $B$ trees, and take a .blue[**majority**] vote.

What is missing is how do we estimate the test error rate (without using the test data) to pick the best $B$ value or number of bagged trees ?

---

# 8.2.1 Bagging

Each bagged tree is trained using .blue[**two-thirds**] of the training data. The other one-third of the training data, known as .blue[**out-of-bag**] or $OOB$ observations are then used.

Given an $OOB$ observation, we first predict its response .blue[**using only bagged trees**] that do not use this $OOB$ observation in its training. An $OOB$ prediction is made by averaging the predicted response (regression) or a majority vote (classification) among these specific bagged trees. 

As each $OOB$ observation has an $OOB$ prediction and a true response, the $OOB$ error can be computed via the overall $OOB$ $RSS$ (regression) or classification error like Gini index. The $B$ value that gives the .blue[**least**] $OOB$ error will be taken.

---

# 8.2.1 Bagging

.pull-left[
Bagging usually improves prediction accuracy at the expense of interpretability. 

It is no longer clear which predictors are more important as we now have more than one tree (which may gives differing views on the importance of a given predictor).

An overall summary of the importance of each predictor can be achieved by recording how much the average $RSS$ or Gini index .blue[**improves (or decreases)**] when each tree is splitted over a given predictor.
]

.pull-right[
```{r echo = FALSE, out.width="110%"}
knitr::include_graphics("img/figure_8_9.jpg")
```
]


---

# 8.2.2 Random Forests

One problem with bagging is that each bagged tree may be .red[**highly similar**] to one another. 

This happens if there is a strong predictor in the data set. By the greedy principle, most of the bagged trees will .red[**use this strong predictor**] in the top split, regardless of changes in the training set led by bootstrapping. 

This may cause the tree model to be .red[**less robust**] to changes in the data over time.

To further improve the robustness of the model, a .blue[**change**] is made during the creation of a split in a given bagged tree on bootstrapped training samples. 

Instead of using the full set of $p$ predictors as split candidates, only .blue[**a random sample**] of $m \approx \sqrt{p}$ predictors are used instead.

This reduces the odds of the strong predictor being selected during the split and provide other predictors a chance.

---

# 8.2.2 Random Forests

See [StatQuest with Josh Starmer](https://www.youtube.com/watch?v=J4Wdy0Wc_xQ) 10 min animation.

.center[
<iframe width="853" height="480" src="https://www.youtube.com/embed/J4Wdy0Wc_xQ" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>]

---

# 8.2.2 Random Forests

```{r echo = FALSE, out.width="50%"}
knitr::include_graphics("img/figure_8_8.jpg")
```


---

# 8.2.3 Boosting

Recall in bagging, each tree is built on a bootstrap training data set. 

On the other hand, in boosting, each tree is built on information that the previous trees are unable to catch. That is, we fit a tree using .blue[**the current residuals**], rather than the response.

This is similar to the partial least square components in which each component is build based on the response's residual left by the model created from the previous partial least square components.

The key is to make each tree .blue[**small**] (with a few terminal nodes) but still "different" from one another.

For simplicity, only boosting regression trees are described in the book.

---

# 8.2.3 Boosting

.pull-left[
Boosting has three tuning parameters.

1. The .blue[**number of trees**] $B$

2. The maximum number of .blue[**splits/depth**] $d$ each tree is allow to have. (to make each tree small)

3. A .blue[**shrinkage parameter or learning rate**] $\lambda$ (to make each tree "different" from one another)
]

.pull-right[
```{r echo = FALSE, out.width="120%"}
knitr::include_graphics("img/algorithm_8_2.jpg")
```
]

---

# 8.2.3 Boosting

See [StatQuest with Josh Starmer](https://www.youtube.com/watch?v=3CC4N4z3GJc) 16 min animation.

.center[
<iframe width="853" height="480" src="https://www.youtube.com/embed/3CC4N4z3GJc" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>]

---

# 8.2.3 Boosting

```{r echo = FALSE, out.width="60%"}
knitr::include_graphics("img/figure_8_11.jpg")
```


---

# 8.2.4 Bayesian Additive Regression Trees (BART)

For simplicity, only BART for regression are described in the book.

- In bagging, each tree is built on a .blue[**random set of training data**].
- In random forest, each tree also uses .blue[**a random set of predictors**].
- In boosting, each tree is built to capture .blue[**signal that is not yet accounted**] for by the current set of trees

BART combines all these elements together.

---

# 8.2.4 Bayesian Additive Regression Trees (BART)


The book define some notation

- $K$ be the total .blue[**number of regression trees**] used to create the model (later by adding them together)
- $B$ be the .blue[**number of iterations**] the BART algorithm will run to improve the additive regression trees.
- $\hat{f}^b_k(x)$ be the .blue[**prediction**] of an observation $x$ for the $k$th regression tree used in the $b$th iteration of the BART algorithm, where $1\le k \le K$ and $1\le b \le B$
- $L$, where $L<B$ be the number of .blue[**burn-in**] iterations.

Like boosting, the tree model is created by summing all regression trees at the end of each iteration.

$$\hat{f}^b(x) = \sum_{k=1}^{K}\hat{f}^b_k(x)$$
for $b=1,...,B$

---

# 8.2.4 Bayesian Additive Regression Trees (BART)

In the first iteration of the BART algorithm, all $K$ trees are initialized to
to predict only .blue[**one**] value: $\hat{f}^1_k(x) = \frac{1}{nK}\sum_{i=1}^{n}y_i$
, the .blue[**mean**] of the response values divided by the total number of trees. 

Thus

$$\hat{f}^1(x) = \sum_{k=1}^K\hat{f}^1_k(x) = \sum_{k=1}^K\frac{1}{nK}\sum_{i=1}^{n}y_i = \frac{1}{n}\sum_{i=1}^{n}y_i$$
---

# 8.2.4 Bayesian Additive Regression Trees (BART)

In subsequent $b$ iteration, where $2\le b \le B$

Given the $k$th tree, we calculate for each observation $i$ where, $1\le i \le n$, a .blue[**partial residual**] $r_i$ is created by subtracting the response $y_i$ with the predictions made from the .blue[**other**] $K-1$ trees denoted as $(k')$. 

The other $K-1$ trees include those that have already been updated by BART $(k'<k)$ and those that will be updated by BART later $(k'>k)$.

$$r_i = y_i - \sum_{k'<k}\hat{f}^b_{k'}(x_i) - \sum_{k'>k}\hat{f}^{b-1}_{k'}(x_i)$$
---

# 8.2.4 Bayesian Additive Regression Trees (BART)

.pull-left[
Next, BART will create .blue[**possible changes/pertubations**] on the $k$th tree and assign a .blue[**probability**] calculated by Bayesian methods (given the current $k$th tree terminal nodes or $\hat{f}^{b-1}_{k}(x)$ and partial residual $r_i$).

- Adding branches 
- Pruning branches
- Change the prediction on terminal nodes
- Remain as they are

Changes that greatly improve the fit to the partial residual will be given a .blue[**higher probability**]. However, penalty scores will be given when adding new branches to prevent the tree from growing too deep. Recall that unlike boosting, we do not have control on the depth
]

.pull-right[
```{r echo = FALSE, out.width="100%"}
knitr::include_graphics("img/figure_8_12.jpg")
```
]

---

# 8.2.4 Bayesian Additive Regression Trees (BART)

BART will then "roll a dice" to .blue[**choose**] a perturbation on the $k$th tree as $\hat{f}^{b}_{k}(x)$. High probability does not mean it will get chosen.

The .blue[**selected**] tree $\hat{f}^{b}_{k}(x)$ will then we used to .blue[**fit**] the partial residual $r_i$

This same process will be done for all $K$ trees. 

Adding up all the regression tree will giving the new tree model as $\hat{f}^b(x) = \sum_{k=1}^{K}\hat{f}^b_k(x)$ for the $b$th iteration.

---

# 8.2.4 Bayesian Additive Regression Trees (BART)

Figures from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6800811/

```{r echo = FALSE, out.width="60%"}
knitr::include_graphics("img/paper1.jpg")
```

---

# 8.2.4 Bayesian Additive Regression Trees (BART)

Figures from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6800811/

.pull-left[
```{r echo = FALSE, out.width="80%"}
knitr::include_graphics("img/paper2.jpg")
```
]

.pull-right[
```{r echo = FALSE, out.width="100%"}
knitr::include_graphics("img/paper3.jpg")
```
]
---

# 8.2.4 Bayesian Additive Regression Trees (BART)

In total there will be a total of $B$ tree models

$$\hat{f}^1(x), ..., \hat{f}^B(x)$$
However, the first few models usually gives a poor prediction. Hence it is better to .blue[**get rid**] of the first few $L$ models. 

To obtain a single prediction, the final tree model involves taking the .blue[**average**] predicted value of the .blue[**remaining**] $B-L$ tree models.

$$\hat{f}(x) = \frac{1}{B-L}\sum_{b=L+1}^{B}\hat{f}^b(x)$$

---

# 8.2.4 Bayesian Additive Regression Trees (BART)

```{r echo = FALSE, out.width="60%"}
knitr::include_graphics("img/algorithm_8_3.jpg")
```