Session 6 · Trees and Ensembles¶

Last session the logit predicted poverty with a test AUC of 0.79, but missed many poor households.

Today's question. Can a more flexible model find the poor households the logit misses, and what do we give up in exchange?

What is new today: decision trees (DecisionTreeClassifier, Gini impurity, max_depth, min_samples_leaf, plot_tree, export_text), bagging and the bootstrap, random forests (RandomForestClassifier, out-of-bag error), gradient boosting (HistGradientBoostingClassifier), cross_val_predict, and permutation importance.

In [1]:

from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.model_selection import train_test_split, StratifiedKFold, cross_validate, cross_val_predict
from sklearn.pipeline import make_pipeline, Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, recall_score, precision_score

DATA_URL = "https://raw.githubusercontent.com/hizocar/machine_learning_egei/main/docs/data/casen_households.csv"
LOCAL = Path("../data/casen_households.csv")
df = pd.read_csv(LOCAL if LOCAL.exists() else DATA_URL)

NUM = ["head_age", "head_schooling_years", "household_size"]
HOURS = ["head_weekly_hours"]
CAT = ["region", "area", "head_sex", "head_activity", "head_job_type", "health_insurance",
       "head_nationality", "head_indigenous", "has_children", "has_elderly", "household_type",
       "housing_tenure", "dwelling_type", "overcrowding", "head_education"]
X, y = df[NUM + HOURS + CAT], df["poor"]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, stratify=y, random_state=42)
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
X_train.shape, X_test.shape

from pathlib import Path import numpy as np import pandas as pd import matplotlib.pyplot as plt from sklearn.model_selection import train_test_split, StratifiedKFold, cross_validate, cross_val_predict from sklearn.pipeline import make_pipeline, Pipeline from sklearn.compose import ColumnTransformer from sklearn.impute import SimpleImputer from sklearn.preprocessing import StandardScaler, OneHotEncoder from sklearn.linear_model import LogisticRegression from sklearn.metrics import roc_auc_score, recall_score, precision_score DATA_URL = "https://raw.githubusercontent.com/hizocar/machine_learning_egei/main/docs/data/casen_households.csv" LOCAL = Path("../data/casen_households.csv") df = pd.read_csv(LOCAL if LOCAL.exists() else DATA_URL) NUM = ["head_age", "head_schooling_years", "household_size"] HOURS = ["head_weekly_hours"] CAT = ["region", "area", "head_sex", "head_activity", "head_job_type", "health_insurance", "head_nationality", "head_indigenous", "has_children", "has_elderly", "household_type", "housing_tenure", "dwelling_type", "overcrowding", "head_education"] X, y = df[NUM + HOURS + CAT], df["poor"] X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, stratify=y, random_state=42) cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42) X_train.shape, X_test.shape

Out[1]:

((11250, 19), (3750, 19))

1. The reference: the Session 5 logit¶

Every model today is compared with this one, on the same folds and the same test set.

In [2]:

pre_logit = ColumnTransformer([
    ("num", make_pipeline(SimpleImputer(strategy="median"), StandardScaler()), NUM),
    ("hours", make_pipeline(SimpleImputer(strategy="constant", fill_value=0), StandardScaler()), HOURS),
    ("cat", make_pipeline(SimpleImputer(strategy="constant", fill_value="missing"),
                          OneHotEncoder(handle_unknown="ignore")), CAT),
])
logit = Pipeline([("pre", pre_logit), ("clf", LogisticRegression(max_iter=3000))])
r = cross_validate(logit, X_train, y_train, cv=cv, scoring="roc_auc")
print(f"Logit CV AUC: {r['test_score'].mean():.3f} ± {r['test_score'].std():.3f}")

