Decision Trees and Random Forests

1. A decision tree is a sequence of questions

A tree might ask “Is age above a certain value?” then “Is BMI above another value?” Each answer moves to another branch. Tree split thresholds are learned from training data. A single deep tree can memorize details and overfit.

Example: Write a tiny tree by hand as case_when() rules and run it on three typed patients. The cutoffs 136.5 and 60.5 are the ones the fitted tree in section 2 will learn.

library(dplyr)
patients <- tibble(age = c(45, 70, 66), sbp = c(120, 150, 140))
patients |>
  mutate(guess = case_when(
    sbp < 136.5  ~ "No",     # first question: is SBP low?
    age >= 60.5  ~ "Yes",    # second question: is age high?
    TRUE         ~ "No"
  ))

Output

# A tibble: 3 x 3
    age   sbp guess
  <dbl> <dbl> <chr>
1    45   120 No   
2    70   150 Yes  
3    66   140 Yes

The first patient has a low SBP and stops at the first question (No). The other two have high SBP, so the second question (age) decides: both are over 60.5 and get Yes. A 45-year-old with SBP 150 would have got No. A fitted tree is just a longer list of questions like this one, with cutoffs chosen from the training rows.

A random forest fits many randomized trees and aggregates their outputs. It can learn nonlinear patterns and interactions; it is not inherently a causal explanation or a better model for every dataset.

2. Prepare the same fair comparison

Both algorithms should be evaluated with the same training/test definition and, ideally, matched CV folds and metrics. We compare models on training rows only.

library(tidymodels)
ehr <- readr::read_csv("data/synthetic_ehr.csv", show_col_types = FALSE) |>
  mutate(diabetes = factor(diabetes, levels = c("Yes", "No")),
         sex = factor(sex), clinic_site = factor(clinic_site),
         family_history = factor(family_history))
set.seed(500)
split <- initial_split(ehr, prop = 0.8, strata = diabetes)
train <- training(split)
rec <- recipe(diabetes ~ age + bmi + sbp + pulse + family_history + sex + clinic_site,
              data = train) |>
  step_impute_median(all_numeric_predictors()) |>
  step_impute_mode(all_nominal_predictors()) |>
  step_dummy(all_nominal_predictors()) |>
  step_zv(all_predictors())
set.seed(500)
folds <- vfold_cv(train, v = 5, strata = diabetes)

Example: Fit one small tree on age, bmi, and sbp using the training rows, and print its rules.

tree_fit <- decision_tree(tree_depth = 3, cost_complexity = 0.02, min_n = 10) |>
  set_engine("rpart") |>
  set_mode("classification") |>
  fit(diabetes ~ age + bmi + sbp, data = train)
extract_fit_engine(tree_fit)

Output

n= 480 

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

1) root 480 45 No (0.09375000 0.90625000)  
  2) sbp>=136.5 70 26 No (0.37142857 0.62857143)  
    4) bmi>=30.55 51 24 No (0.47058824 0.52941176)  
      8) age>=60.5 43 19 Yes (0.55813953 0.44186047) *
      9) age< 60.5 8  0 No (0.00000000 1.00000000) *
    5) bmi< 30.55 19  2 No (0.10526316 0.89473684) *
  3) sbp< 136.5 410 19 No (0.04634146 0.95365854) *

The tree starts with all 480 rows (45 are Yes). Its first question is sbp>=136.5: 70 rows answer yes and 410 answer no, and the 410 low-SBP rows are called No (only 19 of them are Yes). The 70 high-SBP rows are split again by BMI and age. The one leaf marked Yes holds 43 rows (SBP at least 136.5, BMI at least 30.55, age at least 60.5), of which 24 are truly Yes. These cutoffs came from train, and folds is ready for the model comparison in the next sections.

3. Create a forest

forest_spec <- rand_forest(trees = 300, mtry = 3, min_n = 5) |>
  set_engine("ranger", importance = "permutation") |>
  set_mode("classification")
forest_wf <- workflow() |>
  add_recipe(rec) |>
  add_model(forest_spec)

Read the arguments:

trees = 300: build 300 trees.
mtry = 3: consider three randomly selected predictors when splitting a node.
min_n = 5: require a node to hold at least this many observations before it can be split further, helping control complexity.
importance = "permutation": request a variable-importance summary, which is not causal effect.

