Tree-Based Ensembles: Random Forests and Gradient Boosting

In plain language

Tree-based ensembles are a family of machine-learning methods that build hundreds of simple decision trees and combine their predictions into one stable, accurate answer. Random forests build each tree on a different random sample of the data so the trees disagree with each other in useful ways; averaging them out cancels individual errors. Gradient boosting builds trees in sequence, each one correcting the mistakes of the previous one, a bit like a team where each new expert fixes what the last expert got wrong. In health research, these tools are used to flag patients at high risk of hospitalization or complications — but a high risk score does not mean any single factor causes that risk, and the raw score must be calibrated before treating it as a true probability.

How decision trees work — and why single trees fail

A classification and regression tree (CART) partitions the feature space by recursively splitting on the single feature and threshold that maximally reduces a loss criterion — Gini impurity or cross-entropy for classification, mean squared error for regression. Each split creates two child nodes; the tree grows until a stopping rule fires (maximum depth, minimum leaf size, or minimum impurity reduction). Leaf nodes contain the training patients that reached that partition, and the leaf prediction is the majority class or mean outcome in those patients.

Single trees are high-variance estimators: a small change in the training data can route patients down entirely different branches, producing a very different tree — and very different predictions — from nearly identical datasets. This instability is the primary motivation for ensembles, which average over many perturbed trees to smooth away variance while keeping bias roughly comparable to a single large tree.

Random forests: bagging plus random feature subsampling

Breiman's random forest (2001) combines two ideas to build a stable ensemble.

First, bootstrap aggregating (bagging): draw B bootstrap samples of size n with replacement from the training data. Fit one full tree per bootstrap sample. Average the B leaf predictions (for classification, average the predicted probabilities; for regression, average the leaf means). Because each tree sees a different realization of the training data, their errors are partly uncorrelated — the ensemble variance is approximately sigma^2 times (rho + (1-rho)/B), where rho is the pairwise tree correlation. More trees and lower correlation both reduce ensemble variance.

Second, random feature subsampling: at every split, randomly select m features from the total p and find the best split only within that subset (default: m = sqrt(p) for classification, m = p/3 for regression). Without feature subsampling, a single dominant predictor appears in every tree's top splits, making trees near-identical and destroying the variance-reduction benefit of bagging. Feature subsampling decorrelates trees while introducing only a modest increase in individual-tree bias.

Out-of-bag (OOB) error is a free built-in validation estimate. Each bootstrap sample leaves out roughly 37% of training patients. For each training patient, average predictions only from trees where that patient was out-of-bag — then compute error on those OOB predictions. This is an approximately unbiased estimate of generalization error without requiring a separate validation split.

Key hyperparameters and sane defaults for random forests: n_estimators (500-1,000; more trees always reduce variance and never hurt, only increase compute), min_samples_leaf (5-20 for stability with rare events), max_features (sqrt(p) for classification). Random forests are notably close to tuning-free compared to gradient boosting.

Gradient boosting: sequential residual fitting

Friedman's gradient boosting machine (2001) builds trees sequentially, each one fitting the negative gradient of the loss (the pseudo-residuals) left by the current ensemble. The prediction after M trees is a weighted sum: F_M(x) = F_0(x) + learning_rate times the sum of h_m(x) for m = 1 to M, where F_0 is a simple baseline (the global log-odds or mean), learning_rate is the step size (eta, typically 0.01 to 0.1), and each h_m is a shallow tree (depth 2-8) fit to the residuals from the current model. By taking small steps and adding many trees, boosting progressively reduces both bias and variance.

XGBoost (Chen and Guestrin, 2016) extends vanilla gradient boosting with second-order gradient approximations, column and row subsampling per tree, L1/L2 regularization on leaf weights, and an efficient tree pruning criterion. These additions together yield substantially faster training and often higher accuracy than vanilla gradient boosting on tabular data — the regime in which claims and EHR feature matrices fall.

Boosting requires more active hyperparameter management than random forests. The key parameters are learning_rate (lower is safer but requires more trees), n_estimators, max_depth (shallower trees generalize better: depth 3-6 for most tasks), and early stopping, which monitors performance on a validation fold and halts tree addition once performance has not improved for a user-defined patience window. Early stopping is the primary defense against overfitting in gradient boosting and should always be used.

Interpreting the output