pre_logit = ColumnTransformer([ ("num", make_pipeline(SimpleImputer(strategy="median"), StandardScaler()), NUM), ("hours", make_pipeline(SimpleImputer(strategy="constant", fill_value=0), StandardScaler()), HOURS), ("cat", make_pipeline(SimpleImputer(strategy="constant", fill_value="missing"), OneHotEncoder(handle_unknown="ignore")), CAT), ]) logit = Pipeline([("pre", pre_logit), ("clf", LogisticRegression(max_iter=3000))]) r = cross_validate(logit, X_train, y_train, cv=cv, scoring="roc_auc") print(f"Logit CV AUC: {r['test_score'].mean():.3f} ± {r['test_score'].std():.3f}")

Logit CV AUC: 0.781 ± 0.009

2. One question at a time: the Gini impurity¶

A tree splits the households into two groups with a yes/no question, then splits each group again. To choose the question, it measures how mixed each group is. With a share $p$ of poor households, the Gini impurity is

$$G = 1 - p^2 - (1-p)^2 = 2p(1-p)$$

It is 0 when a group is pure (all poor or none poor) and 0.5 when it is a 50/50 mix. A split is good if the weighted average impurity of the two groups is much lower than the impurity before it.

In [3]:

def gini(p):
    return 2 * p * (1 - p)

def split_impurity(x, y, cut):
    """Weighted Gini of the two groups x <= cut and x > cut."""
    left, right = y[x <= cut], y[x > cut]
    n = len(y)
    return len(left) / n * gini(left.mean()) + len(right) / n * gini(right.mean())

d = pd.DataFrame({"x": X_train["head_schooling_years"], "y": y_train}).dropna()
print(f"Before splitting: poor share {d.y.mean():.3f}, Gini {gini(d.y.mean()):.3f}")
cuts = np.arange(0.5, 20, 1.0)
imp = [split_impurity(d.x, d.y, c) for c in cuts]
best = cuts[int(np.argmin(imp))]
print(f"Best cut on schooling: <= {best}  ->  weighted Gini {min(imp):.3f}")

plt.plot(cuts, imp, marker="o")
plt.axhline(gini(d.y.mean()), ls="--", c="grey", label="before the split")
plt.xlabel("cut point: years of schooling"); plt.ylabel("weighted Gini after the split"); plt.legend();

def gini(p): return 2 * p * (1 - p) def split_impurity(x, y, cut): """Weighted Gini of the two groups x <= cut and x > cut.""" left, right = y[x <= cut], y[x > cut] n = len(y) return len(left) / n * gini(left.mean()) + len(right) / n * gini(right.mean()) d = pd.DataFrame({"x": X_train["head_schooling_years"], "y": y_train}).dropna() print(f"Before splitting: poor share {d.y.mean():.3f}, Gini {gini(d.y.mean()):.3f}") cuts = np.arange(0.5, 20, 1.0) imp = [split_impurity(d.x, d.y, c) for c in cuts] best = cuts[int(np.argmin(imp))] print(f"Best cut on schooling: <= {best} -> weighted Gini {min(imp):.3f}") plt.plot(cuts, imp, marker="o") plt.axhline(gini(d.y.mean()), ls="--", c="grey", label="before the split") plt.xlabel("cut point: years of schooling"); plt.ylabel("weighted Gini after the split"); plt.legend();

Before splitting: poor share 0.231, Gini 0.355
Best cut on schooling: <= 12.5  ->  weighted Gini 0.346

The tree does exactly this for every variable and every cut point, and keeps the best one. Then it repeats inside each group. No scaling is needed: a cut "schooling ≤ 8.5" is the same question whether schooling is measured in years or in standardised units.

3. A first tree¶

Trees need numbers too, so we still one-hot encode the categories, but we drop the scaler.

In [4]:

from sklearn.tree import DecisionTreeClassifier, plot_tree, export_text

pre_tree = ColumnTransformer([
    ("num", SimpleImputer(strategy="median"), NUM),
    ("hours", SimpleImputer(strategy="constant", fill_value=0), HOURS),
    ("cat", make_pipeline(SimpleImputer(strategy="constant", fill_value="missing"),
                          OneHotEncoder(handle_unknown="ignore", sparse_output=False)), CAT),
])
tree3 = Pipeline([("pre", pre_tree), ("clf", DecisionTreeClassifier(max_depth=3, random_state=42))])
tree3.fit(X_train, y_train)
names = [n.split("__", 1)[1] for n in tree3.named_steps["pre"].get_feature_names_out()]

