ml-003 · Random forest vs gradient boosting — R

Problem statement

Open this practice in GitHub Codespaces runs in your browser; open this folder under practices/ after it starts

Simulate a nonlinear tabular regression task in-script (no files): n = 2000 rows, 8 features drawn Uniform(-2, 2), outcome y = 4*sin(1.5*f1) + 3*f2*f3 + 2*f4^2 - 1.5*f5 + noise, with f6-f8 pure noise features, seed 42. On one fixed train/validation split (80/20), compare a ranger random forest (near defaults) against xgboost — via the xgb.train() + xgb.DMatrix() API — at eta 0.30 and 0.05 across ensemble sizes 25-400, printing a validation-RMSE table. Then fit xgboost once more with early_stopping_rounds and report the best iteration, and finally compare the two models’ top-3 feature importances. You know it is correct when: both models beat the predict-the-mean baseline, the RF column barely moves across sizes while the xgboost columns move a lot, early stopping halts before the round budget, and f6-f8 sit at the bottom of both importance rankings.

Learning objectives

After completing this practice you should be able to:

• explain bagging vs boosting and predict how each responds to more trees
• fit ranger() and xgb.train() on the same split and compare them fairly with a validation-RMSE curve
• explain why RF is robust to its defaults while boosting requires eta and nrounds to be tuned together
• use xgboost early stopping (evals + early_stopping_rounds) and read the best iteration back with xgb.attr()
• compare ranger’s impurity importance with xgb.importance() gain and explain why they can disagree on interaction terms

Concepts and definitions

See concepts for the shared background. R-specific notes: use the xgb.train() interface with xgb.DMatrix objects, not the xgboost(data, label) convenience wrapper — only xgb.train() exposes the evals watch-list and early stopping cleanly. In xgboost R >= 3.0 the monitored sets are passed as evals = list(valid = dvalid) (the argument was renamed from watchlist), and the best round is read with xgb.attr(fit, "best_iteration"). ranger’s importance is off by default; ask for importance = "impurity" when you need it.

Input / output description

Item	Description
Input	Nothing external — data simulated in-script with set.seed(42): a data.frame with columns f1-f8 (Uniform(-2, 2)) and outcome y, 2000 rows
Output	Printed report: baseline RMSE; a 5-row table of validation RMSE (rf_default, xgb_eta_030, xgb_eta_005) by ensemble size; early-stopping best iteration + RMSE; top-3 importances per model and their overlap count

Before-coding checklist

• I can restate the problem without looking at it
• I know what the output should look like for a tiny example
• I have identified the tricky part
• I know how I will check my answer

My attempt

Work here. Run cells interactively; nothing below overwrites this section.

# Milestone 1: simulate the data (n=2000, 8 uniform features, nonlinear
# signal on f1-f5 only, seed 42) and make an 80/20 train/validation split.

# Milestone 2: baseline RMSE (predict the training mean), then fit ranger
# at num.trees in c(25, 50, 100, 200, 400) and record validation RMSE.

# Milestone 3: same sizes for xgb.train (max_depth = 3) at eta 0.30 and
# 0.05; assemble the three curves into one comparison data.frame.

# Milestone 4: refit xgboost with nrounds = 500, evals = list(valid = ...),
# early_stopping_rounds = 30; report the best iteration and its RMSE.

# Milestone 5: compare top-3 feature importances (ranger impurity vs
# xgb.importance gain) and check how much lands on the noise features.

Hint 1

TipHint 1 — a nudge

Write small functions, not one script: make_data, split_data, rmse, then one function per model that takes the train/valid data frames plus an ensemble size and returns a validation RMSE. Build the comparison table with vapply() over c(25, 50, 100, 200, 400). xgboost wants a numeric matrix — write a feature_matrix() helper that drops the y column and calls as.matrix() — and ranger is happy with the formula interface y ~ ..

Hint 2

TipHint 2 — almost a spoiler

