Lecture 3: Tree-Based Models

Machine Learning for Macro Forecasting and Financial Econometrics (SIdE 2026)

Author

Rodrigo Sarlo

Published

July 1, 2026

Overview

Welcome to the third practical session. Today, we transition from linear and regularized linear models to non-linear, tree-based machine learning methods.

In Part 1, we explore single Decision Trees. We discuss how they capture non-linear relationships and interactive effects without requiring explicit functional form specifications. Additionally, we conclude our introduction to the {tidymodels} framework, demonstrating how to select optimal hyperparameters from tuning objects generated by tune_grid(), visualize tuning diagnostics, and extract predictions for analysis. Looking ahead, we evaluate all tuned models on our holdout test sample in the next class.

In Part 2, we transition to tree ensembles: Random Forests and Gradient Boosted Trees. For boosting, we focus on the widely-used XGBoost implementation, while briefly discussing alternative variants. We examine the bias-variance trade-offs of bagging and boosting, show how to configure the most well-known hyperparameters, and extract Variable Importance Metrics (VIM) using the {vip} package to interpret these models.

Part 1: Decision Trees & Model Selection in Tidymodels

1.1 Mathematical Background

A Decision Tree partitions the predictor space $\mathbb{R}^p$ into a set of $M$ non-overlapping, multi-dimensional rectangles (or regions) $R_1, R_2, \dots, R_M$. For any input $x \in R_m$, the model predicts a constant value $c_m$, which is typically the mean of the target variable $y$ within that region.

The model can be written as:

\[f(x) = \sum_{m=1}^{M} c_m I(x \in R_m)\]

where $I(\cdot)$ represents the indicator function that returns $1$ if $x \in R_m$ and $0$ otherwise.

Recursive Binary Splitting

To construct these regions, we use a greedy, top-down approach known as recursive binary splitting. Beginning at the root node containing all observations, we select the splitting predictor $X_j$ and the split point $s$ that partition the space into $R_1 = \{X | X_j \le s\}$ and $R_2 = \{X | X_j > s\}$ to minimize the sum of squared residuals (or another target loss function):

\[\min_{j, s} \left\{ \sum_{i: x_{ij} \le s} (y_i - \hat{y}_{R_1})^2 + \sum_{i: x_{ij} > s} (y_i - \hat{y}_{R_2})^2 \right\}\]

where $\hat{y}_{R_1}$ is the mean of $y$ in region $R_1$, and $\hat{y}_{R_2}$ is the mean of $y$ in region $R_2$. This splitting process is repeated recursively on the resulting child nodes.

Computational Complexity & Local Structure

While recursive binary splitting is conceptually straightforward, searching for the optimal split $(j, s)$ at every node is computationally intensive, scaling with both the number of predictors $p$ and sample size $T$. This computational burden increases dramatically for models that implement more complex terminal-node structures than a simple fixed constant (the local mean).

For example, the Generalized Random Forest (GRF) framework of Athey, Tibshirani, and Wager (2019) estimates non-parametric models (such as local linear regressions or causal treatment effects) at each leaf. Evaluating splits based on these localized statistical models rather than simple averages requires solving complex optimization problems at every candidate partition, significantly increasing computational time.

Cost-Complexity Pruning

Growing a tree until every terminal node (leaf) has very few observations leads to severe overfitting - a well-known limitation of single decision trees that the ensemble methods discussed in Part 2 seek to address. To mitigate this, we can grow a large tree $T_0$ and prune it back using cost-complexity pruning by minimizing:

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

where $|T|$ is the number of terminal nodes in tree $T$, and $\alpha \ge 0$ is the tuning parameter (called cost_complexity in {tidymodels}) that penalizes tree size.

1.2 Model Specification & Hyperparameter Tuning

In this section, we construct and run our model-tuning pipeline. We first specify a regression decision tree using {parsnip} and flag three key parameters for optimization: cost_complexity, tree_depth, and min_n. We load the pre-configured preprocessing recipe (fredmd_recipe) and rolling validation splits (folds) from Lecture 2, and bundle these components into a unified workflow. Finally, we execute the grid search to generate our downstream performance diagnostics.

Note that there are two main differences between the setup here and the regularized linear models from the previous lecture: (i) model-specific hyperparameters must be optimized; and (ii) since {rpart} handles both regression and classification problems, we must explicitly set the model mode using set_mode().

library(tidyverse)
library(tidymodels)

# Specify the decision tree model for regression
tree_spec <- decision_tree(
  cost_complexity = tune(),
  tree_depth = tune(),
  min_n = tune()
) %>%
  set_engine("rpart") %>%
  # Setting the mode is required when a package can work with different ML problems
  # Options available are:  "censored regression", "classification", "quantile regression", and "regression"
  set_mode("regression")

# Load the recycled objects
load("data/recycled_objects.RData")
treated_data <- read_csv("data/train_val_fredmd.csv")

# Verify that fredmd_recipe has been loaded into the workspace
print(fredmd_recipe)

# Create the workflow
tree_workflow <- workflow() %>%
  add_recipe(fredmd_recipe) %>%
  add_model(tree_spec)