plt.figure(figsize=(16, 6))
plot_tree(tree3.named_steps["clf"], feature_names=names, class_names=["not poor", "poor"],
          filled=True, impurity=False, proportion=True, rounded=True, fontsize=8);

from sklearn.tree import DecisionTreeClassifier, plot_tree, export_text pre_tree = ColumnTransformer([ ("num", SimpleImputer(strategy="median"), NUM), ("hours", SimpleImputer(strategy="constant", fill_value=0), HOURS), ("cat", make_pipeline(SimpleImputer(strategy="constant", fill_value="missing"), OneHotEncoder(handle_unknown="ignore", sparse_output=False)), CAT), ]) tree3 = Pipeline([("pre", pre_tree), ("clf", DecisionTreeClassifier(max_depth=3, random_state=42))]) tree3.fit(X_train, y_train) names = [n.split("__", 1)[1] for n in tree3.named_steps["pre"].get_feature_names_out()] plt.figure(figsize=(16, 6)) plot_tree(tree3.named_steps["clf"], feature_names=names, class_names=["not poor", "poor"], filled=True, impurity=False, proportion=True, rounded=True, fontsize=8);

The first question is not about schooling: it is "does the head work 27.5 hours a week or less?" (which includes heads who do not work at all). Among the households answering yes, a third are poor; among the rest, 13%. Read each path from the top: it is a set of rules anyone can apply.

In [5]:

print(export_text(tree3.named_steps["clf"], feature_names=names, show_weights=False))

print(export_text(tree3.named_steps["clf"], feature_names=names, show_weights=False))

|--- head_weekly_hours <= 27.50
|   |--- head_age <= 64.50
|   |   |--- housing_tenure_owned <= 0.50
|   |   |   |--- class: 1
|   |   |--- housing_tenure_owned >  0.50
|   |   |   |--- class: 0
|   |--- head_age >  64.50
|   |   |--- household_size <= 1.50
|   |   |   |--- class: 0
|   |   |--- household_size >  1.50
|   |   |   |--- class: 0
|--- head_weekly_hours >  27.50
|   |--- has_children_yes <= 0.50
|   |   |--- head_job_type_self_employed <= 0.50
|   |   |   |--- class: 0
|   |   |--- head_job_type_self_employed >  0.50
|   |   |   |--- class: 0
|   |--- has_children_yes >  0.50
|   |   |--- head_schooling_years <= 13.50
|   |   |   |--- class: 0
|   |   |--- head_schooling_years >  13.50
|   |   |   |--- class: 0

4. How deep? Overfitting again¶

Each extra level doubles the number of groups. Deep enough, every leaf holds a handful of households and the tree memorises the training set, like KNN with k = 1.

In [6]:

rows = []
for depth in [1, 2, 3, 4, 5, 6, 8, 10, 12, None]:
    m = Pipeline([("pre", pre_tree), ("clf", DecisionTreeClassifier(max_depth=depth, random_state=42))])
    s = cross_validate(m, X_train, y_train, cv=cv, scoring="roc_auc", return_train_score=True)
    rows.append({"max_depth": "none" if depth is None else depth,
                 "train AUC": s["train_score"].mean(), "CV AUC": s["test_score"].mean()})
depth_table = pd.DataFrame(rows).round(3)
depth_table

rows = [] for depth in [1, 2, 3, 4, 5, 6, 8, 10, 12, None]: m = Pipeline([("pre", pre_tree), ("clf", DecisionTreeClassifier(max_depth=depth, random_state=42))]) s = cross_validate(m, X_train, y_train, cv=cv, scoring="roc_auc", return_train_score=True) rows.append({"max_depth": "none" if depth is None else depth, "train AUC": s["train_score"].mean(), "CV AUC": s["test_score"].mean()}) depth_table = pd.DataFrame(rows).round(3) depth_table

