Elements of Data Science
SDS 322E

H. Sherry Zhang
Department of Statistics and Data Sciences
University of Texas at Austin

Get the repository for today:
library(usethis)
create_from_github("SDS322E-26FALL/1202-trees", fork = FALSE)

Learning objectives

  • Understand the basic idea of decision trees for regression and classification
  • Fit decision trees using tidymodels for both regression and classification tasks
  • Identify the key parameters of decision trees and their effects
  • Understand the idea of random forest and how it improves over decision trees
  • Apply random forest using tidymodels (with hyper-parameter tuning)

Decision trees

Classification tree

Can we predict penguin sex from bill_len and bill_dep with a few simple rules?

A tree is a set of yes/no questions. Start at the top and follow the answers down:

  • Is bill_len < 48? No (right): predict male
  • Yes (left): is bill_dep < 18? …

Each leaf at the bottom is a prediction: the majority class of the training penguins that end up there.

Classification tree

  • Divide/partition the predictor space \(X_1, \ldots, X_p\) into \(J\) regions \(R_1, \ldots, R_J\).

  • For every observation that falls within region \(R_j\), we make the same prediction: the majority class of the training observations in \(R_j\).

  • The goal is to increase the “purity” of each region, so that most regions are dominated by a single class: choose the splits that minimize the impurity (Gini index or entropy, see the appendix).

Each rectangle is one leaf of the tree: the splits are always parallel to the axes.

Regression tree

  • The same partition of the predictor space into regions \(R_1, \ldots, R_J\), but the response is now numeric.

  • For every observation that falls within region \(R_j\), we predict \(\bar{y}_{R_j}\): the mean of the response values for the training observations in \(R_j\) (instead of the majority class).

  • Purity doesn’t apply to a numeric response, so we choose the splits that minimize the residual sum of squares instead: \(RSS = \sum_{j=1}^{J} \sum_{i \in R_j} (y_i - \bar{y}_{R_j})^2\)

Trees for regression and classification

The same idea works for both: split the predictor space into rectangles, and make one prediction per rectangle.

Regression tree Classification tree
Response numeric, e.g. body_mass categorical, e.g. sex
Prediction in a region the mean of the training responses the majority class of the training responses
Choose the splits to minimize the residual sum of squares (RSS) the impurity (Gini index or entropy)
tidymodels decision_tree(mode = "regression") decision_tree(mode = "classification")

The same is true for KNN: average the responses of the \(k\) nearest neighbors for regression, or take their majority vote for classification, with nearest_neighbor(mode = "regression") or nearest_neighbor(mode = "classification").

How big should a tree be?

A tree keeps splitting until the regions are pure or too small. The hyperparameter min_n is the minimum number of penguins a region must have to be split again:

  • Small min_n: many small regions that follow individual penguins: a very flexible tree, like a small \(k\) in KNN
  • Large min_n: a few big regions: a simple tree
  • Like \(k\) in KNN, min_n is a hyperparameter, so we can choose it with cross-validation

Here we set cost_complexity = 0, so that min_n alone controls the size of the tree. cost_complexity and tree_depth are two other ways to limit the size of a tree.

Random forest

A bootstrap sample

A bootstrap sample is a new dataset of the same size, drawn from the original data with replacement: pick a row at random, put it back, and repeat.

For example, with 6 penguins:

Original data 1 2 3 4 5 6
Bootstrap sample 2 2 4 5 5 6
  • Penguins 2 and 5 appear twice; penguins 1 and 3 are left out
  • Each bootstrap sample is a slightly different version of the data, with no new data needed - A repeated penguin stands in for similar penguins we never measured.
  • On average, a bootstrap sample contains about 63% of the original rows; the rest are repeats

One tree is unstable

Fit the same classification tree for sex (min_n = 15) on three bootstrap samples of the 333 penguins:

Small changes in the data give very different trees: a single tree has high variance, just like KNN with \(k = 1\).

Lesson from KNN: averaging over many neighbors gave a much more stable prediction. Can we average over many trees?