For early stopping: xgb.train(params = list(objective = "reg:squarederror", eta = 0.15, max_depth = 3), data = dtrain, nrounds = 500, evals = list(valid = dvalid), early_stopping_rounds = 30, verbose = 0), then as.integer(xgb.attr(fit, "best_iteration")); a later predict(fit, dvalid) already uses the best iteration. For importances, ranger(..., importance = "impurity")$variable.importance is a named vector, and xgb.importance(model = fit) returns a table you can turn into one with setNames(tab$Gain, tab$Feature); sort both, take the top-3 names, and intersect().

Complete solution

NoteOfficial solution — open only after a real attempt

The official, tested script lives at solutions/solution.R. Run it with:

Rscript practices/ml-003-forest-vs-xgboost/solutions/solution.R

Core code (the script adds a main() that prints the full report, plus an autorun guard for the tests):

suppressPackageStartupMessages({
  library(ranger)
  library(xgboost)
})

SEED <- 42L
ENSEMBLE_SIZES <- c(25L, 50L, 100L, 200L, 400L)

make_data <- function(n_rows = 2000L, seed = SEED) {
  set.seed(seed)
  X <- matrix(runif(n_rows * 8L, min = -2, max = 2), ncol = 8L)
  colnames(X) <- paste0("f", 1:8)
  signal <- 4 * sin(1.5 * X[, "f1"]) + 3 * X[, "f2"] * X[, "f3"] +
    2 * X[, "f4"]^2 - 1.5 * X[, "f5"]
  df <- as.data.frame(X)
  df$y <- signal + rnorm(n_rows, mean = 0, sd = 1.5)
  df
}

split_data <- function(df, valid_frac = 0.2, seed = SEED) {
  set.seed(seed + 1L)
  valid_idx <- sample.int(nrow(df), floor(nrow(df) * valid_frac))
  list(train = df[-valid_idx, , drop = FALSE],
       valid = df[valid_idx, , drop = FALSE])
}

rmse <- function(y_true, y_pred) sqrt(mean((y_true - y_pred)^2))

feature_matrix <- function(df) {
  as.matrix(df[, setdiff(names(df), "y"), drop = FALSE])
}

rf_validation_rmse <- function(train_df, valid_df, num_trees) {
  fit <- ranger(y ~ ., data = train_df, num.trees = num_trees, seed = SEED)
  rmse(valid_df$y, predict(fit, data = valid_df)$predictions)
}

xgb_validation_rmse <- function(train_df, valid_df, nrounds, eta) {
  dtrain <- xgb.DMatrix(feature_matrix(train_df), label = train_df$y)
  fit <- xgb.train(
    params = list(objective = "reg:squarederror", eta = eta,
                  max_depth = 3L, seed = SEED),
    data = dtrain, nrounds = nrounds, verbose = 0)
  rmse(valid_df$y, predict(fit, xgb.DMatrix(feature_matrix(valid_df))))
}

validation_curve_table <- function(train_df, valid_df,
                                   sizes = ENSEMBLE_SIZES) {
  data.frame(
    n_estimators = sizes,
    rf_default = vapply(sizes, function(s)
      rf_validation_rmse(train_df, valid_df, s), numeric(1)),
    xgb_eta_030 = vapply(sizes, function(s)
      xgb_validation_rmse(train_df, valid_df, s, 0.30), numeric(1)),
    xgb_eta_005 = vapply(sizes, function(s)
      xgb_validation_rmse(train_df, valid_df, s, 0.05), numeric(1))
  )
}