# Create a regular grid of hyperparameters to evaluate
tree_grid <- grid_regular(
  cost_complexity(range = c(0, 0.01), trans = NULL),
  tree_depth(range = c(3, 10)),
  min_n(range = c(5, 20)),
  levels = 3
)
print(tree_grid)


# Optional: Parallel tuning backend for efficiency
# library(future)
# plan(multisession, workers = parallel::detectCores() - 4)

# Execute grid tuning over the rolling splits
tree_results <- tune_grid(
  tree_workflow,
  resamples = folds,
  grid = tree_grid,
  metrics = metric_set(rmse, mae),
  control = control_grid(save_pred = TRUE, verbose = TRUE)
)

# plan(sequential) # Reset if using future parallel execution

# Saving this tuned model for our final comparison next class
save(tree_results, tree_spec, tree_workflow, file = "data/tuned_trees.RData")

1.3 Model Selection & Diagnostics

Once the grid search is complete, we evaluate the performance of our hyperparameter combinations across the validation splits, select the optimal configurations, and finalize our modeling workflow.

Performance evaluation in {tidymodels} is managed by the {yardstick} package. While we rely on metric_set(rmse, mae), {yardstick} supports many regression metrics, including rsq (R-squared), mape (absolute percentage error), and robust measures like huber_loss.

To inspect results, {tune} provides an implementation of the autoplot() method, which automatically constructs {ggplot2} plots based on the input object’s class. Beyond the default view, we can filter metrics or isolate top-performing parameter combinations:

# 1. Default visualization: plot all metrics across hyperparameters
autoplot(tree_results)

# 2. Isolate a specific metric (e.g., RMSE)
autoplot(tree_results, metric = "rmse")

# 3. Focus visual attention on the top-performing grid combinations
# In the current case, we get the same results from #1
autoplot(tree_results, select_best = TRUE)

To identify the optimal configuration from our grid, {tune} provides specialized selection functions:

select_best(): Returns the single parameter set that minimizes error (or maximizes performance) on the metric of interest.
Parsimonious Selection: Functions like select_by_one_std_err() and select_by_pct_loss() choose simpler models whose performance is within a certain threshold of the best model (e.g., within one standard error or a specified percentage loss) to mitigate overfitting.

We first display the top configurations and then select the absolute best parameters to minimize RMSE:

# Display the top 5 hyperparameter configurations ranked by RMSE
show_best(tree_results, metric = "rmse", n = 5)

# Select the absolute best parameter set (minimizing RMSE)
best_tree_params <- select_best(tree_results, metric = "rmse")
print(best_tree_params)

The finalize_workflow() function is a critical step in the {tidymodels} lifecycle: it updates the model specification (and optionally recipe steps) by replacing the tune() placeholders with the concrete parameter values and an object ready to be fitted. This creates a fully specified, self-contained pipeline object.

# Finalize the workflow
final_tree_wf <- finalize_workflow(tree_workflow, best_tree_params)
print(final_tree_wf)

We fit the finalized workflow on the train/validation dataset and plot the tree structure using {rpart.plot}.

library(rpart.plot)

# Fit model to the train/validation dataset
final_tree_fit <- fit(final_tree_wf, data = tail(treated_data, 360))
print(final_tree_fit)

# Extract and plot the estimated decision tree
tree_obj <- extract_fit_parsnip(final_tree_fit)$fit
rpart.plot(tree_obj)

1.4 Forecasting

We demonstrate how to retrieve the forecasts generated during the rolling validation process, and how to use our finalized model to forecast on the holdout test sample.

1: Predictions from the Validation Period
2: Holdout Forecasts

Since we specified save_pred = TRUE in our tuning controls, we can extract the out-of-sample predictions generated on the assessment folds during validation. We use collect_predictions() and filter by our optimal parameter values using the parameters argument.

# Extract validation predictions for the best model configuration
val_predictions <- collect_predictions(tree_results, parameters = best_tree_params)
print(head(val_predictions))

To generate out-of-sample forecasts, we load the holdout test set (test_fredmd.csv) and apply the standard predict() method to our finalized model fit.

Because the final_tree_fit object contains the entire preprocessing pipeline, {tidymodels} automatically handles the transformations on this new dataset before generating predictions.

# Load the holdout test dataset
holdout_data <- read_csv("data/test_fredmd.csv")

# Predict on the holdout test data using the finalized full fit
holdout_predictions <- predict(final_tree_fit, new_data = holdout_data) %>%
  bind_cols(holdout_data %>% select(date, CPIAUCSL))
print(head(holdout_predictions))

Part 2: Random Forests & Gradient Boosting

2.1 Mathematical Background

A key limitation of single decision trees is their instability: they are highly sensitive to small perturbations in the training data. Because splits are chosen greedily at each node, a minor change in a single observation can result in a completely different initial split, altering the structure of the entire downstream tree.

In terms of predictive performance, this instability translates directly into high estimator variance, which leads to a high out-of-sample Mean Squared Error (MSE). Ensemble methods address this shortfall by combining multiple trees to stabilize predictions and minimize MSE:

Bagging & Random Forests (Parallel Learning): Reduce variance by averaging the outputs of multiple deep, fully grown trees (which are individually highly unstable, but stable when averaged).
Boosting (Sequential Learning): Focuses on sequentially correcting prediction errors by growing shallow, low-variance trees (which are highly stable due to their limited depth) to fit the negative gradient of the loss function with respect to the preceding model’s predictions. Under squared-error loss, this is mathematically equivalent to sequentially fitting the residuals of the current ensemble.

Bootstrap Aggregation (Bagging) & Random Forests

Bagging (bootstrap aggregation) averages predictions across $B$ bootstrapped training sets to decrease prediction variance. In standard bagging and random forests, these individual trees are fully grown to keep bias as low as possible, and their resulting high variance is mitigated by averaging the predictions:

\[\hat{f}_{bag}(x) = \frac{1}{B} \sum_{b=1}^{B} \hat{f}^{*b}(x)\]

Random Forests improve upon bagging by decorrelation. If a dataset contains a single highly dominant predictor, the trees grown on different bootstrap samples will still look very similar, resulting in highly correlated predictions and limiting the variance reduction achieved by averaging. To prevent this, the Random Forest algorithm restricts each split to a random subset of $m \approx \sqrt{p}$ (or $p/3$ for regression) predictors. This decorrelates the individual trees, allowing averaging to reduce ensemble variance much more effectively.

Gradient Boosted Trees

Unlike Random Forests, which grow independent trees in parallel, Boosting grows trees sequentially. The pseudo-code below describes the classical gradient boosting algorithm for regression. Note that the sequential fitting of the residuals ($r_i$) implies that we are assuming a squared-error (MSE) loss function:

Set $\hat{f}(x) = \bar{y}$ and $r_i = y_i-\bar{y}$ for all $i$.
For $b = 1, 2, \dots, B$:
1. Fit a tree $\hat{f}^b$ with $d$ splits to the training data $(X, r)$.
2. Update the ensemble: $\hat{f}(x) \leftarrow \hat{f}(x) + \eta \hat{f}^b(x)$, where $\eta \in (0,1]$ is the learning rate (shrinkage).
3. Update the residuals: $r_i \leftarrow r_i - \eta \hat{f}^b(x_i)$.
The final boosted model is: $\hat{f}(x) = \sum_{b=1}^{B} \eta \hat{f}^b(x)$.

Generality of Ensemble Strategies

While bagging and boosting are most famous for their implementations with decision trees, they are general meta-algorithms that can be applied to any base learner (such as neural networks, support vector machines, or linear models). For example, the {xgboost} implementation is not strictly limited to decision trees; it supports alternative base learners (e.g., regularized linear regressions) as the models fitted in step 2a.

2.2 Random Forest Specification & Tuning

We specify a Random Forest using rand_forest(), setting the engine to ranger and importance = "impurity". This setting is required to calculate the importance scores needed for our downstream {vip} plots in Section 2.4, which default to "none".

Besides the engine’s name, and the model passed through the pipe operator %>%, extra named arguments passed to set_engine() are sent directly to the selected engine.

In this setup, we optimize model performance by focusing on min_n and keeping the remaining parameters fixed:

min_n: The minimum observation threshold to attempt a split, which we tune to control tree complexity.
mtry: The number of randomly sampled predictors at each split, which we fix at floor(p/3) (164 candidate predictors). As noted in the mathematical background, this is a standard rule of thumb for regression trees that saves substantial computational time by avoiding tuning an extra hyperparameter.
trees: The number of trees grown in the ensemble, which we fix at 500. Unlike parameters that control individual tree size, increasing the number of trees does not lead to overfitting. Once the ensemble is large enough to stabilize the variance reduction, adding more trees only raises computational runtime without improving predictive performance.

# Specify Random Forest model
rf_spec <- rand_forest(
  mtry = floor(493/3),
  trees = 500,
  min_n = tune()
) %>%
  set_engine("ranger", importance = "impurity") %>%
  set_mode("regression")

rf_workflow <- workflow() %>%
  add_recipe(fredmd_recipe) %>%
  add_model(rf_spec)

# Define tuning grid
rf_grid <- grid_regular(
  min_n(range = c(5, 50)),
  levels = 5
)

# library(future)
# plan(multisession, workers = parallel::detectCores() - 4)

set.seed(456)

# Tune Random Forest over the folds
rf_results <- tune_grid(
  rf_workflow,
  resamples = folds,
  grid = rf_grid,
  metrics = metric_set(rmse, mae),
  control = control_grid(save_pred = TRUE, verbose = TRUE)
)

# plan(sequential)

# Saving the tuned model
save(rf_results, rf_spec, rf_workflow, file = "data/tuned_rf.RData")

# Quick inspection of the hyperparameter tuning process
autoplot(rf_results)

# Select best parameters
best_rf_params <- select_best(rf_results, metric = "rmse")