Grow a forest

A random forest grows many trees (default 500), each one a little different, and lets them vote:

  1. Random rows: each tree is fit on its own bootstrap sample of the training data (as on the previous slide)
  2. Random predictors: at each split, a tree may only choose from a random subset of mtry predictors

Why step 2? If one predictor is very strong, every tree would split on it first, and all the trees would look alike: averaging similar trees doesn’t reduce the variance much. Hiding it from some splits makes the trees different, so their average is more stable.

To predict, combine the trees:

  • majority vote for classification: rand_forest(mode = "classification") or
  • the average for regression: rand_forest(mode = "regression").

Combining 500 trees: majority vote

A new penguin with bill_len = 41 and bill_dep = 16: ask every tree in the forest. Tree \(k\) is the tree (min_n = 15) on bootstrap sample \(k\)

Vote
Tree 1 female
Tree 2 female
Tree 3 male
Tree 4 male
Tree 5 female
Tree 6 male
… …
Tree 500 male

Final prediction: female, the majority: 287 of the 500 trees vote female, even though some trees disagree

From individual trees to the final prediction

Do this vote at every point of the space: trees 1 to 3 (the same trees as on bootstrap samples 1 to 3), and the final prediction of all 500 trees

  • Each tree is blocky and makes mistakes in different places: each has its own bootstrap sample
  • Far from the boundary, almost every tree agrees. Near it, the trees disagree, and the vote settles it
  • The final prediction is more stable than any single tree: changing a few trees barely changes the vote

Hyperparameters of a random forest

  • min_n: how deep each tree grows, as in a single tree
    • small: deep, flexible trees; large: shallow, simple trees
    • tune it
  • mtry (new): how many predictors each split can choose from
    • small: the trees are very different; large: the trees are very similar
    • tune it
  • trees (new): how many trees to average
    • few: unstable predictions; many: stable, but slower to fit
    • don’t tune it: just use enough (the default, 500)

trees: more trees never hurt

Example (as in the tuning workflow): predict body_mass from the other six penguin variables, and compute the RMSE on the testing data

  • One tree: RMSE 400 g. Averaging more trees quickly brings it down, and after about 100 trees it stays flat at about 311 g
  • Adding trees never makes the forest overfit: it only takes longer to fit. So we don’t tune trees, we just use enough of them (the default, 500, is plenty here)

mtry: how different are the trees?

Which predictor does each of the 500 trees split on first?

  • mtry = 6 (all predictors): every tree picks the strongest predictors (species, flipper_len) first, so the trees are near-copies, and averaging copies doesn’t reduce the variance much
  • Small mtry: the trees are forced to use different predictors, so they are different, and their average is more stable. But if mtry is too small, each tree is too weak
  • The best mtry is somewhere in between, and depends on the data: tune it with cross-validation

Random forest workflow with tuning and CV

  • Step 0: Training-testing split

  • Step 1: Specify the model, with mtry and min_n left to tune

  • Step 2: Specify the pre-processing recipe

  • Step 3: Build a workflow

  • Step 4: Fit the model

    • 4a: Generate the cross validation folds from the training data
    • 4b: Generate the hyperparameter grid
    • 4c: Fit the model on the cv folds and hp grid
    • 4d: Average the metrics over the folds
    • 4e: Find the best model
    • 4f: Plug in the best values and refit on the whole training data

  • Step 5/6: Predict and calculate prediction accuracy metrics on the testing data

Random forest workflow with tuning and CV

Predict body_mass from the other six penguin variables.

Step 0: Training-testing split:

penguins_clean <- as_tibble(datasets::penguins) |> na.omit() |> select(-year)
set.seed(123)
penguins_split <- initial_split(penguins_clean, prop = 0.75)
penguins_train <- training(penguins_split)
penguins_test  <- testing(penguins_split)

Step 1/2/3: Specify the model, with mtry = tune() and min_n = tune(), the pre-processing recipe, and the workflow

