Predictors are selected by recursive elimination on random forest importance: fit a forest on all predictors, record the out-of-bag error, drop the least important 20%, and repeat. The configuration minimising the OOB error is chosen.
from analysis.select_features import load_analysis_data, select_features
wide = load_analysis_data(ROOT / "warehouse.duckdb")
path = select_features(wide)
path[["n_variables", "oob_mse"]]| n_variables | oob_mse | |
|---|---|---|
| 0 | 134 | 0.200 |
| 1 | 107 | 0.194 |
| 2 | 85 | 0.186 |
| 3 | 68 | 0.180 |
| 4 | 54 | 0.179 |
| 5 | 43 | 0.175 |
| 6 | 34 | 0.175 |
| 7 | 27 | 0.172 |
| 8 | 21 | 0.165 |
| 9 | 16 | 0.168 |
| 10 | 12 | 0.173 |
| 11 | 9 | 0.182 |
| 12 | 7 | 0.190 |
| 13 | 5 | 0.213 |
| 14 | 4 | 0.228 |
| 15 | 3 | 0.252 |
| 16 | 2 | 0.285 |
| 17 | 1 | 0.493 |
best = path.loc[path.oob_mse.idxmin()]
fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(path.n_variables, path.oob_mse, s=45, alpha=0.6, color="grey")
ax.scatter(best.n_variables, best.oob_mse, s=110, facecolor="none",
edgecolor="black", linewidth=2)
ax.annotate(f"Minimum ({int(best.n_variables)} predictors)",
xy=(best.n_variables, best.oob_mse),
xytext=(best.n_variables + 25, best.oob_mse + 0.05),
arrowprops=dict(arrowstyle="->"))
ax.invert_xaxis()
ax.set_xlabel("Number of predictors")
ax.set_ylabel("Out-of-bag error (MSE)")
ax.spines[["top", "right"]].set_visible(False)
plt.show()
selected = list(best.variables)
Sixty candidate configurations of trees, mtry and min_n are evaluated by five-fold cross-validation on the training set, selecting on RMSE.
from analysis.tune import stratified_split, tune_forest, tuning_results
split = stratified_split(wide, predictors=selected)
search = tune_forest(split)
tuning_results(search).head()Fitting 5 folds for each of 60 candidates, totalling 300 fits| max_features | min_samples_leaf | n_estimators | rmse | rsq | mae | |
|---|---|---|---|---|---|---|
| 0 | 8 | 4 | 1719 | 0.416 | 0.733 | 0.318 |
| 1 | 19 | 3 | 995 | 0.418 | 0.732 | 0.316 |
| 2 | 15 | 4 | 1425 | 0.419 | 0.730 | 0.318 |
| 3 | 9 | 5 | 1942 | 0.422 | 0.725 | 0.322 |
| 4 | 19 | 5 | 1923 | 0.429 | 0.717 | 0.325 |
from analysis.tune import tune_xgboost
xgb_search = tune_xgboost(split)
tuning_results(xgb_search).head()Fitting 5 folds for each of 60 candidates, totalling 300 fits| gamma | learning_rate | max_depth | min_child_weight | rmse | rsq | mae | |
|---|---|---|---|---|---|---|---|
| 0 | 1.762e-03 | 0.008 | 13 | 18 | 0.411 | 0.740 | 0.321 |
| 1 | 5.864e-05 | 0.005 | 11 | 18 | 0.412 | 0.739 | 0.324 |
| 2 | 1.544e-01 | 0.038 | 6 | 21 | 0.419 | 0.730 | 0.325 |
| 3 | 1.904e-09 | 0.012 | 11 | 32 | 0.425 | 0.719 | 0.325 |
| 4 | 1.100e-07 | 0.003 | 12 | 13 | 0.432 | 0.712 | 0.332 |
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
import numpy as np, pandas as pd
lm = LinearRegression().fit(split.X_train, split.y_train)
models = {
"rf": search.best_estimator_,
"lm": lm,
"xgb": xgb_search.best_estimator_,
}
rows = []
for name, model in models.items():
pred = model.predict(split.X_test)
rows.append({
"model": name,
"rmse": np.sqrt(mean_squared_error(split.y_test, pred)),
"rsq": r2_score(split.y_test, pred),
"mae": mean_absolute_error(split.y_test, pred),
})
comparison = pd.DataFrame(rows).round(3)
comparison| model | rmse | rsq | mae | |
|---|---|---|---|---|
| 0 | rf | 0.412 | 0.748 | 0.315 |
| 1 | lm | 0.484 | 0.653 | 0.373 |
| 2 | xgb | 0.439 | 0.713 | 0.327 |
fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharex=True, sharey=True)
for ax, (name, model) in zip(axes, models.items()):
pred = model.predict(split.X_test)
lo, hi = split.y_test.min(), split.y_test.max()
ax.plot([lo, hi], [lo, hi], color="green", linestyle="--")
ax.scatter(pred, split.y_test, alpha=0.4, color="black")
ax.set_title(name)
ax.set_aspect("equal")
ax.spines[["top", "right"]].set_visible(False)
axes[0].set_ylabel("Observed")
fig.supxlabel("Predicted")
plt.tight_layout()
Permutation importance measures how much the prediction error increases when a predictor’s values are shuffled. Values are reported as a ratio to the original error, so 1.0 means the predictor carries no information. Demographic indicators are shown in bold.
from analysis.importance import feature_importance, plot_importance
X_all = wide[selected]
y_all = wide["SUBJ_LIFE_SAT"]
imp_rf = feature_importance(models["rf"], X_all, y_all)
imp_xgb = feature_importance(models["xgb"], X_all, y_all)
fig, axes = plt.subplots(1, 2, figsize=(13, 6))
plot_importance(imp_rf, ax=axes[0], title="Random Forest")
plot_importance(imp_xgb, ax=axes[1], title="XGBoost")
plt.tight_layout()
imp_rf| feature | importance | lower | upper | |
|---|---|---|---|---|
| 0 | EMP_RA | 2.246 | 2.064 | 2.416 |
| 1 | SUBJ_SOC_SUPP | 1.658 | 1.549 | 1.749 |
| 2 | SR_ELD_RA_T | 1.653 | 1.535 | 1.779 |
| 3 | INCOME_DISP | 1.454 | 1.381 | 1.525 |
| 4 | SUBJ_PERC_CORR | 1.400 | 1.348 | 1.450 |
| 5 | UNEM_RA | 1.261 | 1.204 | 1.319 |
| 6 | SR_TOT_RA_T | 1.183 | 1.151 | 1.208 |
| 7 | YOU_DEP_RA_T | 1.155 | 1.127 | 1.188 |
| 8 | KID_WOM_RA_T | 1.153 | 1.126 | 1.184 |
| 9 | AIR_POL | 1.142 | 1.111 | 1.169 |
| 10 | EDU38_SH | 1.135 | 1.112 | 1.159 |
| 11 | POP_80M_SH_T | 1.120 | 1.097 | 1.144 |
| 12 | DEATH_RA_M | 1.118 | 1.093 | 1.142 |
| 13 | ROOMS_PC | 1.110 | 1.095 | 1.129 |
| 14 | POP_DEN_GR_T | 1.104 | 1.073 | 1.133 |
| 15 | VOTERS_SH | 1.103 | 1.081 | 1.125 |
| 16 | HOMIC_RA | 1.097 | 1.079 | 1.120 |
| 17 | BB_ACC | 1.084 | 1.073 | 1.099 |
| 18 | POP_TOT_GI_T | 1.082 | 1.060 | 1.104 |
| 19 | DEATH_RA_T | 1.067 | 1.056 | 1.076 |
| 20 | SURF_T | 1.063 | 1.046 | 1.076 |
Accumulated Local Effects show the functional form of each predictor. Unlike partial dependence, ALE only evaluates the model on feature combinations that occur in the data — which matters here, since the selected predictors are substantially correlated.
from analysis.ale import ale_1d, plot_ale_1d
order = list(imp_rf.feature)
fig, axes = plt.subplots(4, 4, figsize=(12, 10))
for ax, feature in zip(axes.ravel(), order):
plot_ale_1d(ale_1d(models["rf"], X_all, feature), ax=ax)
fig.supylabel("ALE of subjective well-being")
plt.tight_layout()
from analysis.interactions import interaction_strength, interaction_with, plot_interaction
int_rf = interaction_strength(models["rf"], X_all)
int_xgb = interaction_strength(models["xgb"], X_all)
fig, axes = plt.subplots(1, 2, figsize=(13, 6))
plot_interaction(int_rf, ax=axes[0], title="Random Forest")
plot_interaction(int_xgb, ax=axes[1], title="XGBoost")
plt.tight_layout()
from analysis.ale import ale_2d, plot_ale_2d
emp_rf = interaction_with(models["rf"], X_all, "EMP_RA")
emp_xgb = interaction_with(models["xgb"], X_all, "EMP_RA")
r = ale_2d(models["rf"], X_all, ("EMP_RA", "SR_ELD_RA_T"), grid_size=20)
print("empty:", (r.counts == 0).sum(), "of", r.counts.size)
fig, axes = plt.subplots(1, 2, figsize=(13, 6))
plot_interaction(emp_rf, ax=axes[0], title="Random Forest",
xlabel="Interaction strength with employment rate")
plot_interaction(emp_xgb, ax=axes[1], title="XGBoost",
xlabel="Interaction strength with employment rate")
plt.tight_layout()empty: 166 of 400
from analysis.ale import ale_2d, plot_ale_2d
top_partner = emp_rf.iloc[0].feature
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
plot_ale_2d(ale_2d(models["rf"], X_all, ("EMP_RA", top_partner), grid_size=20), ax=axes[0])
axes[0].set_title(f"Strongest interaction: {top_partner}")
plot_ale_2d(ale_2d(models["rf"], X_all, ("EMP_RA", "SR_ELD_RA_T"), grid_size=20), ax=axes[1])
axes[1].set_title("Manuscript pair: SR_ELD_RA_T")
plt.tight_layout()
fig = plt.figure(figsize=(12, 12))
gs = fig.add_gridspec(2, 2, height_ratios=[1, 1.2], hspace=0.3, wspace=0.4)
ax_a = fig.add_subplot(gs[0, 0])
ax_b = fig.add_subplot(gs[0, 1])
ax_c = fig.add_subplot(gs[1, :])
plot_interaction(int_rf, ax=ax_a)
plot_interaction(emp_rf, ax=ax_b,
xlabel="Interaction strength with employment rate")
plot_ale_2d(ale_2d(models["rf"], X_all, ("EMP_RA", "SR_ELD_RA_T"), grid_size=20), ax=ax_c)
for ax, tag in zip((ax_a, ax_b, ax_c), "ABC"):
ax.text(-0.1, 1.05, tag, transform=ax.transAxes, fontsize=14, fontweight="bold")