# Finalize and fit model
final_rf_wf <- finalize_workflow(rf_workflow, best_rf_params)
set.seed(456)
final_rf_fit <- fit(final_rf_wf, data = tail(treated_data, 360))

2.3 Gradient Boosting

Unlike Random Forests, Gradient Boosted Trees (often simply referred to as Gradient Boosting) are highly susceptible to overfitting if the number of sequential boosting rounds, $B$, is too large. To mitigate this risk, we optimize two key hyperparameters: the learning rate (learn_rate) and the tree depth (tree_depth). Additionally, a widely used strategy in practice is to implement early stopping using a held-out validation sample. Early stopping monitors the model’s performance on this validation set at each boosting iteration and halts training if the validation loss fails to improve for a specified number of consecutive rounds. This combination of highly non-linear base learners and sequential updates makes Gradient Boosting extremely powerful.

Although highly effective, we omit early stopping here to avoid look-ahead bias. In {tidymodels}, the default early-stopping wrapper partitions the validation set randomly rather than chronologically. While we could bypass this wrapper to force a chronological split, doing so is difficult to integrate with the unified resampling pipeline. This highlights a key trade-off of high-level frameworks: we trade engine-specific customization for pipeline simplicity. In what follows, we control overfitting solely by tuning the learning rate, tree depth, and boosting rounds.

# Specify XGBoost model
xgb_spec <- boost_tree(
  trees = tune(),
  tree_depth = tune(),
  learn_rate = tune()
) %>%
  set_engine("xgboost") %>%
  set_mode("regression")

xgb_workflow <- workflow() %>%
  add_recipe(fredmd_recipe) %>%
  add_model(xgb_spec)

# Define tuning grid
xgb_grid <- grid_regular(
  trees(range = c(500, 2000)),
  tree_depth(range = c(1, 6)),
  learn_rate(range = c(0.005, 0.1), trans = NULL),
  levels = 3
)

# plan(multisession, workers = parallel::detectCores() - 4)

set.seed(789)

# Tune XGBoost over the folds
xgb_results <- tune_grid(
  xgb_workflow,
  resamples = folds,
  grid = xgb_grid,
  metrics = metric_set(rmse, mae),
  control = control_grid(save_pred = TRUE, verbose = TRUE)
)

# plan(sequential)

# Saving the tuned model
save(xgb_results, xgb_spec, xgb_workflow, file = "data/tuned_xgb.RData")

# Quick inspection of the hyperparameter tuning process
autoplot(xgb_results)

# Select best parameters
best_xgb_params <- select_best(xgb_results, metric = "rmse")

# Finalize and fit model
final_xgb_wf <- finalize_workflow(xgb_workflow, best_xgb_params)
set.seed(789)
final_xgb_fit <- fit(final_xgb_wf, data = tail(treated_data, 360))

2.4 Variable Importance

Some modeling approaches, such as linear regressions or decision trees, are inherently interpretable as we can easily grasp what the estimated coefficients, representing partial effects, or sequential splitting rules represent for the relationships being studied. This transparency is lost, however, when using more general non-linear models - including tree ensembles like Random Forests and Gradient Boosting that aggregate hundreds of individual trees. To understand what happens inside these estimators, we must use ad-hoc post-estimation methods. Variable Importance (VI) solves this by measuring how much each predictor contributes to reducing the model’s prediction error.

Specifically, tree ensembles typically compute variable importance using either impurity reduction (i.e., reduction in node variance) or gain (measuring improvement in training loss). While these metrics are specific to trees, more general alternatives, such as permutation-based importance, can be used to evaluate feature contribution across any type of machine learning model.

Here we show how to use the {vip} package to extract and plot the top 20 predictors for our finalized models.

Variable Importance Plots

library(vip)
library(gridExtra)

# Plot Random Forest VIM
p_rf <- final_rf_fit %>%
  extract_fit_parsnip() %>%
  vip(
    num_features = 20, 
    geom = "col", 
    aesthetics = list(fill = "darkred")
  ) +
  theme_minimal() +
  labs(
    title = "Random Forest",
    subtitle = "Top 20 predictors (Impurity)",
    y = "Importance",
    x = "Predictor"
  )

# Plot XGBoost VIM
p_xgb <- final_xgb_fit %>%
  extract_fit_parsnip() %>%
  vip(
    num_features = 20, 
    geom = "col", 
    aesthetics = list(fill = "darkorange")
  ) +
  theme_minimal() +
  labs(
    title = "XGBoost",
    subtitle = "Top 20 predictors (Gain)",
    y = "Importance",
    x = "Predictor"
  )

grid.arrange(p_rf, p_xgb, ncol = 2)