rf_tune_spec <- rand_forest(mode = "regression", trees = 500, mtry = tune(), min_n = tune()) |>
  set_engine("ranger")
rf_recipe <- recipe(body_mass ~ ., data = penguins_train)
rf_tune_wf <- workflow() |> add_model(rf_tune_spec) |> add_recipe(rf_recipe)

Random forest workflow with tuning and CV

Step 4a: Create the cross-validation folds from the training data

set.seed(123)
penguins_folds <- vfold_cv(penguins_train, v = 10)
penguins_folds
#  10-fold cross-validation 
# A tibble: 10 × 2
   splits           id    
   <list>           <chr> 
 1 <split [224/25]> Fold01
 2 <split [224/25]> Fold02
 3 <split [224/25]> Fold03
 4 <split [224/25]> Fold04
 5 <split [224/25]> Fold05
 6 <split [224/25]> Fold06
 7 <split [224/25]> Fold07
 8 <split [224/25]> Fold08
 9 <split [224/25]> Fold09
10 <split [225/24]> Fold10
  • Same as for KNN: one row per fold, each fold is fit on 9/10 of the training data and validated on the other 1/10

Random forest workflow with tuning and CV

Step 4b: Set up the hyperparameter grid: now two hyperparameters

rf_grid <- grid_regular(mtry(range = c(1, 6)), min_n(range = c(2, 20)), levels = c(6, 4))
rf_grid |> print(n = 8)
# A tibble: 24 × 2
   mtry min_n
  <int> <int>
1     1     2
2     2     2
3     3     2
4     4     2
5     5     2
6     6     2
7     1     8
8     2     8
# ℹ 16 more rows
  • One column per hyperparameter we set to tune(): mtry and min_n
  • levels = c(6, 4): 6 values of mtry (1 to 6) \(\times\) 4 values of min_n (2, 8, 14, 20) = 24 combinations: every value of mtry is paired with every value of min_n

Random forest workflow with tuning and CV

Step 4c: Fit the forest for every combination on every fold

rf_tune_res <- tune_grid(rf_tune_wf, resamples = penguins_folds, grid = rf_grid)
rf_tune_res
# Tuning results
# 10-fold cross-validation 
# A tibble: 10 × 4
   splits           id     .metrics          .notes          
   <list>           <chr>  <list>            <list>          
 1 <split [224/25]> Fold01 <tibble [48 × 6]> <tibble [0 × 4]>
 2 <split [224/25]> Fold02 <tibble [48 × 6]> <tibble [0 × 4]>
 3 <split [224/25]> Fold03 <tibble [48 × 6]> <tibble [0 × 4]>
 4 <split [224/25]> Fold04 <tibble [48 × 6]> <tibble [0 × 4]>
 5 <split [224/25]> Fold05 <tibble [48 × 6]> <tibble [0 × 4]>
 6 <split [224/25]> Fold06 <tibble [48 × 6]> <tibble [0 × 4]>
 7 <split [224/25]> Fold07 <tibble [48 × 6]> <tibble [0 × 4]>
 8 <split [224/25]> Fold08 <tibble [48 × 6]> <tibble [0 × 4]>
 9 <split [224/25]> Fold09 <tibble [48 × 6]> <tibble [0 × 4]>
10 <split [225/24]> Fold10 <tibble [48 × 6]> <tibble [0 × 4]>

Question: how many trees are we fitting here?

24 combinations \(\times\) 10 folds = 240 random forests, and each forest has 500 trees: 240 \(\times\) 500 = 120,000 trees!

  • .metrics: for each fold, the validation metrics (rmse, rsq) of every combination

Random forest workflow with tuning and CV

Step 4d: Average the metrics over the 10 folds

