Lecture 13 Demo: Trees, Bagging, and Random Forests

Packages and data

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

library(ggplot2)
library(glmnet)
library(randomForest)
library(rpart)

We use the Hitters data from ISLR2, the same dataset used in the Lecture 10 R demo. The response is Salary, measured in thousands of dollars.

The preprocessing is similar to Lecture 10:

• remove rows with missing salary by using na.omit;
• keep the original variables;
• convert factor predictors to dummy variables with model.matrix when a numeric matrix is needed.

For the main tree and forest methods, we do not need to add all two-way interactions manually. Trees and forests can learn interactions through recursive splitting. Later, we also include ridge and Lasso with the Lecture 10 interaction-expanded design as a comparison.

data("Hitters", package = "ISLR2")
hitters <- na.omit(Hitters)

hitters

• Variables such as AtBat, Hits, HmRun, Runs, RBI, and Walks describe recent batting performance.
• Variables beginning with C are career totals.
• PutOuts, Assists, and Errors describe fielding.
• League, Division, and NewLeague are categorical indicators.

Train/test split

set.seed(32950)
n <- nrow(hitters)
train_id <- sample(seq_len(n), size = floor(0.8 * n))
test_id <- setdiff(seq_len(n), train_id)

train <- hitters[train_id, ]
test <- hitters[test_id, ]

length(train_id)

[1] 210

length(test_id)

[1] 53

For regularized linear models, we use a numeric design matrix. The factor variables are converted to dummy variables by model.matrix. The tree methods below use formula interfaces, so they can handle the factor variables directly.

X <- model.matrix(Salary ~ ., data = hitters)[, -1]
y <- hitters$Salary

X_train <- X[train_id, , drop = FALSE]
X_test <- X[test_id, , drop = FALSE]
y_train <- y[train_id]
y_test <- y[test_id]

dim(X_train)

[1] 210  19

Linear baselines for comparison

Before fitting trees, compare against linear methods (OLS, ridge and Lasso) on the same predictors.

lm_fit <- lm(Salary ~ ., data = train)
pred_lm <- predict(lm_fit, newdata = test)
# OLS fit test error
mse_lm <- mean((y_test - pred_lm)^2)

# Ridge fit
set.seed(32950)
cv_ridge <- cv.glmnet(X_train, y_train, alpha = 0, standardize = TRUE)

# Lasso fit
set.seed(32950)
cv_lasso <- cv.glmnet(X_train, y_train, alpha = 1, standardize = TRUE)

pred_ridge <- predict(cv_ridge, newx = X_test, s = "lambda.min")
pred_lasso <- predict(cv_lasso, newx = X_test, s = "lambda.min")

mse_ridge <- mean((y_test - pred_ridge)^2)
mse_lasso <- mean((y_test - pred_lasso)^2)

data.frame(
  Method = c("Linear regression", "Ridge", "Lasso"),
  Test_MSE = c(mse_lm, mse_ridge, mse_lasso)
)

For comparison with the Lecture 10 demo, we also fit ridge and Lasso after adding two-way interactions among the quantitative predictors.

reg_formula_int <- Salary ~ League + Division + NewLeague +
  (AtBat + Hits + HmRun + Runs + RBI + Walks + Years +
   CAtBat + CHits + CHmRun + CRuns + CRBI + CWalks +
   PutOuts + Assists + Errors)^2

X_int <- model.matrix(reg_formula_int, data = hitters)[, -1]

X_int_train <- X_int[train_id, , drop = FALSE]
X_int_test <- X_int[test_id, , drop = FALSE]

dim(X_int_train)

[1] 210 139

set.seed(32950)
cv_ridge_int <- cv.glmnet(X_int_train, y_train, alpha = 0, standardize = TRUE)

set.seed(32950)
cv_lasso_int <- cv.glmnet(X_int_train, y_train, alpha = 1, standardize = TRUE)

pred_ridge_int <- predict(cv_ridge_int, newx = X_int_test, s = "lambda.min")
pred_lasso_int <- predict(cv_lasso_int, newx = X_int_test, s = "lambda.min")

mse_ridge_int <- mean((y_test - pred_ridge_int)^2)
mse_lasso_int <- mean((y_test - pred_lasso_int)^2)

