Computational Statistics – Decision Trees

`tidymodels` syntax

partition the data
build a recipe
select a model
create a workflow
fit the model
validate the model

Decision trees in action

A visual introduction to machine learning by Yee and Chu.

Yee and Chu created a step-by-step build of a recursive binary tree to model the differences between homes in SF and homes in NYC. http://www.r2d3.us/visual-intro-to-machine-learning-part-1/

Classification and Regression Trees (CART)

Basic Classification and Regression Trees (CART) Algorithm:

Start with all observations in one group.
Find the variable/split that best separates the response variable (successive binary partitions based on the different predictors / explanatory variables).
- Evaluation “homogeneity” within each group
- Divide the data into two groups (“leaves”) on that split (“node”).
- Within each split, find the best variable/split that separates the outcomes.
Continue until the groups are too small or sufficiently “pure”.
Prune tree.

Minimize heterogeneity

For every observation that falls into the region $R_{m}$ ,

prediction = the mean of the response values for observations in $R_{m}$ .

$\Rightarrow$ Minimize Residual Sum of Squares (RSS): $R S S = \sum_{m = 1}^{| T |} \sum_{i \in R_{m}} (y_{i} - {\overset{―}{y}}_{R_{m}})^{2}$ where ${\overset{―}{y}}_{R_{m}}$ is the mean response for observations within the $m$ th region.

Recursive binary splitting

Select the predictor $X_{j}$ and the cutpoint $s$ such that splitting the predictor space into the regions ${X | X_{j} < s}$ and ${X | X_{j} \geq s}$ lead to the greatest reduction in RSS.

For any $j$ and $s$ , define the pair of half-planes to be $R_{1} (j, s) = {X | X_{j} < s} and R_{2} (j, s) = {X | X_{j} \geq s}$ Find the value of $j$ and $s$ that minimize the equation: $\sum_{i : x_{i} \in R_{1} (j, s)} (y_{i} - {\overset{―}{y}}_{R_{1}})^{2} + \sum_{i : x_{i} \in R_{2} (j, s)} (y_{i} - {\overset{―}{y}}_{R_{2}})^{2}$

where ${\overset{―}{y}}_{R_{1}}$ is the mean response for observations in $R_{1} (j, s)$ and ${\overset{―}{y}}_{R_{2}}$ is the mean response observations in $R_{2} (j, s)$ .

Trees in action

Measures of impurity

${\hat{p}}_{m k}$ = proportion of observations in the $m$ th region from the $k$ th class.

classification error rate = fraction of observations in the node & not in the most common class:

$E_{m} = 1 - max_{k} ({\hat{p}}_{m k})$

Gini index $G_{m} = \sum_{k = 1}^{K} {\hat{p}}_{m k} (1 - {\hat{p}}_{m k})$
cross-entropy $D_{m} = - \sum_{k = 1}^{K} {\hat{p}}_{m k} \log {\hat{p}}_{m k}$ (Gini index & cross-entropy will both take on a value near zero if the ${\hat{p}}_{m k}$ values are all near zero or all near one.)

Recursive binary splitting

For any $j$ and $s$ , define the pair of half-planes to be $R_{1} (j, s) = {X | X_{j} < s} and R_{2} (j, s) = {X | X_{j} \geq s}$

Seek the value of $j$ and $s$ that minimize the equation: $\begin{aligned} \sum_{i : x_{i} \in R_{1} (j, s)} \sum_{k = 1}^{K} {\hat{p}}_{R_{1} k} (1 - {\hat{p}}_{R_{1} k}) + \sum_{i : x_{i} \in R_{2} (j, s)} \sum_{k = 1}^{K} {\hat{p}}_{R_{2} k} (1 - {\hat{p}}_{R_{2} k}) \\ equivalently: \\ n_{R_{1}} \sum_{k = 1}^{K} {\hat{p}}_{R_{1} k} (1 - {\hat{p}}_{R_{1} k}) + n_{R_{2}} \sum_{k = 1}^{K} {\hat{p}}_{R_{2} k} (1 - {\hat{p}}_{R_{2} k}) \end{aligned}$

Stopping

We can always make the tree more “pure” by continuing the split.

Too many splits will overfit the model to the training data!

Ways to control:

cost_complexity
tree_depth
min_n

Overfitting: http://www.r2d3.us/visual-intro-to-machine-learning-part-2/

Cost complexity

There is a cost to having a larger (more complex!) tree.

Define the cost complexity criterion, $α > 0 :$ $\begin{aligned} numerical: C_{α} (T) & = \sum_{m = 1}^{| T |} \sum_{i \in R_{m}} (y_{i} - {\overset{―}{y}}_{R_{m}})^{2} + α \cdot | T | \\ categorical: C_{α} (T) & = \sum_{m = 1}^{| T |} \sum_{i \in R_{m}} I (y_{i} \neq k (m)) + α \cdot | T | \end{aligned}$ where $k (m)$ is the class with the majority of observations in node $m$ and $| T |$ is the number of terminal nodes in the tree.