Out[6]:

	max_depth	train AUC	CV AUC
0	1	0.642	0.642
1	2	0.699	0.696
2	3	0.738	0.735
3	4	0.757	0.746
4	5	0.776	0.753
5	6	0.795	0.759
6	8	0.837	0.746
7	10	0.885	0.724
8	12	0.929	0.693
9	none	1.000	0.618

In [7]:

dt = depth_table[depth_table["max_depth"] != "none"]
plt.plot(dt["max_depth"], dt["train AUC"], marker="o", label="train")
plt.plot(dt["max_depth"], dt["CV AUC"], marker="o", label="cross-validation")
plt.axhline(r["test_score"].mean(), ls="--", c="grey", label="logit (CV)")
plt.xlabel("max_depth"); plt.ylabel("AUC"); plt.legend(); plt.title("A deeper tree is not a better tree");

dt = depth_table[depth_table["max_depth"] != "none"] plt.plot(dt["max_depth"], dt["train AUC"], marker="o", label="train") plt.plot(dt["max_depth"], dt["CV AUC"], marker="o", label="cross-validation") plt.axhline(r["test_score"].mean(), ls="--", c="grey", label="logit (CV)") plt.xlabel("max_depth"); plt.ylabel("AUC"); plt.legend(); plt.title("A deeper tree is not a better tree");

Training AUC climbs to 1.0 for a tree with no depth limit; its CV AUC collapses to about 0.62. The best depth by CV is around 6, and even then a single tree stays below the logit.

min_samples_leaf is the other brake: it forbids leaves with fewer than that many households.

In [8]:

for leaf in [1, 20, 100]:
    m = Pipeline([("pre", pre_tree), ("clf", DecisionTreeClassifier(min_samples_leaf=leaf, random_state=42))])
    s = cross_validate(m, X_train, y_train, cv=cv, scoring="roc_auc")
    print(f"no depth limit, min_samples_leaf={leaf:<3}  CV AUC {s['test_score'].mean():.3f}")

for leaf in [1, 20, 100]: m = Pipeline([("pre", pre_tree), ("clf", DecisionTreeClassifier(min_samples_leaf=leaf, random_state=42))]) s = cross_validate(m, X_train, y_train, cv=cv, scoring="roc_auc") print(f"no depth limit, min_samples_leaf={leaf:<3} CV AUC {s['test_score'].mean():.3f}")

no depth limit, min_samples_leaf=1    CV AUC 0.618

no depth limit, min_samples_leaf=20   CV AUC 0.749

no depth limit, min_samples_leaf=100  CV AUC 0.763

5. Many trees: bagging and random forests¶

A single deep tree is unstable: change a few households and the questions change. The fix is to grow many trees and average them.

1. Bootstrap: draw a sample of the same size with replacement from the training set (some households appear twice, about a third not at all).
2. Bagging (bootstrap aggregating): fit one tree per bootstrap sample and average their predicted probabilities. Each tree overfits differently; the average cancels much of the noise.
3. Random forest: bagging plus one more trick. At each split, a tree may only look at a random subset of the variables (max_features), so the trees disagree more, and the average improves.

The households a tree never saw (its out-of-bag sample) give a free validation score.

In [9]:

from sklearn.ensemble import RandomForestClassifier

bagging = Pipeline([("pre", pre_tree), ("clf", RandomForestClassifier(
    n_estimators=300, max_features=None, min_samples_leaf=5, n_jobs=-1, random_state=42))])
forest = Pipeline([("pre", pre_tree), ("clf", RandomForestClassifier(
    n_estimators=300, min_samples_leaf=5, oob_score=True, n_jobs=-1, random_state=42))])