Consider a random forest trained on 50,000 Medicare beneficiaries with heart failure, predicting 1-year hospitalization. Patient A receives a predicted risk score of 0.72. The model's c-statistic is 0.81. The top permutation variable importance (VIMP) feature is prior inpatient admissions in the preceding 12 months.

(1) Formal interpretation.

A raw random forest score of 0.72 is NOT equivalent to "this patient has a 72% probability of hospitalization." Random forests are known to be poorly calibrated by default: averaging over many trees compresses predicted probabilities toward the ensemble mean, away from 0 and 1. Before treating 0.72 as a probability, calibration must be assessed on a held-out dataset (calibration slope, calibration-in-the-large, or the integrated calibration index), and recalibration must be applied if calibration is poor — via isotonic regression or Platt scaling (CalibratedClassifierCV in scikit-learn). After confirmed recalibration, 0.72 represents the model's estimated conditional probability of hospitalization given this patient's observed feature values, conditional on those features being complete and the model being transportable to this population.

The c-statistic of 0.81 means that if a pair of patients is drawn at random — one who was hospitalized and one who was not — the model assigns the higher score to the hospitalized patient 81% of the time. This is a measure of rank discrimination, not calibration. A model can have excellent discrimination (high c-statistic) and poor calibration simultaneously.

Prior inpatient admissions appearing at the top of the variable importance ranking reflects its predictive contribution to the model — it does NOT mean that preventing or reducing admissions (or changing the feature value by intervention) would lower the patient's risk. Variable importance measures the average decrease in model accuracy when that feature is permuted in the data; it is a property of the fitted predictive model, not a causal effect estimate. Patients with high prior admissions may simultaneously have unmeasured severity, frailty, and social determinants that drive both the feature and the outcome. Asserting that prior admissions cause future hospitalizations solely from VIMP is an unsupported causal claim and a common analytic error.

(2) Practical interpretation.

A score of 0.72 places this patient roughly in the top decile of predicted risk in this population. The actionable message is: "the model flags this patient as high-risk, primarily based on prior hospitalization history." Prior admissions drive the prediction — that does not mean preventing or reducing those admissions changes the underlying risk. Use the score for care management triage, risk stratification, or as a covariate in downstream analyses, not as evidence that any individual feature is a modifiable cause of future hospitalization.

RWE-specific considerations

High-cardinality claims features. ICD-10 diagnosis codes (70,000+ unique), NDC drug codes (100,000+), and CPT procedure codes present a high-cardinality categorical challenge. Standard one-hot encoding is infeasible at this scale. Practical approaches: aggregate codes to drug class or body-system hierarchy, use ever/never binary flags for curated code groups, or apply target encoding (substituting the within-training-fold mean outcome for each code level). Target encoding is powerful but carries a leakage trap: if encoding is computed on the full dataset before train/test splitting, it leaks outcome information into the features and inflates performance. Target encoding must always be fit on the training fold only, never on the full dataset.

Missing data. Scikit-learn RF requires explicit imputation before training. The ranger and randomForestSRC packages in R implement surrogate splits — when a feature is missing for a prediction, the tree routes through the next-best alternative split on a correlated feature. Surrogate splits are more principled than median imputation but their behavior must be explicitly verified, especially for informative missingness patterns (e.g., lab values missing for the lowest-acuity patients).

Temporal leakage in claims. The most destructive failure in claims-based ML is temporal leakage: including features from claims dated after the index date in what is labeled as a baseline feature matrix. In a hospitalization prediction model indexed on the first heart failure diagnosis date, a claim for a post-diagnosis echocardiogram is future information. A model that sees this claim learns to predict hospitalization from a feature that was itself caused by the outcome — inflating AUC dramatically (sometimes from 0.65 to 0.90+) while producing useless predictions in prospective deployment. Temporal leakage is invisible to performance metrics computed on the same contaminated dataset. The required guard: enforce a strict feature extraction window of [index_date minus lookback, index_date) — with index_date excluded — and audit every feature's construction date before training.

Calibration drift across databases. A random forest trained on commercial claims may be well-calibrated in-sample but systematically mis-predict in a Medicare or Medicaid population. Validate calibration — not just AUC — in every target population before deployment. See the companion concept prediction-model-validation-recalibration-rwe for the external validation workflow. Cross-database calibration checking is especially important for models exported from development cohorts to support regulatory submissions.