Example: Fit this forest once on all of train. ranger also scores each tree on the rows that tree did not see (the out-of-bag rows), which gives a quick error rate.

set.seed(500)
forest_fit <- fit(forest_wf, data = train)
forest_fit |> extract_fit_engine() |> pluck("prediction.error")

Output

[1] 0.07718913

The out-of-bag error is about 0.077, so about 7.7% of the training rows are misclassified by trees that did not see them. Always answering No would be wrong on 9.4% of rows (45 of 480), so the forest is only modestly better than doing nothing at the default cutoff.

4. Cross-validate the forest

set.seed(500)
forest_cv <- fit_resamples(
  forest_wf, resamples = folds,
  metrics = metric_set(roc_auc, pr_auc),
  control = control_resamples(save_pred = TRUE)
)
collect_metrics(forest_cv)

Output

# A tibble: 2 x 6
  .metric .estimator  mean     n std_err .config        
  <chr>   <chr>      <dbl> <int>   <dbl> <chr>          
1 pr_auc  binary     0.425     5  0.0531 pre0_mod0_post0
2 roc_auc binary     0.807     5  0.0717 pre0_mod0_post0

Compare the results with Chapter 16’s logistic model using the same folds. Tiny score differences may be sampling noise.

Example: Read the collect_metrics() table above. The forest’s mean ROC AUC across the five folds is 0.807 and its PR AUC is 0.425. Chapter 16’s logistic model on the same folds had 0.829 and 0.322, so the forest is a bit worse on ROC AUC (gap 0.022) and better on PR AUC (gap 0.103). The forest’s standard errors are 0.072 for ROC AUC and 0.053 for PR AUC. The ROC AUC gap is smaller than that standard error, so it is not convincing; the PR AUC gap is about twice its standard error, which is more notable. Each standard error describes one model’s mean, though, so to call a winner you would compare the two models fold by fold on the same folds.

5. Tune one hyperparameter

Rather than a giant search that is hard to explain, first compare three possible min_n values using training folds.

tuned_spec <- rand_forest(trees = 300, mtry = 3, min_n = tune()) |>
  set_engine("ranger", importance = "permutation") |>
  set_mode("classification")
tuned_wf <- workflow() |>
  add_recipe(rec) |>
  add_model(tuned_spec)
set.seed(500)
tuning <- tune_grid(
  tuned_wf,
  resamples = folds,
  grid = tibble(min_n = c(2L, 5L, 10L)),
  metrics = metric_set(roc_auc, pr_auc)
)
collect_metrics(tuning)

Output

# A tibble: 6 x 7
  min_n .metric .estimator  mean     n std_err .config        
  <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>          
1     2 pr_auc  binary     0.414     5  0.0496 pre0_mod1_post0
2     2 roc_auc binary     0.795     5  0.0770 pre0_mod1_post0
3     5 pr_auc  binary     0.412     5  0.0372 pre0_mod2_post0
4     5 roc_auc binary     0.793     5  0.0756 pre0_mod2_post0
5    10 pr_auc  binary     0.431     5  0.0499 pre0_mod3_post0
6    10 roc_auc binary     0.804     5  0.0816 pre0_mod3_post0

best <- select_best(tuning, metric = "roc_auc")
best

Output

# A tibble: 1 x 2
  min_n .config        
  <int> <chr>          
1    10 pre0_mod3_post0

chosen_forest_workflow <- finalize_workflow(tuned_wf, best)

This is the selection stage, not a final independent evaluation. Tuning too many choices against the same folds eventually overfits them too. We will evaluate only the chosen workflow once on test in the capstone.

Example: Read the tuning table and best. The mean ROC AUC is 0.795 for min_n = 2, 0.793 for 5, and 0.804 for 10, so select_best() picked min_n = 10. The three values differ by 0.011, far less than their standard errors of about 0.08. The choice of 10 is therefore a weak preference, not proof that 10 is truly better. chosen_forest_workflow is the workflow with min_n = 10 filled in, ready to be fitted later.

6. Model interpretation and safety