for name, m in [("bagging (all variables at each split)", bagging), ("random forest (sqrt of the variables)", forest)]:
    s = cross_validate(m, X_train, y_train, cv=cv, scoring="roc_auc")
    print(f"{name:40s} CV AUC {s['test_score'].mean():.3f} ± {s['test_score'].std():.3f}")

forest.fit(X_train, y_train)
oob_p = forest.named_steps["clf"].oob_decision_function_[:, 1]
print(f"Out-of-bag AUC of the forest (no CV needed): {roc_auc_score(y_train, oob_p):.3f}")

from sklearn.ensemble import RandomForestClassifier bagging = Pipeline([("pre", pre_tree), ("clf", RandomForestClassifier( n_estimators=300, max_features=None, min_samples_leaf=5, n_jobs=-1, random_state=42))]) forest = Pipeline([("pre", pre_tree), ("clf", RandomForestClassifier( n_estimators=300, min_samples_leaf=5, oob_score=True, n_jobs=-1, random_state=42))]) for name, m in [("bagging (all variables at each split)", bagging), ("random forest (sqrt of the variables)", forest)]: s = cross_validate(m, X_train, y_train, cv=cv, scoring="roc_auc") print(f"{name:40s} CV AUC {s['test_score'].mean():.3f} ± {s['test_score'].std():.3f}") forest.fit(X_train, y_train) oob_p = forest.named_steps["clf"].oob_decision_function_[:, 1] print(f"Out-of-bag AUC of the forest (no CV needed): {roc_auc_score(y_train, oob_p):.3f}")

bagging (all variables at each split)    CV AUC 0.787 ± 0.004

random forest (sqrt of the variables)    CV AUC 0.797 ± 0.007

Out-of-bag AUC of the forest (no CV needed): 0.798

6. Trees that correct each other: gradient boosting¶

Boosting grows trees in sequence. Each new, small tree is fitted to what the current model still gets wrong, and is added with a small weight (learning_rate). After many rounds, many weak trees add up to a strong model. Too many rounds overfit, so HistGradientBoostingClassifier holds out 10% of the training data and stops early when that held-out score stops improving.

In [10]:

from sklearn.ensemble import HistGradientBoostingClassifier

demo = Pipeline([("pre", pre_tree), ("clf", HistGradientBoostingClassifier(
    early_stopping=True, scoring="roc_auc", max_iter=500, random_state=42))])
demo.fit(X_train, y_train)
gb = demo.named_steps["clf"]
plt.plot(gb.train_score_, label="training data")
plt.plot(gb.validation_score_, label="held-out 10%")
plt.xlabel("boosting rounds (trees added)"); plt.ylabel("AUC"); plt.legend()
plt.title(f"Early stopping after {gb.n_iter_} rounds");

from sklearn.ensemble import HistGradientBoostingClassifier demo = Pipeline([("pre", pre_tree), ("clf", HistGradientBoostingClassifier( early_stopping=True, scoring="roc_auc", max_iter=500, random_state=42))]) demo.fit(X_train, y_train) gb = demo.named_steps["clf"] plt.plot(gb.train_score_, label="training data") plt.plot(gb.validation_score_, label="held-out 10%") plt.xlabel("boosting rounds (trees added)"); plt.ylabel("AUC"); plt.legend() plt.title(f"Early stopping after {gb.n_iter_} rounds");

7. The race, done properly¶