$α$ small: If $α$ is set to be small, we are saying that the risk is more worrisome than the complexity and larger trees are favored because they reduce the risk.
$α$ large: If $α$ is set to be large, then the complexity of the tree is more worrisome and smaller trees are favored.

In practice

Consider $α$ increasing. As $α$ gets bigger, the “best” tree will be smaller.

The training error will be monotonically related to the size of the training tree.

However, the test error will not be monotonically related to the size of the training tree.

A note on $α$

In the text (Introduction to Statistical Learning) and almost everywhere else you might look, the cost complexity is defined as in previous slides.

However, you might notice that in R the cost_complexity value is typically less than 1. From what I can tell, the value of the function that is being minimized in R is the average of the squared errors and the missclassification rate.

$\begin{aligned} numerical: C_{α} (T) & = \frac{1}{n} \sum_{m = 1}^{| T |} \sum_{i \in R_{m}} (y_{i} - {\overset{―}{y}}_{R_{m}})^{2} + α \cdot | T | \\ categorical: C_{α} (T) & = \frac{1}{n} \sum_{m = 1}^{| T |} \sum_{i \in R_{m}} I (y_{i} \neq k (m)) + α \cdot | T | \end{aligned}$

CART algorithm

Algorithm: Building a Regression Tree

Use recursive binary splitting to grow a large tree on the training data, stopping only when each terminal node has fewer than some minimum number of observations.
Apply cost complexity pruning to the large tree in order to obtain a sequence of best subtrees, as a function of $α$ .
Use $V$ -fold cross-validation to choose $α$ . Divide the training observations into $V$ folds, then for each $v = 1, 2, \dots, V$ :
1. Repeat Steps 1 and 2 on all but the $v$ th fold of the training data.
2. Evaluate the mean squared prediction error on the data in the left-out $k$ th fold, as a function of $α$ . For each value of $α$ , average the prediction error (either misclassification or RSS), and pick $α$ to minimize the average error.
Return the subtree from Step 2 that corresponds to the chosen value of $α$ .

CART example w defaults

penguin_cart_recipe <-
  recipe(species ~ . ,
         data = penguin_train) |>
  step_unknown(sex, new_level = "unknown") |>
  step_mutate(year = as.factor(year)) 

penguin_cart_recipe

── Recipe ──────────────────────────────────────────────────────────────────────

── Inputs

Number of variables by role

outcome:   1
predictor: 7

── Operations

• Unknown factor level assignment for: sex

• Variable mutation for: as.factor(year)

penguin_cart <- decision_tree() |>
  set_engine("rpart") |>
  set_mode("classification")

penguin_cart

Decision Tree Model Specification (classification)

Computational engine: rpart

penguin_cart_wflow <- workflow() |>
  add_model(penguin_cart) |>
  add_recipe(penguin_cart_recipe)

penguin_cart_wflow

══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: decision_tree()

── Preprocessor ────────────────────────────────────────────────────────────────
2 Recipe Steps

• step_unknown()
• step_mutate()

── Model ───────────────────────────────────────────────────────────────────────
Decision Tree Model Specification (classification)

Computational engine: rpart

penguin_cart_fit <- penguin_cart_wflow |>
  fit(data = penguin_train)

penguin_cart_fit

══ Workflow [trained] ══════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: decision_tree()

── Preprocessor ────────────────────────────────────────────────────────────────
2 Recipe Steps

• step_unknown()
• step_mutate()

── Model ───────────────────────────────────────────────────────────────────────
n= 258 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

1) root 258 139 Adelie (0.461240310 0.189922481 0.348837209)  
  2) flipper_length_mm< 206.5 164  46 Adelie (0.719512195 0.274390244 0.006097561)  
    4) bill_length_mm< 43.15 118   3 Adelie (0.974576271 0.025423729 0.000000000) *
    5) bill_length_mm>=43.15 46   4 Chinstrap (0.065217391 0.913043478 0.021739130) *
  3) flipper_length_mm>=206.5 94   5 Gentoo (0.010638298 0.042553191 0.946808511)  
    6) bill_depth_mm>=17.15 7   3 Chinstrap (0.142857143 0.571428571 0.285714286) *
    7) bill_depth_mm< 17.15 87   0 Gentoo (0.000000000 0.000000000 1.000000000) *

penguin_cart_fit |>
  predict(new_data = penguin_train) |>
  cbind(penguin_train) |>
  select(.pred_class, species) |>
  table()

           species
.pred_class Adelie Chinstrap Gentoo
  Adelie       115         3      0
  Chinstrap      4        46      3
  Gentoo         0         0     87

library(rpart.plot)
penguins_cart_plot <- 
  penguin_cart_fit |>
  extract_fit_parsnip()

rpart.plot(
  penguins_cart_plot$fit,
  roundint = FALSE)

library(rattle)
penguins_cart_plot <- 
  penguin_cart_fit |>
  extract_fit_parsnip()

fancyRpartPlot(
  penguins_cart_plot$fit,
  sub = NULL,
  palettes = "RdPu")

