Linear E[y|x] on Boston — coefficients (R)
(Intercept) lstat rm age
-1.175311495 -0.668513052 5.019133481 0.009091327 Trees and Ensembles: A Machine Learning Toolkit for Econometrics
Decision Trees, Random Forests & Gradient Boosting
using R, Python & Stata
28 April 2026
The Linearity Assumption — Relaxing It With Trees
Trees relax linearity entirely. A regression tree estimates \(m(\mathbf{x}) = \mathbb{E}[y\mid\mathbf{x}]\) by partitioning the covariate space into rectangles and fitting a constant in each — capturing interactions, thresholds, and non-linearities without the researcher pre-specifying them.
\[\hat{m}(\mathbf{x}) = \sum_{j=1}^{J} \bar{y}_{R_j}\,\mathbf{1}\{\mathbf{x}\in R_j\}.\]
The picture below shows the idea on two Hitters variables: the plane is cut into boxes, and the prediction is the average outcome inside each box.
Hansen (2022, §29.15) calls a regression tree a “0th-order spline with free knots” — a step function whose break-points are chosen by the data.
Recursive Partitioning — How Would You Estimate \(\mathbb{E}[y\mid\mathbf{x}]\)?
Suppose you may not assume a functional form. How would you estimate the conditional mean? A kernel or series estimator is one answer you have met. A tree is another: split the covariate space \(\mathcal{X}\) into \(J\) non-overlapping rectangles \(R_1,\dots,R_J\) and predict the region mean in each:
\[\hat{m}(\mathbf{x}) = \sum_{j=1}^{J} \bar{y}_{R_j}\,\mathbf{1}\{\mathbf{x}\in R_j\}, \qquad \bar{y}_{R_j} = \frac{1}{|R_j|}\sum_{i:\,\mathbf{x}_i\in R_j} y_i.\]
The formula, in code. “Average \(y\) within each region, then assign that mean to every row in the region” is one grouped line in each language:
Region means from a small tree on Boston (R)
leaf n mean_medv
1 3 175 14.95600
2 4 255 23.34980
3 6 46 32.11304
4 7 30 45.09667Region means from a small tree on Boston (Python)
import pandas as pd
from sklearn.tree import DecisionTreeRegressor
boston = pd.read_csv("../data/boston.csv")
X = boston[["lstat", "rm"]].to_numpy(float); y = boston["medv"].to_numpy(float)
t0 = DecisionTreeRegressor(max_depth=2, random_state=SEED).fit(X, y)
boston = boston.assign(leaf=t0.apply(X)) # leaf id per row
print(boston.groupby("leaf")["medv"].agg(["size", "mean"]).round(2)) size mean
leaf
2 255 23.35
3 175 14.96
5 46 32.11
6 30 45.10Region means by group (Stata)
(encoding automatically selected: ISO-8859-1)
(13 vars, 506 obs)
-------------------------------
| Frequency Mean
--------+----------------------
leaf |
1 | 91 23.67582
2 | 263 17.20837
3 | 128 33.59688
4 | 24 17.5375
Total | 506 22.53281
-------------------------------A tree is therefore nothing exotic: it is a rule for choosing the regions, after which prediction is just a group average.
Recursive Partitioning — CART Growth
How are the regions chosen? Searching all partitions jointly is infeasible, so CART (Breiman et al. 1984) grows the tree greedily by recursive binary splitting. At each node it picks the variable \(v\) and cut-point \(s\) that most reduce the within-node sum of squares:
\[\min_{v,\,s}\;\Big[\underbrace{\sum_{i:\,x_{iv}\le s}(y_i-\bar y_L)^2}_{\text{left RSS}} \;+\; \underbrace{\sum_{i:\,x_{iv}> s}(y_i-\bar y_R)^2}_{\text{right RSS}}\Big].\]
The objective, in code. A single candidate split’s score is a few lines — this is exactly what rpart/sklearn evaluate internally for every \((v,s)\):
RSS of a candidate split, evaluated on Boston (R)
# RSS if we split variable x at cut-point s (smaller = better split)
split_rss <- function(x, y, s) {
left <- y[x <= s]; right <- y[x > s]
sum((left - mean(left))^2) + sum((right - mean(right))^2)
}
# compare two candidate cut-points on lstat — the lower RSS is the better split
c(split_at_10 = split_rss(boston$lstat, boston$medv, 10),
split_at_15 = split_rss(boston$lstat, boston$medv, 15))split_at_10 split_at_15
24111.08 27797.71 RSS of a candidate split, evaluated on Boston (Python)
import pandas as pd
boston = pd.read_csv("../data/boston.csv")
def split_rss(x, y, s):
left, right = y[x <= s], y[x > s]
return ((left - left.mean())**2).sum() + ((right - right.mean())**2).sum()
x, y = boston["lstat"], boston["medv"]
print({f"split_at_{s}": round(split_rss(x, y, s), 1) for s in (10, 15)}){'split_at_10': np.float64(24111.1), 'split_at_15': np.float64(27797.7)}RSS of a candidate split on Boston (Stata)
import delimited "../data/boston.csv", clear
* RSS at lstat <= 10 : Var*(N-1) is the sum of squared deviations within each side
summarize medv if lstat <= 10
scalar rss = r(Var) * (r(N) - 1)
summarize medv if lstat > 10
scalar rss = rss + r(Var) * (r(N) - 1)
display as text "total within-node RSS at split lstat=10: " as result %9.0f rss(encoding automatically selected: ISO-8859-1)
(13 vars, 506 obs)
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
medv | 219 29.47443 8.910574 11.9 50
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
medv | 287 17.23589 4.87689 5 31
total within-node RSS at split lstat=10: 24111Each split is a two-regime sample-split regression (Hansen 2022, §29.15): \(y = \mu_L\,\mathbf{1}\{x_v\le s\} + \mu_R\,\mathbf{1}\{x_v> s\} + e\). The tree is a step function whose break-points are chosen by the data.
Growing and Pruning — Why Prune at All?
A fully grown tree keeps splitting until each leaf is almost pure. It fits the training data nearly perfectly — and therefore overfits: it has memorised noise that will not repeat out of sample.
The remedy is the same bias–variance trade-off you know from Lasso’s \(\lambda\): grow a large tree, then prune it back to the size that predicts best on held-out data.
- Too few leaves → high bias (underfit), the step function is too coarse.
- Too many leaves → high variance (overfit), the steps chase noise.
- The right size sits at the trough of the cross-validation curve.
The next slide writes this trade-off as a penalised objective; the slide after estimates it on real data.
Growing and Pruning — Choosing \(\alpha\) by Cross-Validation
There is no reliable plug-in for \(\alpha\) (unlike some Lasso rules) — cross-validation is the workhorse. We follow the five steps below on Boston.
- Data:
Boston, 506 tracts, targetmedv, 12 predictors (all numeric). - Prepare: none needed — trees handle raw scales.
- Estimate: grow a deliberately large tree with
cp = 0, askingrpartfor a 10-fold CV path (xval = 10). - Read:
rpartreturns acptablewith, per tree size, the CV errorxerrorand its SExstd. The 1-SE rule picks the smallest tree within one SE of the minimum. - Plot: CV error against tree size — a U-shape.
Code — cost-complexity pruning on Boston (R)
set.seed(SEED)
# STEP 3 — grow a large tree and read its 10-fold CV pruning path
big_tree <- rpart(medv ~ ., data = boston, method = "anova",
control = rpart.control(cp = 0, minsplit = 10, xval = 10))
cpt <- as_tibble(big_tree$cptable) # columns: CP, nsplit, rel error, xerror, xstd
# STEP 4 — locate the minimum-CV tree and the 1-SE tree
imin <- which.min(cpt$xerror)
thr <- cpt$xerror[imin] + cpt$xstd[imin] # min CV error + 1 SE
i1se <- which(cpt$xerror <= thr)[1] # smallest tree under that threshold
# STEP 5 — plot CV error vs tree size
ggplot(cpt) +
aes(nsplit + 1L, xerror) +
geom_errorbar(aes(ymin = xerror - xstd, ymax = xerror + xstd),
width = 0.4, colour = col_muted, alpha = 0.6) +
geom_line(colour = col_main, linewidth = 0.9) +
geom_point(colour = col_main, size = 1.6) +
geom_hline(yintercept = thr, linetype = "dashed", colour = col_accent) +
geom_vline(xintercept = cpt$nsplit[i1se] + 1L, linetype = "dotted",
colour = col_ok, linewidth = 0.8) +
labs(x = "Tree size (number of leaves)", y = "Cross-validated relative error",
title = "Boston: cost-complexity pruning by 10-fold CV",
subtitle = "Dashed = min CV error + 1 SE · dotted = the 1-SE-rule tree size") +
theme_lecture
Reading the figure: the curve is U-shaped in tree size. The 1-SE rule (dotted line) picks a parsimonious tree near the trough — the same logic that governed \(\lambda\) for Lasso and Ridge. xerror is relative error (1.0 = a root-only tree), so a value of 0.3 means the tree cut CV error to 30% of the no-split baseline.
Application: the Hitters Tree — Step 5, Plot
Finally, draw the tree. A tree diagram is the model — the splits are the rules, the leaves are the predictions.
Step 5 — plot the tree (Python, sklearn)
import matplotlib.pyplot as plt
from sklearn.tree import plot_tree
fig, ax = plt.subplots(figsize=(9.5, 4.2))
# assign the return value to _ so the list of artists is NOT printed
_ = plot_tree(tree, feature_names=["Years", "Hits"], filled=True, rounded=True,
precision=2, fontsize=9, impurity=False, ax=ax)
_ = ax.set_title("Hitters: regression tree for log(Salary) (sklearn)",
fontsize=11, fontweight="bold")
plt.tight_layout(); plt.show()
Step 5 — a readable tree in Stata
* Stata has no native CART plot. Fit through pystacked and read the rules, or use
* the community command -crtrees-. Here we show the leaf means directly.
import delimited "../data/hitters.csv", clear case(lower)
capture drop lsalary
drop if missing(salary)
generate double lsalary = ln(salary)
* Reproduce the ISLR leaves by hand from the discovered split points:
generate byte leaf = 1 if years < 4.5
replace leaf = 2 if years >= 4.5 & hits < 117.5
replace leaf = 3 if years >= 4.5 & hits >= 117.5
table leaf, statistic(mean lsalary) statistic(frequency)(encoding automatically selected: ISO-8859-1)
(21 vars, 263 obs)
(0 observations deleted)
(173 missing values generated)
(90 real changes made)
(83 real changes made)
-------------------------------
| Mean Frequency
--------+----------------------
leaf |
1 | 5.10679 90
2 | 5.99838 90
3 | 6.739687 83
Total | 5.927222 263
-------------------------------Reading the tree (ISLR §8.1.1). Read top to bottom: the first split is Years < 4.5; experienced players are then split at Hits ≈ 117.5. Each leaf box shows its predicted log(Salary).
Bagging — Out-of-Bag Error for Free
How do you estimate test error without a test set? Use the observations each tree did not see. A bootstrap resample omits on average \(\approx 1/e \approx 37\%\) of the rows — the out-of-bag (OOB) sample. Predict each \(\mathbf{x}_i\) using only the trees that excluded it:
\[\hat{f}^{\text{OOB}}(\mathbf{x}_i)=\frac{1}{|\{b:\,i\notin b\}|}\sum_{b:\,i\notin b} T_b(\mathbf{x}_i).\]
The OOB MSE approximates leave-one-out cross-validation at no extra cost. In code it is a single argument or attribute:
OOB error on Boston (Python)
import pandas as pd
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
boston = pd.read_csv("../data/boston.csv")
feat = [c for c in boston.columns if c != "medv"]
X = boston[feat].to_numpy(float); y = boston["medv"].to_numpy(float)
rf = RandomForestRegressor(n_estimators=N_TREES, oob_score=True,
n_jobs=N_CORES, random_state=SEED).fit(X, y)
print(f"OOB MSE: {mean_squared_error(y, rf.oob_prediction_):.2f}")OOB scoring requested from pystacked (Stata)
(encoding automatically selected: ISO-8859-1)
(13 vars, 506 obs)
Error loading Python Script for pystacked.
unrecognized command
r(199);
r(199);We exploit exactly these OOB predictions as a cheap cross-fitting device in Part V — the bridge from prediction to inference. References: Breiman (1996, 2001); Hansen §29.16; Causal ML Ch. 8.
Random Forests — Why Bagging Is Not Enough
If averaging is so powerful, why not stop at bagging? Because the trees are correlated. If one predictor dominates — lstat in Boston — almost every bootstrapped tree splits on it first, so the trees look alike and their errors do not cancel.
Write the variance of the average via the per-tree variance \(\sigma^2\) and the average pairwise correlation \(\bar\rho\):
\[\operatorname{Var}\!\Big(\tfrac{1}{B}\textstyle\sum_{b=1}^B T_b\Big) = \bar\rho\,\sigma^2 + \frac{1-\bar\rho}{B}\,\sigma^2 \;\xrightarrow{\;B\to\infty\;}\;\bar\rho\,\sigma^2.\]
The second term vanishes as \(B\) grows — but the first does not. The variance floor is \(\bar\rho\,\sigma^2\). Adding trees buys you only the second term; to break through the floor you must lower \(\bar\rho\).
Code — the variance formula, plotted (not a simulation)
# This plots the algebra Var = rho + (1-rho)/B (with sigma^2 = 1), NOT data.
var_curve <- function(rho) tibble(B = 1:300, variance = rho + (1 - rho) / (1:300), rho = rho)
var_df <- bind_rows(var_curve(0.05), var_curve(0.30), var_curve(0.80)) %>%
mutate(label = factor(rho, levels = c(0.05, 0.30, 0.80),
labels = c("rho-bar = 0.05 (forest: decorrelated)",
"rho-bar = 0.30 (bagging: somewhat alike)",
"rho-bar = 0.80 (one dominant predictor)")))
ggplot(var_df) +
aes(B, variance, colour = label) +
geom_line(linewidth = 1) +
scale_colour_manual(values = c(col_ok, col_main, col_accent)) +
scale_y_continuous(limits = c(0, 1)) +
labs(x = "Number of trees B", colour = NULL,
y = expression("Variance of tree-average (per-tree "*sigma^2 == 1*")"),
title = "Lowering between-tree correlation beats adding more trees",
subtitle = "Each curve flattens to its floor rho-bar; only decorrelation pushes the floor down") +
theme_lecture +
theme(legend.position = c(0.62, 0.85))
Past a few hundred trees, more trees barely help — what helps is making the trees different from each other.
Application: Boston — Step 5, Variable Importance Plot
Which predictors does the forest lean on? We plot permutation importance (the honest measure).
Plot permutation importance (R)
tibble(Variable = names(fit_rf$variable.importance),
Importance = fit_rf$variable.importance) %>%
ggplot() +
aes(reorder(Variable, Importance), Importance) +
geom_col(fill = col_main, alpha = 0.85) +
coord_flip() +
labs(x = NULL, y = "Permutation importance (rise in OOB MSE when shuffled)",
title = "Boston: which predictors does the forest rely on?",
subtitle = "lstat (% lower-status pop.) and rm (rooms/dwelling) dominate — cf. ISLR Fig. 8.9") +
theme_lecture
Plot permutation importance (Python)
import matplotlib.pyplot as plt
from sklearn.inspection import permutation_importance
perm = permutation_importance(rf, Xte, yte, n_repeats=10,
random_state=SEED, n_jobs=N_CORES)
order = np.argsort(perm.importances_mean)[::-1]
fig, ax = plt.subplots(figsize=(8, 3.4))
_ = ax.barh([feat[i] for i in order[::-1]], perm.importances_mean[order[::-1]],
color="#1a6ea8", alpha=0.85)
_ = ax.set_xlabel("Permutation importance (rise in test MSE when shuffled)")
_ = ax.set_title("Boston: which predictors does the forest rely on? (lstat & rm dominate)",
fontsize=10.5, fontweight="bold")
ax.grid(True, axis="x", color="#e8e8e8")
plt.tight_layout(); plt.show()
Drop-column importance (Stata, pystacked)
* No native permutation importance in Stata. A clean, runnable analogue is
* leave-one-covariate-out (drop-column) importance: refit without each predictor
* and measure the rise in held-out MSE. Needs: ssc install pystacked
import delimited "../data/boston.csv", clear
set seed 14159
generate byte train = runiform() < 0.5
local X crim zn indus chas nox rm age dis rad tax ptratio lstat
* baseline forest and its held-out MSE
quietly pystacked medv `X' if train, type(reg) methods(rf) pyseed(14159)
predict double f_all
quietly generate double e_all = (medv - f_all)^2 if !train
quietly summarize e_all if !train
scalar mse_all = r(mean)
* refit dropping each predictor in turn; a larger MSE rise = more important
foreach v of local X {
local rest : list X - v
quietly pystacked medv `rest' if train, type(reg) methods(rf) pyseed(14159)
quietly predict double f_`v'
quietly generate double e_`v' = (medv - f_`v')^2 if !train
quietly summarize e_`v' if !train
display as text "`v': drop-column importance = " as result %6.3f (r(mean) - mse_all)
}(encoding automatically selected: ISO-8859-1)
(13 vars, 506 obs)
unrecognized command
r(199);
r(199);Across all three languages, lstat and rm are the two predictors the forest leans on most — consistent with economic intuition (neighbourhood wealth and dwelling size). The Stata tab uses a drop-column measure (refit-without-each-variable) rather than permutation, so with correlated predictors the two can disagree on the smaller variables — but the top of the ranking is robust.
Gradient Boosting — Learning From Mistakes
A forest grows trees independently and averages them. What if, instead, each new tree focused on the errors the current model still makes? That is boosting (Friedman 2001): an additive model built sequentially,
\[\hat{f}^{(m)}(\mathbf{x}) = \hat{f}^{(m-1)}(\mathbf{x}) + \nu\,T_m(\mathbf{x}),\]
where \(T_m\) is a shallow tree fit to the negative gradient of the loss (for squared error, the current residuals) and \(\nu\in(0,1)\) is a learning rate taking small steps so no single tree dominates.
The formula, in code. Plain gradient boosting for squared-error loss is a short loop — fit a stump to the residuals, take a small step, update:
Gradient boosting by hand on Boston (R)
y <- boston$medv
preds <- boston %>% select(-medv) # the predictor columns
f_hat <- rep(mean(y), length(y)); nu <- 0.1 # start at the mean; small learning rate
for (m in 1:50) {
d <- preds %>% mutate(r = y - f_hat) # current residuals = neg. gradient
stump <- rpart(r ~ ., data = d, control = rpart.control(maxdepth = 2))
f_hat <- f_hat + nu * predict(stump, d) # small step in the residual's direction
}
c(start_RMSE = sqrt(mean((y - mean(y))^2)),
boosted_RMSE = sqrt(mean((y - f_hat)^2))) # training error falls as trees are added start_RMSE boosted_RMSE
9.188012 2.779872 Gradient boosting by hand on Boston (Python)
import pandas as pd, numpy as np
from sklearn.tree import DecisionTreeRegressor
boston = pd.read_csv("../data/boston.csv")
feat = [c for c in boston.columns if c != "medv"]
X = boston[feat].to_numpy(float); y = boston["medv"].to_numpy(float)
f_hat = np.full(len(y), y.mean()); nu = 0.1 # start at the mean; small learning rate
for m in range(50):
r = y - f_hat # current residuals = neg. gradient
stump = DecisionTreeRegressor(max_depth=2).fit(X, r)
f_hat += nu * stump.predict(X) # small step in the residual's direction
print(f"start RMSE: {np.sqrt(((y - y.mean())**2).mean()):.2f} | "
f"boosted RMSE: {np.sqrt(((y - f_hat)**2).mean()):.2f}")The same loop, done for you by pystacked’s gradient boosting (Stata)
(encoding automatically selected: ISO-8859-1)
(13 vars, 506 obs)
Error loading Python Script for pystacked.
unrecognized command
r(199);
r(199);xgboost (R/Python) and pystacked’s gradboost (Stata) do this loop with a second-order step and regularisation — but the skeleton above is the algorithm.
Application: Boston XGBoost — Estimate With Early Stopping
Same Boston train/test split. We fit XGBoost with a small eta, a depth-3 base learner, and early stopping, then read the round at which the held-out error stopped improving.
Fit XGBoost + read best round (R)
set.seed(SEED)
xvars <- setdiff(names(boston), "medv")
n <- nrow(boston); train <- sample(n, floor(n / 2))
Xtr <- as.matrix(boston[train, xvars]); ytr <- boston$medv[train]
Xte <- as.matrix(boston[-train, xvars]); yte <- boston$medv[-train]
dtrain <- xgb.DMatrix(Xtr, label = ytr) # XGBoost's native data container
dtest <- xgb.DMatrix(Xte, label = yte)
params <- list(objective = "reg:squarederror", max_depth = 3, eta = 0.05,
subsample = 0.8, colsample_bytree = 0.8, lambda = 1)
fit_xgb <- xgb.train(params, dtrain, nrounds = 1000,
evals = list(train = dtrain, test = dtest), # 'evals' replaces 'watchlist' (xgboost >= 2.1)
early_stopping_rounds = 30, verbose = 0)
mse_xgb <- mean((predict(fit_xgb, dtest) - yte)^2)
# In xgboost >= 2.1 the booster is an altrep object; read the log from R attributes
# and take the best round from the log itself (robust to API changes).
elog <- as_tibble(attributes(fit_xgb)$evaluation_log)
best_it <- elog$iter[which.min(elog$test_rmse)] # round with the lowest held-out RMSE
cat(sprintf("XGBoost: test MSE %.2f (RMSE %.2f) at %d rounds\n",
mse_xgb, sqrt(mse_xgb), best_it))XGBoost: test MSE 15.55 (RMSE 3.94) at 196 roundsFit XGBoost + read best round (R)
# Learning curve: train vs test RMSE per round (tidy reshape, no $ assigns)
eval_long <- elog %>%
pivot_longer(c(train_rmse, test_rmse), names_to = "set", values_to = "rmse") %>%
mutate(set = recode(set, train_rmse = "Train", test_rmse = "Test"))
ggplot(eval_long) +
aes(x = iter, y = rmse, colour = set) +
geom_line(linewidth = 0.8) +
geom_vline(xintercept = best_it, linetype = "dashed", colour = "grey40") +
annotate("text", x = best_it, y = max(eval_long$rmse),
label = paste0(" best round = ", best_it),
hjust = 0, vjust = 1, size = 3.2, colour = "grey30") +
scale_colour_manual(values = c(Train = col_main, Test = col_accent)) +
labs(x = "Boosting round", y = "RMSE", colour = NULL,
title = "XGBoost learning curve on Boston",
subtitle = "Train RMSE keeps falling; test RMSE flattens then rises — early stopping halts at the dashed line") +
theme_lecture
Fit XGBoost + read best round (Python)
import pandas as pd, numpy as np, matplotlib.pyplot as plt
import xgboost as xgb
from sklearn.model_selection import train_test_split
boston = pd.read_csv("../data/boston.csv")
feat = [c for c in boston.columns if c != "medv"]
X = boston[feat].to_numpy(dtype=float); y = boston["medv"].to_numpy(dtype=float)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.5, random_state=SEED)
# early_stopping_rounds + eval_metric live in the constructor (xgboost >= 1.6)
model = xgb.XGBRegressor(objective="reg:squarederror", max_depth=3, learning_rate=0.05,
subsample=0.8, colsample_bytree=0.8, reg_lambda=1.0,
n_estimators=1000, early_stopping_rounds=30, eval_metric="rmse",
random_state=SEED, n_jobs=N_CORES)
model.fit(Xtr, ytr, eval_set=[(Xtr, ytr), (Xte, yte)], verbose=False)
print(f"XGBoost test MSE: {np.mean((model.predict(Xte) - yte)**2):.2f} "
f"(best round {model.best_iteration})")
res = model.evals_result()
fig, ax = plt.subplots(figsize=(8, 3.4))
# assign every artist-returning call to _ so no <Line2D ...> repr leaks to output
_ = ax.plot(res["validation_0"]["rmse"], color="#1a6ea8", lw=1.2, label="Train")
_ = ax.plot(res["validation_1"]["rmse"], color="#e8521a", lw=1.2, label="Test")
_ = ax.axvline(model.best_iteration, color="grey", ls="--", lw=1.0,
label=f"best round = {model.best_iteration}")
_ = ax.set_xlabel("Boosting round"); _ = ax.set_ylabel("RMSE")
_ = ax.set_title("XGBoost learning curve on Boston\ntrain keeps falling; test flattens then rises",
fontsize=10.5, fontweight="bold")
_ = ax.legend(fontsize=9); ax.grid(True, color="#e8e8e8")
plt.tight_layout(); plt.show()
Gradient boosting via pystacked (Stata)
* pystacked wraps scikit-learn's gradient boosting. Needs: ssc install pystacked
import delimited "../data/boston.csv", clear
set seed 14159
generate byte train = runiform() < 0.5
pystacked medv crim zn indus chas nox rm age dis rad tax ptratio lstat ///
if train, type(reg) methods(gradboost) pyseed(14159)
predict double mhat
generate double sq_err = (medv - mhat)^2 if !train
quietly summarize sq_err if !train
display as text "Gradient-boosting test MSE (Stata/pystacked): " as result %6.2f r(mean)> ked
(encoding automatically selected: ISO-8859-1)
(13 vars, 506 obs)
Error loading Python Script for pystacked.
unrecognized command
r(199);
r(199);Reading the learning curve. The gap between train and test RMSE is overfitting. Train RMSE keeps falling; test RMSE flattens then rises. Early stopping halts at the trough (best_iteration) — that round is the chosen n_estimators.
XGBoost — SHAP and Partial Dependence on Boston
The plot code implements the PDP formula directly: a loop over a grid, setting one column to each grid value and averaging the predictions (no map_ helpers — a plain loop, so the operation is transparent).
SHAP summary (R) — reuses fit_xgb, Xtr
shap_long <- shap.prep(xgb_model = fit_xgb, X_train = Xtr) # reuse fit from earlier slide
shap.plot.summary(shap_long, dilute = 2) +
labs(title = "XGBoost SHAP summary on Boston",
subtitle = "Each dot is a tract; colour = feature value (red high, blue low); x = effect on predicted medv") +
theme_lecture
Partial dependence for rm and lstat (R)
# Implements PDP_j(v) = mean over rows of f(x with column j set to v)
pdp_curve <- function(model, X, var, grid_n = 40) {
vals <- seq(min(X[, var]), max(X[, var]), length.out = grid_n)
yhat <- numeric(grid_n); Xtmp <- X
for (k in seq_along(vals)) {
Xtmp[, var] <- vals[k] # set column j to v
yhat[k] <- mean(predict(model, xgb.DMatrix(Xtmp))) # average prediction
}
tibble(value = vals, yhat = yhat, variable = var)
}
bind_rows(pdp_curve(fit_xgb, Xtr, "rm"), pdp_curve(fit_xgb, Xtr, "lstat")) %>%
ggplot() +
aes(value, yhat) +
geom_line(colour = col_main, linewidth = 1.1) +
facet_wrap(~ variable, scales = "free_x") +
labs(x = "Predictor value", y = "Average predicted medv ($000s)",
title = "Partial-dependence plots: average marginal effect on home value",
subtitle = "Predicted price rises with rm (rooms) and falls with lstat (% lower-status) — cf. ISLR §8.3.4") +
theme_lecture
SHAP + partial dependence (Python) — reuses model, Xtr, Xte, feat
import numpy as np, matplotlib.pyplot as plt, shap
explainer = shap.TreeExplainer(model) # reuse fitted model from earlier slide
sv = explainer(Xte)
mean_abs = np.abs(sv.values).mean(axis=0)
order = np.argsort(mean_abs)[::-1]
def pdp_curve(mdl, Xbg, j, grid_n=40): # PDP_j(v): set col j to v, average
grid = np.linspace(Xbg[:, j].min(), Xbg[:, j].max(), grid_n)
out = np.empty(grid_n); Xtmp = Xbg.copy()
for k, v in enumerate(grid):
Xtmp[:, j] = v
out[k] = mdl.predict(Xtmp).mean()
return grid, out
fig, axes = plt.subplots(1, 2, figsize=(11, 3.4))
_ = axes[0].barh([feat[i] for i in order[::-1]], mean_abs[order[::-1]],
color="#1a6ea8", alpha=0.85)
_ = axes[0].set_xlabel("Mean |SHAP value|")
_ = axes[0].set_title("Feature importance by SHAP", fontsize=10.5, fontweight="bold")
axes[0].grid(True, axis="x", color="#e8e8e8")
for name, col in [("rm", "#1a6ea8"), ("lstat", "#e8521a")]:
g, yh = pdp_curve(model, Xtr, feat.index(name))
_ = axes[1].plot(g, yh, lw=1.6, color=col, label=name)
_ = axes[1].set_xlabel("Predictor value"); _ = axes[1].set_ylabel("Avg predicted medv")
_ = axes[1].set_title("Partial dependence: rm up, lstat down", fontsize=10.5, fontweight="bold")
_ = axes[1].legend(fontsize=9); axes[1].grid(True, color="#e8e8e8")
plt.tight_layout(); plt.show()
Critical reading. SHAP and PDP describe how the model predicts, not how the world works. A PDP that slopes up for rm says the fitted function increases in rooms holding other features at their observed values — not a causal “effect of adding a room.” For causal questions we need the econometrics of Part V.
Beyond the Average — Individual Conditional Expectation (ICE)
A PDP shows the average marginal effect — but averaging can hide heterogeneity. If rm matters steeply for some tracts and barely for others, the PDP shows only the mean and looks deceptively smooth. ICE curves (Goldstein et al. 2015) disaggregate the PDP into one curve per observation:
\[\text{ICE}_i(v) = \hat f\big(x_{i1},\dots,x_{ij}=v,\dots,x_{ip}\big), \qquad \text{PDP}_j(v) = \frac{1}{n}\sum_{i=1}^{n}\text{ICE}_i(v).\]
The PDP is literally the average of the ICE curves. When the ICE curves are parallel, the effect is homogeneous and the PDP tells the whole story; when they fan out or cross, the feature interacts with the others — the same heterogeneity that motivates causal forests in Part V.
ICE + PDP for rooms (R) — reuses fit_xgb, Xtr
set.seed(SEED)
# ICE: vary one feature over a grid for a sample of rows, keep EACH row's curve
ice_curves <- function(model, X, var, n_ice = 40, grid_n = 30) {
rows <- sample(nrow(X), min(n_ice, nrow(X)))
vals <- seq(min(X[, var]), max(X[, var]), length.out = grid_n)
Xsub <- X[rows, , drop = FALSE]
parts <- vector("list", length(vals)) # explicit loop, no map_
for (k in seq_along(vals)) {
Xtmp <- Xsub; Xtmp[, var] <- vals[k]
parts[[k]] <- tibble(id = rows, value = vals[k],
yhat = predict(model, xgb.DMatrix(Xtmp)))
}
bind_rows(parts)
}
ice_df <- ice_curves(fit_xgb, Xtr, "rm")
pdp_df <- ice_df %>% summarise(yhat = mean(yhat), .by = value) # PDP = mean of ICE
ggplot(ice_df) +
aes(value, yhat, group = id) +
geom_line(colour = col_muted, alpha = 0.30) +
geom_line(data = pdp_df, aes(value, yhat), inherit.aes = FALSE,
colour = col_accent, linewidth = 1.3) +
labs(x = "rooms per dwelling (rm)", y = "predicted medv ($000s)",
title = "ICE curves (grey) and their average — the PDP (orange)",
subtitle = "Each grey line is one tract; diverging lines would signal an interaction") +
theme_lecture
ICE + PDP for rooms (Python) — reuses model, Xtr, feat
import matplotlib.pyplot as plt
from sklearn.inspection import PartialDependenceDisplay
fig, ax = plt.subplots(figsize=(7.5, 3.4))
# kind="both" overlays the individual ICE curves and their average (the PDP)
_ = PartialDependenceDisplay.from_estimator(
model, Xtr, [feat.index("rm")], kind="both",
subsample=40, random_state=SEED, ax=ax)
_ = ax.set_title("ICE curves and their average (PDP) for rooms (rm)",
fontsize=10.5, fontweight="bold")
plt.tight_layout(); plt.show()
Reading note: parallel ICE curves ⇒ a homogeneous effect; curves that fan out or cross ⇒ an interaction the PDP conceals. ICE is the bridge from average description to the heterogeneous causal effects of Part V.
The Forecast
The model tracks turning points it has never seen — because the patterns in the lagged changes recur, even when the levels do not.
How to Evaluate — Rolling-Origin CV
One train/test split is one draw. The honest protocol re-fits at successive origins and averages the forecast errors — \(k\)-fold’s time-respecting cousin:
\[\text{CV} = \frac{1}{H}\sum_{h=1}^{H} L\!\left(y_{t_h+1},\ \hat y_{t_h+1\mid t_h}\right), \qquad t_1 < t_2 < \cdots < t_H.\]
- Tune hyperparameters (\(\eta\), depth, lags \(p\)) on rolling-origin error, never on the final test window.
- Then one DM test on the held-out period settles XGBoost vs the baseline.
NYSE — Predicted vs Actual
Volume is highly persistent, so lag-1 is already a tough benchmark; XGBoost adds the non-linear lag interactions and the day-of-week effect on top.