data.frame(
  Method = c("Ridge with two-way interactions", "Lasso with two-way interactions"),
  Test_MSE = c(mse_ridge_int, mse_lasso_int)
)

One regression tree

We first fit one large regression tree and then use cost-complexity pruning, as in Lecture 12. The CV-selected tree is used for prediction. Separately, we make a smaller tree for visualization so that the splits are readable in class.

set.seed(32950)
tree_big <- rpart(
  Salary ~ .,
  data = train,
  method = "anova",
  control = rpart.control(cp = 0, minsplit = 10, xval = 10) # minsplit: the minimum number of observations in a node for splitting; xval: number of CVs
)

In Lecture 12, cost-complexity pruning chooses a subtree by minimizing

\[ \sum_{m=1}^{|T|} \sum_{i:x_i \in R_m} (y_i - \bar y_m)^2 + \alpha |T|. \]

The rpart output uses the same idea, but it calls the tuning parameter CP instead of \(\alpha\). The CP values are scaled relative to the root-node training error, so they are not on the raw RSS scale. This scaling changes the units, not the logic: larger CP means a stronger penalty and a smaller tree.

The pruning table is the computed cost-complexity path. Each row corresponds to one candidate subtree:

• nsplit: number of internal splits in the tree;
• CP: rpart’s scaled version of the pruning penalty \(\alpha\);
• rel_error: training error relative to the root-only tree;
• xerror: cross-validation error relative to the root-only tree;
• xstd: standard error of the cross-validation error.

So rpart is not using a different pruning idea. It is implementing the weakest-link cost-complexity pruning path, then using cross-validation (xval = 10) to estimate prediction error for each subtree.

Using CV error for pruning the tree

cp_table <- as.data.frame(tree_big$cptable)
cp_table$nsplit <- as.integer(cp_table$nsplit)
names(cp_table) <- c("alpha_scaled_CP", "splits", "train_relative_error",
                     "cv_relative_error", "cv_standard_error")
cp_table

plotcp(tree_big)

In this example, the CV error does not rise much as the number of splits increases. This does not mean that larger trees are always better. It means that, for this particular training split and 10-fold CV run, CV does not find strong evidence that the smaller trees predict better than some larger trees. The xstd column is useful here: if several subtrees have CV errors within roughly one standard error of the minimum, their predictive performance is not clearly distinguishable.

In practice, a common alternative is the one-standard-error rule: choose the smallest tree whose CV error is within one standard error of the minimum. That gives a simpler tree, but it adds one more selection rule. To keep this demo focused, we use the minimum-CV tree for test MSE and use a separate small tree only for visualization.

The prediction tree uses the minimum CV error over the full cost-complexity path.

best_cp_cv <- tree_big$cptable[which.min(tree_big$cptable[, "xerror"]), "CP"]
tree_cv <- prune(tree_big, cp = best_cp_cv)

tree_cv_summary <- data.frame(
  alpha_scaled_CP = best_cp_cv,
  Splits = sum(tree_cv$frame$var != "<leaf>"),
  Leaves = sum(tree_cv$frame$var == "<leaf>"),
  CV_Relative_Error = min(tree_big$cptable[, "xerror"]),
  Test_MSE = NA_real_
)


pred_tree_cv <- predict(tree_cv, newdata = test)
mse_tree_cv <- mean((y_test - pred_tree_cv)^2)

tree_cv_summary$Test_MSE <- mse_tree_cv
tree_cv_summary

The single pruned tree has noticeably larger test MSE than the linear baselines here. This is not a mistake: a single tree is easy to interpret, but its piecewise-constant predictions can be unstable and relatively crude. The ensemble methods below are designed to improve this weakness by averaging many trees.

Visualization of a smaller tree

For visualization only, choose the best CV tree among cost-complexity subtrees with at most five splits. This is easier to read, but it is not the tree used for the test MSE comparison.

cptable_display <- tree_big$cptable[tree_big$cptable[, "nsplit"] <= 5, , drop = FALSE]
best_cp_display <- cptable_display[which.min(cptable_display[, "xerror"]), "CP"]
tree_display <- prune(tree_big, cp = best_cp_display)