Variables assigned high importance are not necessarily causal, fair, stable, or even available at deployment. Correlated predictors share information; rankings can change depending on the data. Check leakage, calibration, subgroup performance, and errors. Neither tree ensembles nor logistic regression should be deployed as a clinical decision aid without external validation, privacy review, and oversight.

Example: Print the permutation importance of the forest we fitted in section 3, largest first.

forest_fit |>
  extract_fit_engine() |>
  pluck("variable.importance") |>
  sort(decreasing = TRUE) |>
  round(4)

Output

              sbp               bmi               age clinic_site_North 
           0.0135            0.0070            0.0050            0.0022 
family_history_X1             pulse clinic_site_South             sex_M 
           0.0015            0.0012            0.0010           -0.0018

Shuffling sbp hurts the forest’s predictions most (0.0135), then bmi (0.0070) and age (0.0050). Shuffling sex_M has a negative value (-0.0018), which means shuffling it did not hurt at all: its importance is near zero. A high value says the forest leaned on that column for prediction; it does not say that changing the measurement would change the outcome.

7. Where do PCA and clustering fit?

PCA compresses correlated measurements into combinations; K-means groups observations by similarity. These are unsupervised techniques and are not substitutes for supervised validation. Part IV opens with them in Chapter 23, after you have seen labels, folds, and inference risk in Parts II and III. They also require scale/missingness decisions. See the governance appendix.

Example: Run PCA on the four vital-sign columns, with the rows that have no missing values and every column scaled to the same spread.

vitals <- ehr |> select(age, bmi, sbp, pulse) |> tidyr::drop_na()
nrow(vitals)

Output

[1] 517

round(summary(prcomp(vitals, scale. = TRUE))$importance, 3)

Output

                         PC1   PC2   PC3   PC4
Standard deviation     1.354 1.015 0.988 0.399
Proportion of Variance 0.458 0.258 0.244 0.040
Cumulative Proportion  0.458 0.716 0.960 1.000

517 of the 600 rows have all four values. The first component captures 45.8% of the variation and the first two together 71.6%. No diabetes label was used anywhere, so the result describes how the measurements vary together, not who has the condition.

Try it

Compare CV ROC AUC for fixed min_n=5 with other tried sizes.

Show solution

Answer: The fixed forest with min_n = 5 has a CV ROC AUC of 0.807. The tuning runs gave 0.795 (min_n = 2), 0.793 (5), and 0.804 (10). The differences are about 0.01, which is noise. The logistic model on the same folds (0.829) is the highest, though not clearly better.

Step by step

Collect the ROC AUC of the fixed forest and of the three tuned values.

fixed <- collect_metrics(forest_cv) |>
  filter(.metric == "roc_auc") |>
  transmute(model = "fixed forest, min_n = 5", mean, std_err)
tried <- collect_metrics(tuning) |>
  filter(.metric == "roc_auc") |>
  transmute(model = paste("tuning, min_n =", min_n), mean, std_err)
bind_rows(fixed, tried)

Output

# A tibble: 4 x 3
  model                    mean std_err
  <chr>                   <dbl>   <dbl>
1 fixed forest, min_n = 5 0.807  0.0717
2 tuning, min_n = 2       0.795  0.0770
3 tuning, min_n = 5       0.793  0.0756
4 tuning, min_n = 10      0.804  0.0816

2. The min_n = 5 forest appears twice, with 0.807 and 0.793. Same setting, different numbers: a forest uses random draws. Both runs used the same folds and the same seed, but tune_grid() fits min_n = 2 first, so the min_n = 5 forest most likely gets different random draws than in fit_resamples() (we did not test this). That gap (0.014) shows how much random variation to expect. 3. For reference, fit the logistic model on the same folds.

logistic_wf <- workflow() |>
  add_recipe(rec) |>
  add_model(logistic_reg() |> set_engine("glm") |> set_mode("classification"))
fit_resamples(logistic_wf, resamples = folds, metrics = metric_set(roc_auc)) |>
  collect_metrics() |>
  select(.metric, mean, std_err)

Output

# A tibble: 1 x 3
  .metric  mean std_err
  <chr>   <dbl>   <dbl>
1 roc_auc 0.829  0.0447

4. Logistic regression has 0.829, the highest of the five numbers. Every standard error is 0.04 to 0.08, so none of these models is clearly better; the simpler logistic model is a reasonable choice.