--- title: "Lecture 3: Tree-Based Models" subtitle: "Machine Learning for Macro Forecasting and Financial Econometrics (SIdE 2026)" author: "Rodrigo Sarlo" date: "July 2026" format: html: toc: true toc-float: true code-fold: true code-summary: "Show the Code" code-tools: true --- ## Overview Welcome to the third practical session. Today, we transition from linear and regularized linear models to non-linear, tree-based machine learning methods. In **Part 1**, we explore single Decision Trees. We discuss how they capture non-linear relationships and interactive effects without requiring explicit functional form specifications. Additionally, we conclude our introduction to the `{tidymodels}` framework, demonstrating how to select optimal hyperparameters from tuning objects generated by `tune_grid()`, visualize tuning diagnostics, and extract predictions for analysis. Looking ahead, we evaluate all tuned models on our holdout test sample in the next class. In **Part 2**, we transition to tree ensembles: Random Forests and Gradient Boosted Trees. For boosting, we focus on the widely-used XGBoost implementation, while briefly discussing alternative variants. We examine the bias-variance trade-offs of bagging and boosting, show how to configure the most well-known hyperparameters, and extract Variable Importance Metrics (VIM) using the `{vip}` package to interpret these models. --- ## Part 1: Decision Trees & Model Selection in Tidymodels ### 1.1 Mathematical Background A Decision Tree partitions the predictor space $\mathbb{R}^p$ into a set of $M$ non-overlapping, multi-dimensional rectangles (or regions) $R_1, R_2, \dots, R_M$. For any input $x \in R_m$, the model predicts a constant value $c_m$, which is typically the mean of the target variable $y$ within that region. The model can be written as: $$f(x) = \sum_{m=1}^{M} c_m I(x \in R_m)$$ where $I(\cdot)$ represents the indicator function that returns $1$ if $x \in R_m$ and $0$ otherwise. #### Recursive Binary Splitting To construct these regions, we use a greedy, top-down approach known as recursive binary splitting. Beginning at the root node containing all observations, we select the splitting predictor $X_j$ and the split point $s$ that partition the space into $R_1 = \{X | X_j \le s\}$ and $R_2 = \{X | X_j > s\}$ to minimize the sum of squared residuals (or another target loss function): $$\min_{j, s} \left\{ \sum_{i: x_{ij} \le s} (y_i - \hat{y}_{R_1})^2 + \sum_{i: x_{ij} > s} (y_i - \hat{y}_{R_2})^2 \right\}$$ where $\hat{y}_{R_1}$ is the mean of $y$ in region $R_1$, and $\hat{y}_{R_2}$ is the mean of $y$ in region $R_2$. This splitting process is repeated recursively on the resulting child nodes. ##### Computational Complexity & Local Structure While recursive binary splitting is conceptually straightforward, searching for the optimal split $(j, s)$ at every node is computationally intensive, scaling with both the number of predictors $p$ and sample size $T$. This computational burden increases dramatically for models that implement more complex terminal-node structures than a simple fixed constant (the local mean). > For example, the Generalized Random Forest (GRF) framework of Athey, Tibshirani, and Wager (2019) estimates non-parametric models (such as local linear regressions or causal treatment effects) at each leaf. Evaluating splits based on these localized statistical models rather than simple averages requires solving complex optimization problems at every candidate partition, significantly increasing computational time. #### Cost-Complexity Pruning Growing a tree until every terminal node (leaf) has very few observations leads to severe overfitting - a well-known limitation of single decision trees that the ensemble methods discussed in Part 2 seek to address. To mitigate this, we can grow a large tree $T_0$ and prune it back using cost-complexity pruning by minimizing: $$\sum_{m=1}^{|T|} \sum_{i: x_i \in R_m} (y_i - c_m)^2 + \alpha |T|$$ where $|T|$ is the number of terminal nodes in tree $T$, and $\alpha \ge 0$ is the tuning parameter (called `cost_complexity` in `{tidymodels}`) that penalizes tree size. --- ### 1.2 Model Specification & Hyperparameter Tuning In this section, we construct and run our model-tuning pipeline. We first specify a regression decision tree using `{parsnip}` and flag three key parameters for optimization: `cost_complexity`, `tree_depth`, and `min_n`. We load the pre-configured preprocessing recipe (`fredmd_recipe`) and rolling validation splits (`folds`) from Lecture 2, and bundle these components into a unified workflow. Finally, we execute the grid search to generate our downstream performance diagnostics. Note that there are two main differences between the setup here and the regularized linear models from the previous lecture: (i) model-specific hyperparameters must be optimized; and (ii) since `{rpart}` handles both regression and classification problems, we must explicitly set the model mode using `set_mode()`. ::: {.panel-tabset} ##### 1: Specification ```r library(tidyverse) library(tidymodels) # Specify the decision tree model for regression tree_spec <- decision_tree( cost_complexity = tune(), tree_depth = tune(), min_n = tune() ) %>% set_engine("rpart") %>% # Setting the mode is required when a package can work with different ML problems # Options available are: "censored regression", "classification", "quantile regression", and "regression" set_mode("regression") ``` ##### 2: Workflow Bundle ```r # Load the recycled objects load("data/recycled_objects.RData") treated_data <- read_csv("data/train_val_fredmd.csv") # Verify that fredmd_recipe has been loaded into the workspace print(fredmd_recipe) # Create the workflow tree_workflow <- workflow() %>% add_recipe(fredmd_recipe) %>% add_model(tree_spec) ``` ##### 3: Grid Design ```r # Create a regular grid of hyperparameters to evaluate tree_grid <- grid_regular( cost_complexity(range = c(0, 0.01), trans = NULL), tree_depth(range = c(3, 10)), min_n(range = c(5, 20)), levels = 3 ) print(tree_grid) ``` ##### 4: Grid Execution ```r # Optional: Parallel tuning backend for efficiency # library(future) # plan(multisession, workers = parallel::detectCores() - 4) # Execute grid tuning over the rolling splits tree_results <- tune_grid( tree_workflow, resamples = folds, grid = tree_grid, metrics = metric_set(rmse, mae), control = control_grid(save_pred = TRUE, verbose = TRUE) ) # plan(sequential) # Reset if using future parallel execution # Saving this tuned model for our final comparison next class save(tree_results, tree_spec, tree_workflow, file = "data/tuned_trees.RData") ``` ::: --- ### 1.3 Model Selection & Diagnostics Once the grid search is complete, we evaluate the performance of our hyperparameter combinations across the validation splits, select the optimal configurations, and finalize our modeling workflow. ::: {.panel-tabset} ##### 1: Tuning Diagnostics Performance evaluation in `{tidymodels}` is managed by the `{yardstick}` package. While we rely on `metric_set(rmse, mae)`, `{yardstick}` supports many regression metrics, including `rsq` (R-squared), `mape` (absolute percentage error), and robust measures like `huber_loss`. To inspect results, `{tune}` provides an implementation of the `autoplot()` method, which automatically constructs `{ggplot2}` plots based on the input object's class. Beyond the default view, we can filter metrics or isolate top-performing parameter combinations: ```r # 1. Default visualization: plot all metrics across hyperparameters autoplot(tree_results) # 2. Isolate a specific metric (e.g., RMSE) autoplot(tree_results, metric = "rmse") # 3. Focus visual attention on the top-performing grid combinations # In the current case, we get the same results from #1 autoplot(tree_results, select_best = TRUE) ``` ##### 2: Select Best To identify the optimal configuration from our grid, `{tune}` provides specialized selection functions: * `select_best()`: Returns the single parameter set that minimizes error (or maximizes performance) on the metric of interest. * Parsimonious Selection: Functions like `select_by_one_std_err()` and `select_by_pct_loss()` choose simpler models whose performance is within a certain threshold of the best model (e.g., within one standard error or a specified percentage loss) to mitigate overfitting. We first display the top configurations and then select the absolute best parameters to minimize RMSE: ```r # Display the top 5 hyperparameter configurations ranked by RMSE show_best(tree_results, metric = "rmse", n = 5) # Select the absolute best parameter set (minimizing RMSE) best_tree_params <- select_best(tree_results, metric = "rmse") print(best_tree_params) ``` ##### 3: Finalize Workflow The `finalize_workflow()` function is a critical step in the `{tidymodels}` lifecycle: it updates the model specification (and optionally recipe steps) by replacing the `tune()` placeholders with the concrete parameter values and an object ready to be fitted. This creates a fully specified, self-contained pipeline object. ```r # Finalize the workflow final_tree_wf <- finalize_workflow(tree_workflow, best_tree_params) print(final_tree_wf) ``` ##### 4: Full Fit & Visualization We fit the finalized workflow on the train/validation dataset and plot the tree structure using `{rpart.plot}`. ```r library(rpart.plot) # Fit model to the train/validation dataset final_tree_fit <- fit(final_tree_wf, data = tail(treated_data, 360)) print(final_tree_fit) # Extract and plot the estimated decision tree tree_obj <- extract_fit_parsnip(final_tree_fit)$fit rpart.plot(tree_obj) ``` ::: --- ### 1.4 Forecasting We demonstrate how to retrieve the forecasts generated during the rolling validation process, and how to use our finalized model to forecast on the holdout test sample. ::: {.panel-tabset} ##### 1: Predictions from the Validation Period Since we specified `save_pred = TRUE` in our tuning controls, we can extract the out-of-sample predictions generated on the assessment folds during validation. We use `collect_predictions()` and filter by our optimal parameter values using the `parameters` argument. ```r # Extract validation predictions for the best model configuration val_predictions <- collect_predictions(tree_results, parameters = best_tree_params) print(head(val_predictions)) ``` ##### 2: Holdout Forecasts To generate out-of-sample forecasts, we load the holdout test set (`test_fredmd.csv`) and apply the standard `predict()` method to our finalized model fit. Because the `final_tree_fit` object contains the entire preprocessing pipeline, `{tidymodels}` automatically handles the transformations on this new dataset before generating predictions. ```r # Load the holdout test dataset holdout_data <- read_csv("data/test_fredmd.csv") # Predict on the holdout test data using the finalized full fit holdout_predictions <- predict(final_tree_fit, new_data = holdout_data) %>% bind_cols(holdout_data %>% select(date, CPIAUCSL)) print(head(holdout_predictions)) ``` ::: --- ## Part 2: Random Forests & Gradient Boosting ### 2.1 Mathematical Background A key limitation of single decision trees is their instability: they are highly sensitive to small perturbations in the training data. Because splits are chosen greedily at each node, a minor change in a single observation can result in a completely different initial split, altering the structure of the entire downstream tree. In terms of predictive performance, this instability translates directly into high estimator variance, which leads to a high out-of-sample Mean Squared Error (MSE). Ensemble methods address this shortfall by combining multiple trees to stabilize predictions and minimize MSE: 1. **Bagging & Random Forests** (Parallel Learning): Reduce variance by averaging the outputs of multiple deep, fully grown trees (which are individually highly unstable, but stable when averaged). 2. **Boosting** (Sequential Learning): Focuses on sequentially correcting prediction errors by growing shallow, low-variance trees (which are highly stable due to their limited depth) to fit the negative gradient of the loss function with respect to the preceding model's predictions. Under squared-error loss, this is mathematically equivalent to sequentially fitting the residuals of the current ensemble. #### Bootstrap Aggregation (Bagging) & Random Forests Bagging (bootstrap aggregation) averages predictions across $B$ bootstrapped training sets to decrease prediction variance. In standard bagging and random forests, these individual trees are fully grown to keep bias as low as possible, and their resulting high variance is mitigated by averaging the predictions: $$\hat{f}_{bag}(x) = \frac{1}{B} \sum_{b=1}^{B} \hat{f}^{*b}(x)$$ Random Forests improve upon bagging by decorrelation. If a dataset contains a single highly dominant predictor, the trees grown on different bootstrap samples will still look very similar, resulting in highly correlated predictions and limiting the variance reduction achieved by averaging. To prevent this, the Random Forest algorithm restricts each split to a random subset of $m \approx \sqrt{p}$ (or $p/3$ for regression) predictors. This decorrelates the individual trees, allowing averaging to reduce ensemble variance much more effectively. #### Gradient Boosted Trees Unlike Random Forests, which grow independent trees in parallel, Boosting grows trees sequentially. The pseudo-code below describes the classical gradient boosting algorithm for regression. Note that the sequential fitting of the residuals ($r_i$) implies that we are assuming a squared-error (MSE) loss function: 1. Set $\hat{f}(x) = \bar{y}$ and $r_i = y_i-\bar{y}$ for all $i$. 2. For $b = 1, 2, \dots, B$: a. Fit a tree $\hat{f}^b$ with $d$ splits to the training data $(X, r)$. b. Update the ensemble: $\hat{f}(x) \leftarrow \hat{f}(x) + \eta \hat{f}^b(x)$, where $\eta \in (0,1]$ is the learning rate (shrinkage). c. Update the residuals: $r_i \leftarrow r_i - \eta \hat{f}^b(x_i)$. 3. The final boosted model is: $\hat{f}(x) = \sum_{b=1}^{B} \eta \hat{f}^b(x)$. #### Generality of Ensemble Strategies While bagging and boosting are most famous for their implementations with decision trees, they are general meta-algorithms that can be applied to any base learner (such as neural networks, support vector machines, or linear models). For example, the `{xgboost}` implementation is not strictly limited to decision trees; it supports alternative base learners (e.g., regularized linear regressions) as the models fitted in step 2a. --- ### 2.2 Random Forest Specification & Tuning We specify a Random Forest using `rand_forest()`, setting the engine to `ranger` and `importance = "impurity"`. This setting is required to calculate the importance scores needed for our downstream `{vip}` plots in Section 2.4, which default to `"none"`. > Besides the engine's name, and the model passed through the pipe operator `%>%`, extra named arguments passed to `set_engine()` are sent directly to the selected engine. In this setup, we optimize model performance by focusing on `min_n` and keeping the remaining parameters fixed: * `min_n`: The minimum observation threshold to attempt a split, which we tune to control tree complexity. * `mtry`: The number of randomly sampled predictors at each split, which we fix at `floor(p/3)` (164 candidate predictors). As noted in the mathematical background, this is a standard rule of thumb for regression trees that saves substantial computational time by avoiding tuning an extra hyperparameter. * `trees`: The number of trees grown in the ensemble, which we fix at 500. Unlike parameters that control individual tree size, increasing the number of trees does not lead to overfitting. Once the ensemble is large enough to stabilize the variance reduction, adding more trees only raises computational runtime without improving predictive performance. ::: {.panel-tabset} ##### 1: Specification & Grid ```r # Specify Random Forest model rf_spec <- rand_forest( mtry = floor(493/3), trees = 500, min_n = tune() ) %>% set_engine("ranger", importance = "impurity") %>% set_mode("regression") rf_workflow <- workflow() %>% add_recipe(fredmd_recipe) %>% add_model(rf_spec) # Define tuning grid rf_grid <- grid_regular( min_n(range = c(5, 50)), levels = 5 ) ``` ##### 2: Model Tuning ```r # library(future) # plan(multisession, workers = parallel::detectCores() - 4) set.seed(456) # Tune Random Forest over the folds rf_results <- tune_grid( rf_workflow, resamples = folds, grid = rf_grid, metrics = metric_set(rmse, mae), control = control_grid(save_pred = TRUE, verbose = TRUE) ) # plan(sequential) # Saving the tuned model save(rf_results, rf_spec, rf_workflow, file = "data/tuned_rf.RData") ``` ##### 3: Finalize & Fit ```r # Quick inspection of the hyperparameter tuning process autoplot(rf_results) # Select best parameters best_rf_params <- select_best(rf_results, metric = "rmse") # Finalize and fit model final_rf_wf <- finalize_workflow(rf_workflow, best_rf_params) set.seed(456) final_rf_fit <- fit(final_rf_wf, data = tail(treated_data, 360)) ``` ::: --- ### 2.3 Gradient Boosting Unlike Random Forests, Gradient Boosted Trees (often simply referred to as Gradient Boosting) are highly susceptible to overfitting if the number of sequential boosting rounds, $B$, is too large. To mitigate this risk, we optimize two key hyperparameters: the learning rate (`learn_rate`) and the tree depth (`tree_depth`). Additionally, a widely used strategy in practice is to implement early stopping using a held-out validation sample. Early stopping monitors the model's performance on this validation set at each boosting iteration and halts training if the validation loss fails to improve for a specified number of consecutive rounds. This combination of highly non-linear base learners and sequential updates makes Gradient Boosting extremely powerful. Although highly effective, we omit early stopping here to avoid look-ahead bias. In `{tidymodels}`, the default early-stopping wrapper partitions the validation set randomly rather than chronologically. While we could bypass this wrapper to force a chronological split, doing so is difficult to integrate with the unified resampling pipeline. This highlights a key trade-off of high-level frameworks: we trade engine-specific customization for pipeline simplicity. In what follows, we control overfitting solely by tuning the learning rate, tree depth, and boosting rounds. ::: {.panel-tabset} ##### 1: Specification & Grid ```r # Specify XGBoost model xgb_spec <- boost_tree( trees = tune(), tree_depth = tune(), learn_rate = tune() ) %>% set_engine("xgboost") %>% set_mode("regression") xgb_workflow <- workflow() %>% add_recipe(fredmd_recipe) %>% add_model(xgb_spec) # Define tuning grid xgb_grid <- grid_regular( trees(range = c(500, 2000)), tree_depth(range = c(1, 6)), learn_rate(range = c(0.005, 0.1), trans = NULL), levels = 3 ) ``` ##### 2: Model Tuning ```r # plan(multisession, workers = parallel::detectCores() - 4) set.seed(789) # Tune XGBoost over the folds xgb_results <- tune_grid( xgb_workflow, resamples = folds, grid = xgb_grid, metrics = metric_set(rmse, mae), control = control_grid(save_pred = TRUE, verbose = TRUE) ) # plan(sequential) # Saving the tuned model save(xgb_results, xgb_spec, xgb_workflow, file = "data/tuned_xgb.RData") ``` ##### 3: Finalize & Fit ```r # Quick inspection of the hyperparameter tuning process autoplot(xgb_results) # Select best parameters best_xgb_params <- select_best(xgb_results, metric = "rmse") # Finalize and fit model final_xgb_wf <- finalize_workflow(xgb_workflow, best_xgb_params) set.seed(789) final_xgb_fit <- fit(final_xgb_wf, data = tail(treated_data, 360)) ``` ::: --- ### 2.4 Variable Importance Some modeling approaches, such as linear regressions or decision trees, are inherently interpretable as we can easily grasp what the estimated coefficients, representing partial effects, or sequential splitting rules represent for the relationships being studied. This transparency is lost, however, when using more general non-linear models - including tree ensembles like Random Forests and Gradient Boosting that aggregate hundreds of individual trees. To understand what happens inside these estimators, we must use ad-hoc post-estimation methods. Variable Importance (VI) solves this by measuring how much each predictor contributes to reducing the model's prediction error. Specifically, tree ensembles typically compute variable importance using either impurity reduction (i.e., reduction in node variance) or gain (measuring improvement in training loss). While these metrics are specific to trees, more general alternatives, such as permutation-based importance, can be used to evaluate feature contribution across any type of machine learning model. Here we show how to use the `{vip}` package to extract and plot the top 20 predictors for our finalized models. ::: {.panel-tabset} ##### Variable Importance Plots ```r library(vip) library(gridExtra) # Plot Random Forest VIM p_rf <- final_rf_fit %>% extract_fit_parsnip() %>% vip( num_features = 20, geom = "col", aesthetics = list(fill = "darkred") ) + theme_minimal() + labs( title = "Random Forest", subtitle = "Top 20 predictors (Impurity)", y = "Importance", x = "Predictor" ) # Plot XGBoost VIM p_xgb <- final_xgb_fit %>% extract_fit_parsnip() %>% vip( num_features = 20, geom = "col", aesthetics = list(fill = "darkorange") ) + theme_minimal() + labs( title = "XGBoost", subtitle = "Top 20 predictors (Gain)", y = "Importance", x = "Predictor" ) grid.arrange(p_rf, p_xgb, ncol = 2) ``` :::