tree_display_summary <- data.frame(
  alpha_scaled_CP = best_cp_display,
  Splits = sum(tree_display$frame$var != "<leaf>"),
  Leaves = sum(tree_display$frame$var == "<leaf>"),
  CV_Relative_Error = min(cptable_display[, "xerror"])
)

pred_tree_display <- predict(tree_display, newdata = test)
mse_tree_display <- mean((y_test - pred_tree_display)^2)

tree_display_summary$Test_MSE <- mse_tree_display
tree_display_summary

NA

You can see that the display tree here actually have smaller MSE on the test data than the CV selected tree. Actually, that is not suprising, as the CV error are very similar and a smaller tree is less likely to overfit.

The plot below is the small display tree. In each terminal node, the number shown is the predicted salary for observations that fall into that leaf.

rpart.plot::rpart.plot(
    tree_display,
    type = 2,
    extra = 101,
    fallen.leaves = TRUE,
    box.palette = "Blues",
    main = "Small display tree for salary"
)

The same tree can be read as a small set of terminal-node rules. Each row is one leaf:

• N: number of training observations in that leaf;
• Prediction: average training salary in that leaf;
• Rule: conditions needed to reach that leaf.

leaf_nodes <- as.numeric(row.names(tree_display$frame[tree_display$frame$var == "<leaf>", ]))
tree_paths <- path.rpart(tree_display, nodes = leaf_nodes, print.it = FALSE)

tree_rules <- data.frame(
  Leaf = leaf_nodes,
  N = tree_display$frame[as.character(leaf_nodes), "n"],
  Prediction = round(tree_display$frame[as.character(leaf_nodes), "yval"], 1),
  Rule = sapply(tree_paths, function(z) paste(z[-1], collapse = " & ")),
  row.names = NULL
)

tree_rules

A single tree is indeed very interpretative, though its prediction strength is not ideal.

Bagging and random forest

We now use the randomForest package to perform bagging and random forest.

• Bagging: bootstrap trees with mtry = p, so each split may consider all predictors.
• Random forest: bootstrap trees with mtry < p, so each split considers only a random subset of predictors.

The package reports an out-of-bag (OOB) MSE after each number of trees. OOB error is computed from observations that were not included in a tree’s bootstrap sample.

Bagged trees

Bagging grows many large trees on bootstrap samples and averages their predictions. Here mtry = p, so each split can use any predictor.

set.seed(32950)
p_tree <- ncol(train) - 1

bag_fit <- randomForest(
  Salary ~ .,
  data = train,
  ntree = 500,
  mtry = p_tree,
  importance = TRUE
)

bag_fit

Call:
 randomForest(formula = Salary ~ ., data = train, ntree = 500,      mtry = p_tree, importance = TRUE) 
               Type of random forest: regression
                     Number of trees: 500
No. of variables tried at each split: 19

          Mean of squared residuals: 82670.99
                    % Var explained: 58.45

bag_oob_mse <- bag_fit$mse[bag_fit$ntree]
bag_oob_mse

[1] 82670.99

We can visualize how OOB MSE change with the number of trees.

 plot(bag_fit$mse, type = "l",
       xlab = "Number of trees",
       ylab = "OOB MSE",
       main = "Bagging: OOB error vs number of trees")

We see here that OOB MSE already stabilizes at \(B = 500\). Increasing the number of trees in bagging or random forests typically does not overfit. We typically start with a large \(B\) and increase \(B\) only if OOB MSE does not stabilize.

pred_bag <- predict(bag_fit, newdata = test)
mse_bag <- mean((y_test - pred_bag)^2)
mse_bag

[1] 71601.36

The prediction test error is now smaller than linear methods, indicating that considering interactions of predictors help. Notice that in the R demo of lecture 10, the best test data MSE we get is 60647 from Lasso with interaction terms, consistent with the results here.

Random forest

For regression, a common default is \(m_{\mathrm{try}} \approx p/3\). This keeps each tree strong while reducing correlation across trees.

set.seed(32950)
mtry_rf <- max(1, floor(p_tree / 3))
mtry_rf

[1] 6