Role in causal pipelines. Random forests and gradient boosting serve as nuisance estimators — outcome models and propensity score estimators — inside doubly robust causal methods (TMLE and DML), as described in the parent concept predictive-and-causal-ml-models-rwe. When used this way, the ensemble delivers conditional expectations to a causal estimator; it is a means to an end, not the deliverable. The ensemble's predictive accuracy matters for efficiency; the causal estimate's validity comes from the doubly robust structure, not from the ensemble alone.

Pros, cons, and trade-offs

Pros. Ensembles automatically capture non-linear relationships and feature interactions without manual feature engineering. Random forests are nearly tuning-free and robust to outliers (individual tree predictions from extreme values are averaged away). Gradient boosting achieves best-in-class accuracy on structured tabular data — precisely the format of claims and EHR feature matrices. OOB error gives random forests a built-in internal validation estimate at no additional compute cost. Variable importance metrics identify the most predictive features for subsequent biological or clinical review. Both methods handle feature matrices with thousands of binary indicators (one ICD flag per code) without explicit variable selection. Both scale well to large cohort sizes common in Medicare and commercial claims.

Cons. Neither random forests nor gradient boosting are calibrated out of the box — raw predicted probabilities compress toward the mean and require explicit recalibration before use as risk estimates. Gradient boosting is sensitive to hyperparameters (especially learning rate and n_estimators) and always requires early stopping to prevent overfitting. Both methods are functionally black-box: predictions are not decomposable into a simple equation, limiting local interpretability compared to logistic regression (though SHAP values and partial dependence plots partially address this). Variable importance rankings are biased toward high-cardinality and correlated features (e.g., many correlated ICD codes for the same condition inflate that condition's apparent importance). External validation and recalibration are mandatory before clinical or regulatory use.

When to use

Use tree-based ensembles when: (1) the goal is risk stratification, case-finding, or covariate adjustment and a logistic regression model is likely to underfit because of non-linearity or unmeasured interactions; (2) the feature space is high-dimensional — hundreds to thousands of ICD/NDC/CPT indicator flags — and automatic feature selection is needed; (3) a built-in internal validation estimate is desired (OOB error for RF) or early stopping is feasible (gradient boosting); (4) ensembles serve as nuisance estimators inside TMLE or DML for causal inference where the parent concept provides the estimator structure; (5) a pre-specified variable importance ranking is required for exploratory reporting, with the explicit caveat — documented in the analysis plan — that importance is not causation.

When NOT to use

Do NOT use tree-based ensembles when: (1) the deliverable is an interpretable causal effect estimate with a confidence interval on a well-defined estimand — use regression, propensity methods, or doubly robust causal estimators instead; (2) the sample is very small (fewer than a few hundred patients) or the event is rare and careful tuning and evaluation are infeasible; (3) the use case requires reporting a hazard ratio, odds ratio, or rate ratio with a confidence interval to regulators or payers — regression models deliver this directly and transparently; (4) variable importance scores are being used as evidence of causal importance or proposed as biomarkers without external causal validation — VIMP measures predictive contribution only; (5) the model has not been externally validated and recalibrated in the target population — never deploy a model developed on commercial claims to Medicare beneficiaries without cross-population calibration assessment. Do not use temporal-leakage-contaminated features as a shortcut to inflated performance metrics.

Worked example

Scenario

A health analytics team trains a random forest with five trees to predict 1-year hospitalization in a small heart failure cohort. They want to compute the ensemble's predicted probability for one patient (P-001) and estimate the model's out-of-bag error on 100 OOB patients. All five trees had P-001 in their out-of-bag set for this illustration, so each tree provides an independent prediction.

Dataset

Predicted probabilities from five bootstrap trees for patient P-001 (all five trees had this patient in their OOB set). A separate OOB pool of 100 patients had 25 classification errors across the forest.

Steps

• Five bootstrap trees provide predicted probabilities for patient P-001: tree T1 = 0.6, T2 = 0.8, T3 = 0.4, T4 = 0.7, T5 = 0.5. These five values represent five independent looks at the same patient, each from a tree trained on a different bootstrap sample.

• Random forest ensemble probability for P-001 = (0.6+0.8+0.4+0.7+0.5)/5 = 3.0/5 = 0.6. Because 0.6 exceeds the classification threshold of 0.5, the model flags this patient as high risk for 1-year hospitalization.

• For OOB error estimation: each of the 20 patients per tree who were not in that tree's bootstrap sample contribute one OOB prediction. Pooling across all five trees gives 100 OOB predictions total. Of these 100 OOB predictions, 25 disagreed with the true hospitalization label. OOB classification error = 25/100 = 0.25, meaning the forest misclassifies 25% of OOB patients.

• The raw score of 0.6 for P-001 is not yet a calibrated probability of hospitalization. Random forests compress predicted probabilities toward the ensemble mean; a calibration check on held-out data (calibration slope, integrated calibration index) is required. If the calibration slope is below 1.0, the model over-predicts at low risk and under-predicts at high risk, and isotonic recalibration must be applied before the score can be reported as a probability.

Result

Ensemble mean score = (0.6+0.8+0.4+0.7+0.5)/5 = 3.0/5 = 0.6 (patient P-001 flagged as high risk). OOB classification error = 25/100 = 0.25. The raw 0.6 score requires calibration assessment before being reported as a probability of hospitalization.

Runnable example

import numpy as np
from sklearn.ensemble import RandomForestClassifier, HistGradientBoostingClassifier
from sklearn.calibration import CalibratedClassifierCV, calibration_curve
from sklearn.model_selection import train_test_split
from sklearn.inspection import permutation_importance
from sklearn.metrics import roc_auc_score, brier_score_loss
import xgboost as xgb

# Assumes: X = feature matrix (pre-index features only — no temporal leakage),
#          y = binary outcome (1=event, 0=no event), feature_names = list of column names.

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=42)