fit_xgb_early_stopping <- function(train_df, valid_df, eta = 0.15,
                                   max_rounds = 500L, patience = 30L) {
  dtrain <- xgb.DMatrix(feature_matrix(train_df), label = train_df$y)
  dvalid <- xgb.DMatrix(feature_matrix(valid_df), label = valid_df$y)
  fit <- xgb.train(
    params = list(objective = "reg:squarederror", eta = eta,
                  max_depth = 3L, seed = SEED),
    data = dtrain, nrounds = max_rounds,
    evals = list(valid = dvalid),
    early_stopping_rounds = patience, verbose = 0)
  list(model = fit,
       best_iteration = as.integer(xgb.attr(fit, "best_iteration")),
       valid_rmse = rmse(valid_df$y, predict(fit, dvalid)))
}

top_features <- function(named_importance, k = 3L) {
  names(sort(named_importance, decreasing = TRUE))[
    seq_len(min(k, length(named_importance)))]
}

importance_overlap <- function(imp_a, imp_b, k = 3L) {
  length(intersect(top_features(imp_a, k), top_features(imp_b, k)))
}

Line-by-line explanation

NoteWhat each part does

• make_data — a single set.seed(seed) covers both runif (features) and rnorm (noise) because the calls happen in a fixed order; the signal mixes a sine, a pure interaction f2*f3, a square, and a linear term, while f6-f8 never enter — the importance ground truth.
• split_data — seeds with seed + 1L so the split permutation is not the same stream position as the data draw; negative indexing (-valid_idx) gives the complement without any overlap by construction.
• feature_matrix — xgboost only accepts numeric matrices; dropping y by name (not position) keeps it safe if column order ever changes.
• rf_validation_rmse — ranger(y ~ ., ...) with only num.trees and seed set: the exercise is what defaults buy you. predict() on a ranger fit returns a list; the numbers live in $predictions.
• xgb_validation_rmse — xgb.train with an explicit params list; max_depth = 3 gives classic shallow weak learners and verbose = 0 keeps the console clean over 10 fits.
• validation_curve_table — vapply(..., numeric(1)) over the sizes for each column; vapply over sapply so a wrong return shape fails loudly.
• fit_xgb_early_stopping — evals = list(valid = dvalid) names the monitored set (the name appears in the evaluation log); early_stopping_rounds = 30 stops after 30 non-improving rounds and the best round is stored as a model attribute, retrieved with xgb.attr(fit, "best_iteration"). A later predict(fit, dvalid) automatically uses the best iteration, so the reported RMSE matches the log’s minimum.
• top_features / importance_overlap — sort a named vector descending, take the first k names, count the intersection; min(k, length(...)) guards against xgboost dropping never-used features from its table.

Design decisions

• One fixed split, not CV. The lesson is about shapes of curves, not error bars; a single 80/20 split keeps every fit comparable and the runtime in seconds.
• xgb.train + xgb.DMatrix, never xgboost(data, label). The wrapper hides the watch-list, and constructing the DMatrix once per fit makes the train/valid roles explicit — same reason the Python page uses eval_set.
• ranger at (near) defaults vs xgboost at two etas. The asymmetry is the point: the forest gets no tuning on purpose, so the table shows its worst case against boosting’s sensitivity. (seed and num.threads are set for reproducibility and polite CPU use, not accuracy.)
• Refit per ensemble size. O(sizes) redundant work, but each table cell is exactly “this model at this size”; total cost is still ~2 seconds.
• Early stopping at eta 0.15, budget 500, patience 30 — chosen so the demo genuinely stops early (~iteration 360) rather than exhausting the budget, which is the behavior the section teaches.

Common mistakes

• Using evals’ older name watchlist on xgboost R >= 3.0 (or evals on an old 1.x install) — check names(formals(xgb.train)) when the error says the argument is unused.
• Passing a data.frame straight to xgb.DMatrix — it wants a numeric matrix; forgetting to drop y before as.matrix() silently leaks the label into the features and produces a suspiciously perfect RMSE.
• Judging boosting rounds on training RMSE (monotone by construction) instead of a held-out set.
• Forgetting importance = "impurity" in the final ranger call and finding $variable.importance empty (ranger computes no importance by default).
• Comparing the models at a single ensemble size instead of the whole curve — at 25 rounds eta 0.05 looks hopeless, at 400 it looks fine.
• Expecting f2 and f3 to top both importance lists: interaction credit is split between them, so the two models can legitimately disagree.