set.seed(32950)
rf_fit <- randomForest(
  Salary ~ .,
  data = train,
  ntree = 500,
  mtry = mtry_rf,
  importance = TRUE
)

rf_fit

Call:
 randomForest(formula = Salary ~ ., data = train, ntree = 500,      mtry = mtry_rf, importance = TRUE) 
               Type of random forest: regression
                     Number of trees: 500
No. of variables tried at each split: 6

          Mean of squared residuals: 80668.42
                    % Var explained: 59.46

rf_oob_mse <- rf_fit$mse[rf_fit$ntree]
rf_oob_mse

[1] 80668.42

plot(rf_fit$mse, type = "l",
       xlab = "Number of trees",
       ylab = "OOB MSE",
       main = "Random forest: OOB error vs number of trees")

pred_rf <- predict(rf_fit, newdata = test)
mse_rf <- mean((y_test - pred_rf)^2)
mse_rf

[1] 68668.24

We further drop the OOB MSE and test data MSE using random forest.

Variable importance for the random forest

Permutation importance measures how much prediction error increases after breaking the relationship between one predictor and the response. The randomForest package computes this using OOB observations.

importance_df <- data.frame(
  Variable = row.names(importance(rf_fit, type = 1)),
  Increase_in_MSE = importance(rf_fit, type = 1)[, 1],
  row.names = NULL
)

importance_top <- head(importance_df[order(-importance_df$Increase_in_MSE), ], 10)
importance_top

Compare test-set performance across methods

results <- data.frame(
  Method = c("Linear regression", "Ridge", "Lasso",
             "Ridge with two-way interactions", "Lasso with two-way interactions",
             "CV-pruned tree", "Bagged trees", "Random forest"),
  Test_MSE = c(mse_lm, mse_ridge, mse_lasso,
               mse_ridge_int, mse_lasso_int,
               mse_tree_cv, mse_bag, mse_rf)
)

results[order(results$Test_MSE), ]

pred_df <- data.frame(
  Salary = rep(y_test, 8),
  Prediction = c(pred_lm, as.numeric(pred_ridge), as.numeric(pred_lasso),
                 as.numeric(pred_ridge_int), as.numeric(pred_lasso_int),
                 pred_tree_cv, pred_bag, pred_rf),
  Method = rep(c("Linear regression", "Ridge", "Lasso",
                 "Ridge with interactions", "Lasso with interactions",
                 "CV-pruned tree", "Bagged trees", "Random forest"),
               each = length(y_test))
)

ggplot(pred_df, aes(x = Salary, y = Prediction, color = Method)) +
  geom_point(alpha = 0.7) +
  geom_abline(slope = 1, intercept = 0, linetype = 2) +
  facet_wrap(~ Method) +
  labs(
    title = "Observed versus predicted salary on the test set",
    x = "Observed salary",
    y = "Predicted salary"
  ) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "none")

In this test split, the two highest-salary observations are fitted better by the linear methods than by the tree-based methods. One reason is that these salaries are much larger than most other observations. Regression trees, bagging, and random forests make predictions by averaging training responses within leaves or across trees, so predictions for extreme observations are often pulled back toward more typical salaries. Linear methods can produce more extreme fitted values through linear combinations of predictors, especially when interaction terms are included. Because MSE squares errors, underpredicting a few very high-salary players can noticeably increase the test MSE for tree-based methods, which may partly explain why the interaction Lasso performs well here.

For many lower-salary observations, the scatter plot suggests that the tree-based methods predict better than the linear methods. This is consistent with the averaging nature of trees: in the denser lower-salary range, local averaging can stabilize predictions and reduce large deviations. Thus, the same averaging behavior can help tree-based methods for typical observations while hurting them for extreme high-salary observations.

Takeaway

• The same Hitters preprocessing from Lecture 10 works here: remove missing values and encode factors when a numeric matrix is needed.
• Linear methods give useful baselines.
• A small pruned tree is interpretable but can have high prediction variance.
• Bagging averages many bootstrap trees and usually improves test error.
• Random forests add random predictor selection at each split, reducing correlation among trees.
• Variable importance summarizes a forest, but it is predictive rather than causal.