# ── 1. Random forest with OOB error and variable importance ──
rf = RandomForestClassifier(
    n_estimators=500,         # more trees = lower variance; 500 is a safe default
    max_features="sqrt",      # Breiman default for classification
    min_samples_leaf=10,      # stability guard for rare events
    oob_score=True,           # enables the free OOB accuracy estimate
    n_jobs=-1,
    random_state=42,
)
rf.fit(X_train, y_train)
print(f"RF OOB accuracy: {rf.oob_score_:.3f}")  # free internal validation
print(f"RF val AUC: {roc_auc_score(y_val, rf.predict_proba(X_val)[:, 1]):.3f}")

# Permutation importance on validation set (less biased than MDI for correlated features)
perm_imp = permutation_importance(rf, X_val, y_val, n_repeats=10, random_state=42)
top5_idx = perm_imp.importances_mean.argsort()[-5:][::-1]
print("Top 5 features (permutation VIMP):")
for i in top5_idx:
    print(f"  {feature_names[i]}: {perm_imp.importances_mean[i]:.4f}")
print("CAUTION: VIMP = predictive contribution, NOT causal effect.")

# Calibration check — random forests are miscalibrated by default
prob_true, prob_pred = calibration_curve(y_val, rf.predict_proba(X_val)[:, 1], n_bins=10)
print(f"RF Brier score (val): {brier_score_loss(y_val, rf.predict_proba(X_val)[:,1]):.4f}")
# If calibration curve deviates from the diagonal, recalibrate:
rf_cal = CalibratedClassifierCV(rf, method="isotonic", cv="prefit")
rf_cal.fit(X_val, y_val)  # isotonic recalibration on the validation set

# ── 2. Gradient boosting with early stopping (scikit-learn) ──
gbm = HistGradientBoostingClassifier(
    learning_rate=0.05,
    max_iter=1000,            # upper bound; early stopping will stop earlier
    max_depth=4,
    min_samples_leaf=20,
    early_stopping=True,
    validation_fraction=0.15,
    n_iter_no_change=20,      # stop if no improvement for 20 rounds
    random_state=42,
)
gbm.fit(X_train, y_train)
print(f"GBM trees used (early stopping): {gbm.n_iter_}")
print(f"GBM val AUC: {roc_auc_score(y_val, gbm.predict_proba(X_val)[:, 1]):.3f}")