Explain the difference between a tree and a forest.

Show solution

Answer: A tree is one set of yes/no questions that sends each row to a single answer. A forest builds many such trees, each on a random sample of rows and with a random subset of predictors at each split, and combines their votes. One tree is easy to read but can memorize the training data; the forest averages out that memorization.

Step by step

A single tree, tree_fit from section 2, gives each row one leaf. Count how many leaves it has by listing its predicted probabilities.

tree_prob <- predict(tree_fit, new_data = train, type = "prob")$.pred_Yes
n_distinct(tree_prob)

Output

[1] 4

2. The tree can output only 4 different probabilities, one per leaf, however many rows it scores. 3. The forest averages 300 trees, so its scores are much finer.

forest_prob <- predict(forest_fit, new_data = train, type = "prob")$.pred_Yes
n_distinct(forest_prob)

Output

[1] 409

4. The forest gives many different probabilities for the 480 rows. The price is that you can no longer read one set of rules from it.

State why high feature importance does not establish causality.

Show solution

Answer: Importance only measures how much the forest’s predictions get worse when a column is shuffled. A column can matter for prediction because it is correlated with the real cause, or with another predictor, without causing anything itself.

Step by step

Print the importance of the forest we fitted.

forest_fit |>
  extract_fit_engine() |>
  pluck("variable.importance") |>
  sort(decreasing = TRUE) |>
  round(4)

Output

              sbp               bmi               age clinic_site_North 
           0.0135            0.0070            0.0050            0.0022 
family_history_X1             pulse clinic_site_South             sex_M 
           0.0015            0.0012            0.0010           -0.0018

2. sbp is largest (0.0135). That tells us the forest relied on sbp; it does not tell us that lowering a patient’s SBP would lower their chance of Yes. 3. Predictors that carry similar information split the credit between them, and a column can be useful only because it travels with something we did not measure. Even in this fabricated data, the numbers describe the fitted model, not what happens if someone changes a measurement.

Explain why it would be wrong to pick min_n based on the final test score.

Show solution

Answer: The test set is supposed to give an honest score for a model that was chosen without it. If you choose min_n by test score, the test set takes part in the choice, so the winning score is partly luck that was selected for, and it is too optimistic for new patients.

Step by step

The tuning table shows how we choose: using training folds only.

collect_metrics(tuning) |>
  filter(.metric == "roc_auc") |>
  select(min_n, mean, std_err)

Output

# A tibble: 3 x 3
  min_n  mean std_err
  <int> <dbl>   <dbl>
1     2 0.795  0.0770
2     5 0.793  0.0756
3    10 0.804  0.0816

2. The three means differ by about 0.01 while their standard errors are about 0.08, so which min_n “wins” depends on the sample. The test set has only 120 rows (9 positives), so it is even noisier. 3. If we tried all three on the test set and kept the best, we would be reporting the largest of several noisy numbers. Use the folds to choose min_n (here 10), then score the chosen workflow once on the test set.

Name a setting where PCA or clustering is useful but does not provide diagnosis labels.

Show solution

Answer: Exploring a set of correlated vital signs. PCA can reduce age, BMI, SBP, and pulse to a couple of summary scores for plotting or checking which measurements move together, without using any outcome label. Another example is clustering clinic visits by their pattern of measurements to look for groups for further review.

Step by step

Run PCA on the four measurements, with no outcome column.

vitals <- ehr |> select(age, bmi, sbp, pulse) |> tidyr::drop_na()
pca <- prcomp(vitals, scale. = TRUE)
round(summary(pca)$importance[2:3, ], 3)

Output

                         PC1   PC2   PC3  PC4
Proportion of Variance 0.458 0.258 0.244 0.04
Cumulative Proportion  0.458 0.716 0.960 1.00

2. The first component holds 45.8% of the variation and the second 25.8%, so two scores summarize about 72% of how these four columns vary. 3. Nothing in this analysis says which patients have Yes. It finds structure in the measurements only. Any link to a diagnosis would need a labeled outcome and a validated, supervised analysis.

Checkpoint

You can describe trees, fit a random forest, tune an instructional parameter on folds, and explain why model comparison must precede the final test.