Alternative solution

NoteAnother way to do it

xgb.cv() replaces the manual split for choosing the round count: it runs k-fold CV internally with the same params, supports early_stopping_rounds, and returns the per-round CV RMSE in $evaluation_log, from which you take the best mean — more honest than a single validation set, at k times the cost. On the forest side, randomForest::randomForest() is the classic drop-in for ranger (same idea, slower on wide data), and ranger(..., importance = "permutation") gives a permutation importance that is more comparable to xgboost’s gain than impurity is. Prefer the CV variant when the data set is small enough that a single 20% validation split is noisy.

Extensions and harder variants

• Replace the single split with xgb.cv() (and manual k-fold for ranger) and add standard errors to the table.
• Sweep eta over {0.01, 0.05, 0.1, 0.3} with early stopping and plot best iteration vs eta — verify the roughly inverse relationship.
• Add a strongly correlated duplicate of f1 and watch both importance rankings redistribute.
• Raise the noise sd from 1.5 to 4 and re-run: does boosting’s edge shrink when the signal-to-noise ratio drops?
• Extract attributes(fit)$evaluation_log from the early-stopping fit and plot validation RMSE by round — mark the best iteration and the patience window after it.

Review questions

Answer out loud before expanding.

CautionQ1. Why does the RF column barely change from 25 to 400 trees while the xgboost columns change a lot?

RF trees are i.i.d. draws averaged together: after a few dozen trees the average has essentially converged, so extra trees only shave a little variance. In boosting each round adds capacity — the ensemble is still moving toward (or past) the optimum — so validation RMSE depends strongly on nrounds until early stopping or overfitting intervenes.

CautionQ2. What exactly does evals = list(valid = dvalid) do, and why is it required for early stopping?

It registers named data sets whose metric xgboost evaluates after every round into the evaluation log. Early stopping watches the last set in that list and stops when its metric has not improved for early_stopping_rounds rounds. Without evals there is no out-of-sample signal to watch — training loss never stalls — so early stopping has nothing to act on.

CautionQ3. xgboost at eta 0.05 with 25 rounds has RMSE near 4.9 — overfitting or underfitting, and how do you tell?

Underfitting: 25 steps of size 0.05 have moved the model only a fraction of the way down the loss, so it is still close to predicting the mean (baseline RMSE about 6). The tell: both training and validation error are high and validation error keeps improving as rounds are added; an overfit model shows low training error with worsening validation error.

CautionQ4. Why is the RMSE at the early-stopping best iteration a slightly optimistic estimate of out-of-sample error?

The validation set chose the stopping round: the reported value is the minimum over ~500 correlated evaluations on that same set — a winner’s-curse effect. Fine for model selection, biased for reporting; use a third untouched test set (or nested CV) for an honest number.

CautionQ5. ranger reports importance for all 8 features but xgb.importance sometimes returns fewer rows. Why?

xgb.importance only lists features that were actually used in at least one split; with shallow trees and early stopping, a pure-noise feature may never be chosen and simply has no row (importance exactly zero). ranger’s deep trees are grown to purity, so every feature — noise included — almost surely gets used somewhere and accumulates a small positive impurity score. The helper’s min(k, length(...)) guard exists precisely for this.

Personal reflection

What was hard? What would you do differently? Date your entries.

References and related practices

• concepts — shared background for this practice
• python.qmd — the same exercise with sklearn + XGBRegressor
• ml-001 (train/validation split) — prerequisite; the split logic here is the same idea
• ml-002 (as-of price join) — sibling demo practice in this lab
• Hastie, Tibshirani, Friedman, The Elements of Statistical Learning, ch. 10 (boosting) and ch. 15 (random forests)
• ?xgb.train and ?ranger — the argument lists worth actually reading