# ── 3. XGBoost with early stopping and calibration ──
dtrain = xgb.DMatrix(X_train, label=y_train)
dval   = xgb.DMatrix(X_val,   label=y_val)
params = {
    "objective": "binary:logistic",
    "eval_metric": "auc",
    "learning_rate": 0.05,
    "max_depth": 4,
    "subsample": 0.8,
    "colsample_bytree": 0.8,
    "seed": 42,
}
xgb_model = xgb.train(
    params, dtrain,
    num_boost_round=1000,
    evals=[(dval, "val")],
    early_stopping_rounds=20,
    verbose_eval=False,
)
xgb_preds = xgb_model.predict(dval)
print(f"XGB val AUC: {roc_auc_score(y_val, xgb_preds):.3f}")
print(f"XGB Brier score (val): {brier_score_loss(y_val, xgb_preds):.4f}")
# Always check calibration on the target population before deployment.

library(ranger); library(xgboost); library(pROC)

# Assumes: dat = data.frame with columns [outcome_col, feature_cols...]
# All feature_cols must be from the pre-index baseline window (no temporal leakage).

set.seed(42)
n <- nrow(dat)
val_idx   <- sample(n, size = floor(0.2 * n))
train_dat <- dat[-val_idx, ]
val_dat   <- dat[ val_idx, ]

# ── 1. Random forest via ranger ──
rf_form <- reformulate(feature_cols, response = outcome_col)
rf_fit <- ranger(
  formula       = rf_form,
  data          = train_dat,
  num.trees     = 500,
  mtry          = floor(sqrt(length(feature_cols))),   # Breiman default
  min.node.size = 10,
  probability   = TRUE,          # return predicted probabilities
  importance    = "permutation", # permutation VIMP (less biased than impurity)
  seed          = 42
)
cat(sprintf("RF OOB prediction error: %.3f\n", rf_fit$prediction.error))

# Variable importance (VIMP) -- predictive, NOT causal
vimp <- sort(rf_fit$variable.importance, decreasing = TRUE)
cat("Top 5 features (permutation VIMP):\n")
print(head(vimp, 5))
cat("NOTE: VIMP reflects predictive contribution only -- not causal effect.\n")

# Validation AUC
val_probs <- predict(rf_fit, data = val_dat)$predictions[, "1"]
rf_auc    <- auc(roc(val_dat[[outcome_col]], val_probs, quiet = TRUE))
cat(sprintf("RF validation AUC: %.3f\n", rf_auc))

# Calibration check: plot calibration curve; recalibrate if slope differs from 1.0
# Use the rms package (val.prob) or a loess smoother for a calibration curve in practice.

# ── 2. Gradient boosting via xgboost with early stopping ──
X_train <- as.matrix(train_dat[, feature_cols])
y_train <- train_dat[[outcome_col]]
X_val   <- as.matrix(val_dat[, feature_cols])
y_val   <- val_dat[[outcome_col]]

dtrain <- xgb.DMatrix(X_train, label = y_train)
dval   <- xgb.DMatrix(X_val,   label = y_val)

params <- list(
  objective        = "binary:logistic",
  eval_metric      = "auc",
  eta              = 0.05,   # learning rate (lower = safer, needs more rounds)
  max_depth        = 4,
  subsample        = 0.8,
  colsample_bytree = 0.8
)

# Cross-validated tuning to select n_rounds (use in model development phase)
cv_res <- xgb.cv(
  params    = params,
  data      = dtrain,
  nrounds   = 500,
  nfold     = 5,
  early_stopping_rounds = 20,
  verbose   = 0
)
best_rounds <- cv_res$best_iteration
cat(sprintf("XGB optimal rounds (CV early stopping): %d\n", best_rounds))

# Final model on full training data at the optimal rounds
xgb_fit <- xgb.train(
  params     = params,
  data       = dtrain,
  nrounds    = best_rounds,
  watchlist  = list(val = dval),
  verbose    = 0
)
xgb_preds <- predict(xgb_fit, dval)
xgb_auc   <- auc(roc(y_val, xgb_preds, quiet = TRUE))
cat(sprintf("XGB validation AUC: %.3f\n", xgb_auc))

# ── 3. Using as nuisance in causal inference (sketch) ──
# For TMLE or DML, pass the ranger or xgboost learner as the outcome or propensity
# nuisance. See predictive-and-causal-ml-models-rwe for the full cross-fitting workflow.
# library(grf); causal_forest(X, Y, W, num.trees=2000) uses RF internally.
# library(DoubleML); lrn("classif.xgboost") wraps xgboost as a mlr3 learner for DML.