CART example w CV

penguin_cart_tune_recipe <-
  recipe(sex ~ body_mass_g + bill_length_mm + species,
         data = penguin_train)

set.seed(470)
penguin_vfold <- vfold_cv(penguin_train,
                          v = 5, strata = sex)

cart_grid <- expand.grid(
  cost_complexity = c(0, 10^(seq(-5,-1,1))),
  tree_depth = seq(1,6, by = 1))

cart_grid

   cost_complexity tree_depth
1            0e+00          1
2            1e-05          1
3            1e-04          1
4            1e-03          1
5            1e-02          1
6            1e-01          1
7            0e+00          2
8            1e-05          2
9            1e-04          2
10           1e-03          2
11           1e-02          2
12           1e-01          2
13           0e+00          3
14           1e-05          3
15           1e-04          3
16           1e-03          3
17           1e-02          3
18           1e-01          3
19           0e+00          4
20           1e-05          4
21           1e-04          4
22           1e-03          4
23           1e-02          4
24           1e-01          4
25           0e+00          5
26           1e-05          5
27           1e-04          5
28           1e-03          5
29           1e-02          5
30           1e-01          5
31           0e+00          6
32           1e-05          6
33           1e-04          6
34           1e-03          6
35           1e-02          6
36           1e-01          6

penguin_cart_tune <- 
  decision_tree(cost_complexity = tune(),
                tree_depth = tune()) |>
  set_engine("rpart") |>
  set_mode("classification")

penguin_cart_wflow_tune <- workflow() |>
  add_model(penguin_cart_tune) |>
  add_recipe(penguin_cart_tune_recipe)

penguin_tuned <- penguin_cart_wflow_tune |>
  tune_grid(resamples = penguin_vfold, 
           grid = cart_grid) 

penguin_tuned |> collect_metrics() |>
  filter(.metric == "accuracy") |> 
  arrange(desc(mean))

# A tibble: 36 × 8
   cost_complexity tree_depth .metric  .estimator  mean     n std_err .config   
             <dbl>      <dbl> <chr>    <chr>      <dbl> <int>   <dbl> <chr>     
 1         0                5 accuracy binary     0.847     5  0.0203 Preproces…
 2         0.00001          5 accuracy binary     0.847     5  0.0203 Preproces…
 3         0.0001           5 accuracy binary     0.847     5  0.0203 Preproces…
 4         0.001            5 accuracy binary     0.847     5  0.0203 Preproces…
 5         0                6 accuracy binary     0.847     5  0.0203 Preproces…
 6         0.00001          6 accuracy binary     0.847     5  0.0203 Preproces…
 7         0.0001           6 accuracy binary     0.847     5  0.0203 Preproces…
 8         0.001            6 accuracy binary     0.847     5  0.0203 Preproces…
 9         0                4 accuracy binary     0.843     5  0.0216 Preproces…
10         0.00001          4 accuracy binary     0.843     5  0.0216 Preproces…
# ℹ 26 more rows

Parameter choice

(Hard to see because tree_depth = 6 is on top of tree_depth = 5, either would be fine to choose.)

penguin_tuned |>
  collect_metrics() |>
  filter(.metric == "accuracy") |> 
  ggplot(aes(x = cost_complexity, y = mean, 
             color = as.factor(tree_depth))) + 
  geom_line()

penguin_tuned |> 
  select_best(metric = "accuracy")

# A tibble: 1 × 3
  cost_complexity tree_depth .config              
            <dbl>      <dbl> <chr>                
1               0          5 Preprocessor1_Model25

Best model – using CV params

penguin_opt <- decision_tree(cost_complexity = 0, 
                             tree_depth = 5) |>
  set_engine("rpart") |>
  set_mode("classification")

penguin_opt

Decision Tree Model Specification (classification)

Main Arguments:
  cost_complexity = 0
  tree_depth = 5

Computational engine: rpart

penguin_final_opt <-
  workflow() |>
  add_model(penguin_opt) |>
  add_recipe(penguin_cart_tune_recipe) |>
  fit(data = penguin_train)

model	high bias (too simple)	high variance (overfit training data)
linear model	one (zero!) predictor	many many predictors
$k$ -NN	$k = n$	$k = 1$
CART	one split	single observation in every terminal node

Decision Trees

Agenda 10/28/24

`tidymodels` syntax

Decision trees in action

Classification and Regression Trees (CART)

Minimize heterogeneity

Recursive binary splitting

Trees in action

Measures of impurity

Recursive binary splitting

Stopping

Cost complexity

In practice

A note on $α$

CART algorithm

CART example w defaults

Plotting the tree (not tidy)

CART example w CV

Parameter choice

Best model – using CV params

Best model – using `finalize_model()`

Best Model Predictions

Best Model Predictions

Bias-Variance Tradeoff

Bias-Variance Tradeoff

Bias-Variance Tradeoff

Reflecting on Model Building

Reflecting on Model Building

Reflecting on Model Building