---
title: "Machine learning"
subtitle: "Tree-based models"
author: "Ron Yurko"
date: "07/13/2020"
---

```{r setup, echo = FALSE}
knitr::opts_chunk$set(fig.retina = 3, warning = FALSE, message = FALSE)
```

## What is Machine Learning?

--
The short version:

- Machine learning (ML) is a subset of statistical learning that focuses on prediction.


--
The longer version:

- ML is the idea of constructing data-driven algorithms that *learn* the mapping between predictor variables and response variable(s). 

 - We do not assume a parametric form for the mapping *a priori*, even if technically one can write one down *a posteriori* (e.g., by translating a tree model to a indicator-variable mathematical expression)

- Linear regression, for instance, is not a ML algorithm since we can write down the linear equation ahead of time, but random forest is a ML algorithm since we've no idea what the number of splits will end up being in each of its individual trees.

---

## Which algorithm is best?

__That's not the right question to ask.__

(And the answer is *not* deep learning. Because if the underlying relationship between your predictors and your response is truly linear, *you do not need to apply deep learning*! Just do linear regression. Really. It's OK.)


--
The right question is ask is: __why should I try different algorithms?__


--
The answer to that is that without superhuman powers, you cannot visualize the distribution of predictor variables in their native space. 

  - Of course, you can visualize these data *in projection*, for instance when we perform EDA
  
  - And the performance of different algorithms will depend on how predictor data are distributed...

---

## Data geometry

```{r out.width='70%', echo = FALSE, fig.align='center'}
knitr::include_graphics("http://www.stat.cmu.edu/~pfreeman/data_geometry.png")
```

The picture above shows data for which there are two predictor variables (along the x-axis and the y-axis) and for which the response variable is binary: x's and o's. An algorithm that utilizes linear boundaries or segments the plane into rectangles will do well given the data to the left, whereas an algorithm that utilizes circular boundaries will fare better given the data to the right.

"do well/fare better": will do a better job at predicting whether a new observation is actually an x or an o.

---

## Decision trees

Decision trees partition training data into __homogenous nodes / subgroups__ with similar response values


--
The subgroups are found __recursively using binary partitions__ 

  - i.e. asking a series of yes-no questions about the predictor variables

We stop splitting the tree once a __stopping criteria__ has been reached (e.g. maximum depth allowed)


--
For each subgroup / node predictions are made with:

- Regression tree: __the average of the response values__ in the node

- Classification tree: __the most popular class__ in the node


--
Most popular approach is Leo Breiman's __C__lassification __A__nd __R__egression __T__ree (CART) algorithm

---

## Decision tree structure

```{r out.width='100%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://bradleyboehmke.github.io/HOML/images/decision-tree-terminology.png")
```

---

## Decision tree structure

We make a prediction for an observation by __following its path along the tree__


```{r out.width='100%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://bradleyboehmke.github.io/HOML/images/exemplar-decision-tree.png")
```


--

- Decision trees are __very easy to explain__ to non-statisticians.

- Easy to visualize and thus easy to interpret __without assuming a parametric form__

---

## Binary recursive partitioning

CART uses __recursive__ splits where each _split / rule_ depends on the previous split / rule _above_ it


--
__Objective at each split__: find the __best__ variable to partition the data into one of two regions, $R_1$ & $R_2$, to __minimize the error__ between the actual response, $y_i$, and the node's predicted constant, $c_i$

- For regression we minimize the sum of squared errors (SSE):

$$S S E=\sum_{i \in R_{1}}\left(y_{i}-c_{1}\right)^{2}+\sum_{i \in R_{2}}\left(y_{i}-c_{2}\right)^{2}$$


--
- For classification trees we minimize the node's _impurity_ the __Gini index__

  - where $p_k$ is the proportion of observations in the node belonging to class $k$ out of $K$ total classes
  
  - want to minimize $Gini$: small values indicate a node has primarily one class (_is more pure_)

$$Gini = 1 - \sum_k^K p_k^2$$


--
Splits yield __locally optimal__ results, so we are NOT guaranteed to train a model that is globally optimal


--
_How do we control the complexity of the tree?_

---

## Tune the __maximum tree depth__ or __minimum node size__


```{r out.width='60%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://bradleyboehmke.github.io/HOML/07-decision-trees_files/figure-html/dt-early-stopping-1.png")
```

---

## Prune the tree by tuning __cost complexity__

Can grow a very large complicated tree, and then __prune__ back to an optimal __subtree__ using a __cost complexity__ parameter $\alpha$ (like $\lambda$ for elastic net) 

- $\alpha$ penalizes objective as a function of the number of __terminal nodes__

- e.g., we want to minimize $SSE + \alpha \cdot (\# \text{  of terminal nodes})$ 

```{r out.width='80%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://bradleyboehmke.github.io/HOML/07-decision-trees_files/figure-html/pruned-tree-1.png")
```

---

## Example data: MLB 2019 batting statistics

Downloaded MLB 2019 batting statistics leaderboard from [Fangraphs](https://www.fangraphs.com/leaders.aspx?pos=all&stats=bat&lg=all&qual=y&type=8&season=2019&month=0&season1=2019&ind=0)

```{r load-data, warning = FALSE, message = FALSE}
library(tidyverse)
mlb_data <- read_csv("http://www.stat.cmu.edu/cmsac/sure/materials/data/fg_batting_2019.csv") %>%
  janitor::clean_names() %>%
  mutate_at(vars(bb_percent:k_percent), parse_number)
head(mlb_data)
```


---

## Regression tree example with the [`rpart` package](https://cran.r-project.org/web/packages/rpart/vignettes/longintro.pdf)

Revisit the modeling of `w_oba` from the [KNN slides](http://www.stat.cmu.edu/cmsac/sure/materials/lectures/slides/15-knn-kernel.html#21)

```{r init-rpart}
library(rpart)
init_mlb_tree <- rpart(formula = w_oba ~ bb_percent + k_percent + iso,
                       data = mlb_data, method  = "anova")
init_mlb_tree
```


---

## Display the tree with [`rpart.plot`](plhttp://www.milbo.org/rpart-plot/)



.pull-left[

```{r plot-tree, fig.align ='center', fig.height=6}
library(rpart.plot)
rpart.plot(init_mlb_tree)
```

]

.pull-right[

- `rpart()` runs 10-fold CV to tune $\alpha$ for pruning 

- Selects # terminal nodes via 1 SE rule

```{r plot-complexity, fig.align ='center', fig.height=5}
plotcp(init_mlb_tree)
```


]


---

## What about the full tree? (check out `rpart.control`)

.pull-left[

```{r plot-full-tree, fig.align ='center', fig.height=4}
full_mlb_tree <- rpart(formula = w_oba ~ bb_percent + k_percent + iso,
                       data = mlb_data, method  = "anova", 
                       control = list(cp = 0, xval = 10))
rpart.plot(full_mlb_tree)
```

]

.pull-right[


```{r plot-full-complexity, fig.align ='center', fig.height=4}
plotcp(full_mlb_tree)
```


]

---

## Train with `caret`

```{r caret-tree, fig.align ='center', fig.height=5}
library(caret)
caret_mlb_tree <- train(w_oba ~ bb_percent + k_percent + iso + avg + obp + slg + war,
                        data = mlb_data, method = "rpart",
                        trControl = trainControl(method = "cv", number = 10),
                        tuneLength = 20)
ggplot(caret_mlb_tree) + theme_bw()
```


---

## Display the final model

```{r, fig.align ='center', fig.height=6}
rpart.plot(caret_mlb_tree$finalModel)
```

---

## Preview of summarizing variables in tree-based models

.pull-left[
__Variable importance__ - based on reduction in SSE (_notice anything odd?_)

```{r var-imp, fig.align ='center', fig.height=4}
library(vip)
vip(caret_mlb_tree, geom = "point") + theme_bw()
```


]
.pull-right[

- Summarize single variable's relationship with __partial dependence plot__

```{r pdp, fig.align ='center', fig.height=4}
library(pdp)
partial(caret_mlb_tree, pred.var = "obp") %>% autoplot() + theme_bw()

```


]