To compare models fairly we need the same steps for every model, so we write them once, as a function (that is today's assistant skill). The threshold is chosen with cross_val_predict: out-of-fold probabilities on the training set, never the test set.

In [11]:

def cost_threshold(p, y_true, fn_cost=3, fp_cost=1):
    """Threshold on a 0.01 grid that minimises fn_cost*FN + fp_cost*FP."""
    ts = np.round(np.arange(0.02, 0.90, 0.01), 2)
    cost = [fn_cost * ((p < t) & (y_true == 1)).sum() + fp_cost * ((p >= t) & (y_true == 0)).sum() for t in ts]
    return ts[int(np.argmin(cost))]


def evaluate(model, X_train, y_train, X_test, y_test, cv):
    """CV AUC, a cost-based threshold chosen on training folds, and test-set results."""
    s = cross_validate(model, X_train, y_train, cv=cv, scoring="roc_auc", return_train_score=True)
    oof = cross_val_predict(model, X_train, y_train, cv=cv, method="predict_proba")[:, 1]
    t = cost_threshold(oof, y_train.to_numpy())
    model.fit(X_train, y_train)
    p = model.predict_proba(X_test)[:, 1]
    yt = y_test.to_numpy()
    return {"train AUC": s["train_score"].mean(), "CV AUC": s["test_score"].mean(),
            "CV sd": s["test_score"].std(), "test AUC": roc_auc_score(yt, p), "threshold": t,
            "recall": recall_score(yt, p >= t), "precision": precision_score(yt, p >= t),
            "cost 3FN+FP": 3 * ((p < t) & (yt == 1)).sum() + ((p >= t) & (yt == 0)).sum()}


models = {
    "Logit (Session 5)": logit,
    "Tree (depth 4)": Pipeline([("pre", pre_tree), ("clf", DecisionTreeClassifier(max_depth=4, random_state=42))]),
    "Random forest": Pipeline([("pre", pre_tree), ("clf", RandomForestClassifier(
        n_estimators=300, min_samples_leaf=5, n_jobs=-1, random_state=42))]),
    "Gradient boosting": Pipeline([("pre", pre_tree), ("clf", HistGradientBoostingClassifier(random_state=42))]),
}
results = pd.DataFrame({name: evaluate(m, X_train, y_train, X_test, y_test, cv) for name, m in models.items()}).T
results.round(3)

def cost_threshold(p, y_true, fn_cost=3, fp_cost=1): """Threshold on a 0.01 grid that minimises fn_cost*FN + fp_cost*FP.""" ts = np.round(np.arange(0.02, 0.90, 0.01), 2) cost = [fn_cost * ((p < t) & (y_true == 1)).sum() + fp_cost * ((p >= t) & (y_true == 0)).sum() for t in ts] return ts[int(np.argmin(cost))] def evaluate(model, X_train, y_train, X_test, y_test, cv): """CV AUC, a cost-based threshold chosen on training folds, and test-set results.""" s = cross_validate(model, X_train, y_train, cv=cv, scoring="roc_auc", return_train_score=True) oof = cross_val_predict(model, X_train, y_train, cv=cv, method="predict_proba")[:, 1] t = cost_threshold(oof, y_train.to_numpy()) model.fit(X_train, y_train) p = model.predict_proba(X_test)[:, 1] yt = y_test.to_numpy() return {"train AUC": s["train_score"].mean(), "CV AUC": s["test_score"].mean(), "CV sd": s["test_score"].std(), "test AUC": roc_auc_score(yt, p), "threshold": t, "recall": recall_score(yt, p >= t), "precision": precision_score(yt, p >= t), "cost 3FN+FP": 3 * ((p < t) & (yt == 1)).sum() + ((p >= t) & (yt == 0)).sum()} models = { "Logit (Session 5)": logit, "Tree (depth 4)": Pipeline([("pre", pre_tree), ("clf", DecisionTreeClassifier(max_depth=4, random_state=42))]), "Random forest": Pipeline([("pre", pre_tree), ("clf", RandomForestClassifier( n_estimators=300, min_samples_leaf=5, n_jobs=-1, random_state=42))]), "Gradient boosting": Pipeline([("pre", pre_tree), ("clf", HistGradientBoostingClassifier(random_state=42))]), } results = pd.DataFrame({name: evaluate(m, X_train, y_train, X_test, y_test, cv) for name, m in models.items()}).T results.round(3)

Out[11]:

	train AUC	CV AUC	CV sd	test AUC	threshold	recall	precision	cost 3FN+FP
Logit (Session 5)	0.789	0.781	0.009	0.794	0.26	0.680	0.454	1540.0
Tree (depth 4)	0.757	0.746	0.009	0.767	0.24	0.682	0.446	1560.0
Random forest	0.903	0.797	0.007	0.814	0.24	0.745	0.460	1421.0
Gradient boosting	0.910	0.798	0.007	0.815	0.20	0.788	0.439	1422.0

The ensembles gain about 0.02 AUC over the logit, and at their own cost-minimising thresholds the policy cost falls by roughly 8%. That is real, but it is not a revolution. The depth-4 tree is slightly worse than the logit, but anyone can read it.

8. What does the forest rely on? Permutation importance¶

Impurity importance (feature_importances_) adds up how much each one-hot column reduced the Gini. It is computed on training data and favours variables with many possible cut points. Permutation importance asks a cleaner question on the test set: if I shuffle this column, so it carries no information, how much AUC do we lose?

In [12]:

from sklearn.inspection import permutation_importance

rf = models["Random forest"]
perm = permutation_importance(rf, X_test, y_test, scoring="roc_auc", n_repeats=5, random_state=42, n_jobs=-1)
perm_imp = pd.Series(perm.importances_mean, index=X_test.columns).sort_values(ascending=False)
impurity = pd.Series(rf.named_steps["clf"].feature_importances_, index=names).sort_values(ascending=False)

fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
perm_imp.head(10)[::-1].plot.barh(ax=ax[0], title="Permutation importance (AUC lost, test set)")
impurity.head(10)[::-1].plot.barh(ax=ax[1], title="Impurity importance (one-hot columns, train)")
plt.tight_layout()
perm_imp.head(8).round(4)

from sklearn.inspection import permutation_importance rf = models["Random forest"] perm = permutation_importance(rf, X_test, y_test, scoring="roc_auc", n_repeats=5, random_state=42, n_jobs=-1) perm_imp = pd.Series(perm.importances_mean, index=X_test.columns).sort_values(ascending=False) impurity = pd.Series(rf.named_steps["clf"].feature_importances_, index=names).sort_values(ascending=False) fig, ax = plt.subplots(1, 2, figsize=(13, 4.5)) perm_imp.head(10)[::-1].plot.barh(ax=ax[0], title="Permutation importance (AUC lost, test set)") impurity.head(10)[::-1].plot.barh(ax=ax[1], title="Impurity importance (one-hot columns, train)") plt.tight_layout() perm_imp.head(8).round(4)

Out[12]:

head_age             0.0373
head_weekly_hours    0.0310
has_children         0.0291
housing_tenure       0.0284
household_type       0.0189
has_elderly          0.0175
head_activity        0.0160
head_job_type        0.0137
dtype: float64

Importance says what the model uses to predict, not what causes poverty. Working fewer hours is associated with poverty, partly because poverty and unemployment share causes. And when two variables carry the same information (schooling years and education level), shuffling one costs little because the other covers for it.

9. Accuracy vs. explanation¶

A household excluded from a programme can ask why. With the depth-3 tree the answer is a short chain of questions. With 300 trees there is no short answer. Whether 0.02 of AUC is worth that is a decision for the programme, and your job is to put both numbers on the table.

10. With your assistant: refactor a script into functions¶

Here is my notebook code that evaluates four classifiers (paste sections 1-7).
Refactor it into two functions: cost_threshold(p, y, fn_cost, fp_cost) and
evaluate(model, X_train, y_train, X_test, y_test, cv) -> dict.
Keep the results identical. The threshold must be chosen with out-of-fold predictions on the
training set only. Add a docstring to each function and nothing else.

Check: run the old and the new code and compare the table number by number. A refactor that changes a result is a bug, even if the new number looks better.

11. Exercises¶

1. Change the cost ratio to 5:1 in cost_threshold. Which model wins now, and how many more poor households does it reach?
2. Fit a random forest with n_estimators = 10, 50, 300. How do the CV AUC and the run time change?
3. Write the depth-3 tree as a one-page "scorecard" a municipal officer could apply with pen and paper. Which household in the test set would be hardest to explain?