collect_metrics(rf_tune_res) |> filter(.metric == "rmse")
# A tibble: 24 × 8
    mtry min_n .metric .estimator  mean     n std_err .config         
   <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>           
 1     1     2 rmse    standard    305.    10    10.9 pre0_mod01_post0
 2     1     8 rmse    standard    305.    10    10.0 pre0_mod02_post0
 3     1    14 rmse    standard    311.    10    10.1 pre0_mod03_post0
 4     1    20 rmse    standard    311.    10    10.3 pre0_mod04_post0
 5     2     2 rmse    standard    286.    10    11.0 pre0_mod05_post0
 6     2     8 rmse    standard    286.    10    10.8 pre0_mod06_post0
 7     2    14 rmse    standard    288.    10    10.9 pre0_mod07_post0
 8     2    20 rmse    standard    291.    10    10.8 pre0_mod08_post0
 9     3     2 rmse    standard    287.    10    11.1 pre0_mod09_post0
10     3     8 rmse    standard    287.    10    11.6 pre0_mod10_post0
# ℹ 14 more rows
  • One row per combination of mtry and min_n: mean is the average validation RMSE over the n = 10 folds

Random forest workflow with tuning and CV

Step 4d: Average the metrics over the 10 folds

Plot the cross validation RMSE for each combination 
best_rf <- select_best(rf_tune_res, metric = "rmse")
rf_tune_res |>
  collect_metrics() |>
  filter(.metric == "rmse") |>
  ggplot(aes(x = mtry, y = mean, color = factor(min_n))) +
  geom_line() +
  geom_point() +
  geom_point(data = \(d) filter(d, mtry == best_rf$mtry, min_n == best_rf$min_n),
             size = 5, shape = 21, stroke = 1.5, color = "black") +
  scale_x_continuous(breaks = 1:6) +
  labs(x = "mtry (predictors considered at each split)", y = "Cross validation RMSE (g)",
       color = "min_n")

Lower RMSE is better: the circled combination has the lowest average RMSE.

Random forest workflow with tuning and CV

Step 4e: Find the best combination

show_best(rf_tune_res, metric = "rmse", n = 5)
# A tibble: 5 × 8
   mtry min_n .metric .estimator  mean     n std_err .config         
  <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>           
1     2     2 rmse    standard    286.    10    11.0 pre0_mod05_post0
2     2     8 rmse    standard    286.    10    10.8 pre0_mod06_post0
3     3    14 rmse    standard    287.    10    11.2 pre0_mod11_post0
4     3     2 rmse    standard    287.    10    11.1 pre0_mod09_post0
5     4     2 rmse    standard    287.    10    11.6 pre0_mod13_post0
best_rf <- select_best(rf_tune_res, metric = "rmse")
best_rf
# A tibble: 1 × 3
   mtry min_n .config         
  <int> <int> <chr>           
1     2     2 pre0_mod05_post0
  • show_best(): the top combinations, sorted by mean RMSE (lowest first)
  • select_best(): a one-row tibble with the best mtry and min_n, ready to plug into the workflow

Random forest workflow with tuning and CV

Step 4f: Plug in the best values and refit on the whole training data

final_rf_wf <- finalize_workflow(rf_tune_wf, best_rf)
final_rf_wf
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: rand_forest()

── Preprocessor ────────────────────────────────────────────────────────────────
0 Recipe Steps

── Model ───────────────────────────────────────────────────────────────────────
Random Forest Model Specification (regression)

Main Arguments:
  mtry = 2
  trees = 500
  min_n = 2

Computational engine: ranger
  • The tune() placeholders are replaced by mtry = 2 and min_n = 2: now an ordinary workflow that we can fit

Random forest workflow with tuning and CV

Step 4f: Plug in the best values and refit on the whole training data

set.seed(1)
final_rf_fit <- fit(final_rf_wf, data = penguins_train)
  • The 240 forests in step 4c were only used to choose the hyperparameters: the final forest is fit once on all of the training data

Random forest workflow with tuning and CV

Step 5/6: Predict and calculate prediction accuracy metrics on the testing data, once

res_final <- augment(final_rf_fit, penguins_test)
metrics(res_final, truth = body_mass, estimate = .pred)
# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard     313.   
2 rsq     standard       0.842
3 mae     standard     255.

last_fit(final_rf_wf, split = penguins_split) does the refit and the testing-data evaluation in one step.