6.8 Explaining models

A flexible model can be probed even when it cannot be read: shuffle a predictor and see how much worse it gets, vary a predictor and watch the prediction, or split a single prediction fairly among the predictors. These explain the model, which is not the same as explaining the world.

import itertools
import math

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.axes import Axes
from sklearn import datasets
from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor
from sklearn.inspection import permutation_importance
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split

from statly import style

style.use()
rng = np.random.default_rng(1)

diabetes = datasets.load_diabetes(
    as_frame=True, scaled=False
)  # 442 patients, disease progression
X, y = diabetes.data, diabetes.target.to_numpy()
NAMES = list(X.columns)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
boosted = GradientBoostingRegressor(
    n_estimators=200, max_depth=3, learning_rate=0.03, random_state=0
).fit(X_train, y_train)

6.8.1 Which predictors matter? Permutation importance

A model’s own importance scores (for trees, how much each split reduced the error, Chapter 6.6) are measured on the training data and favour predictors with many possible split points. Permutation importance asks a cleaner question on held-out data: if the values of one predictor were shuffled among the patients, breaking its link with the outcome while keeping everything else, how much worse would the predictions be?

result = permutation_importance(
    boosted, X_test, y_test, scoring="neg_mean_squared_error", n_repeats=30, random_state=0
)
order = np.argsort(result.importances_mean)

fig, ax = plt.subplots(figsize=(8, 3.6))
ax.barh(np.array(NAMES)[order], result.importances_mean[order], color=style.BLUE, alpha=0.7)
ax.errorbar(
    result.importances_mean[order],
    range(len(NAMES)),
    xerr=[
        result.importances_mean[order] - result.importances.min(axis=1)[order],
        result.importances.max(axis=1)[order] - result.importances_mean[order],
    ],
    fmt="none",
    color=style.INK,
    lw=1,
)
ax.axvline(0, color=style.INK, lw=0.8)
ax.set(xlabel="increase in mean squared error when shuffled")
plt.show()

Figure 6.8.1: Permutation importance for boosted trees predicting disease progression, on 133 held-out patients: the increase in mean squared error when each predictor is shuffled (average and range over 30 shuffles). s5 (a blood serum measurement) and bmi matter most; several predictors could be shuffled with almost no loss.

A caution: correlated predictors share credit. If two predictors carry the same information, shuffling either one alone hurts little, because the model can lean on the other, and both look unimportant.

X_copy_train, X_copy_test = (
    X_train.assign(bmi_copy=X_train.bmi),
    X_test.assign(bmi_copy=X_test.bmi),
)
with_copy = GradientBoostingRegressor(
    n_estimators=200, max_depth=3, learning_rate=0.03, random_state=0, max_features=0.5
).fit(X_copy_train, y_train)
copy_result = permutation_importance(
    with_copy,
    X_copy_test,
    y_test,
    scoring="neg_mean_squared_error",
    n_repeats=30,
    random_state=0,
)
alone = result.importances_mean[NAMES.index("bmi")]
shared = copy_result.importances_mean[
    [list(X_copy_test.columns).index(c) for c in ("bmi", "bmi_copy")]
]

fig, ax = plt.subplots(figsize=(7, 2.4))
ax.barh(
    ["bmi alone", "bmi (with copy)", "bmi copy"],
    [alone, *shared],
    color=[style.BLUE, style.GOLD, style.GOLD],
    alpha=0.75,
)
ax.invert_yaxis()
ax.set(xlabel="permutation importance")
plt.show()

Figure 6.8.2: The same model refitted with an exact copy of bmi added as an extra column. The importance that bmi had alone (left) is now split between the two copies (right), and each looks less important than it really is.

6.8.2 How does a predictor act? Partial dependence and ICE

Partial dependence shows the average prediction as one predictor is moved across its range, with every patient’s other values kept as they are. Individual conditional expectation (ICE) curves show the same for each patient separately: where they run parallel, the predictor acts the same way for everyone; where they cross, it interacts with something else.

def ice_curves(feature: str, patients: int = 60) -> tuple[np.ndarray, np.ndarray]:
    """Predictions for a sample of patients as one predictor moves over a grid of values."""
    grid = np.linspace(*np.percentile(X[feature], [2, 98]), 40)
    sample = X_test.sample(patients, random_state=0)
    curves = np.empty((patients, len(grid)))
    for k, value in enumerate(grid):
        curves[:, k] = boosted.predict(sample.assign(**{feature: value}))
    return grid, curves


fig, axes = style.subplots(2, height=3.2, width=10)
for ax, feature in zip(axes, ("bmi", "s5")):
    grid, curves = ice_curves(feature)
    ax.plot(grid, curves.T, color=style.GREY, lw=0.6, alpha=0.6)
    ax.plot(grid, curves.mean(axis=0), color=style.CORAL, lw=3)
    ax.set(xlabel=feature, ylabel="predicted progression")
fig.tight_layout()
plt.show()

Figure 6.8.3: Partial dependence (thick red line) and ICE curves for 60 patients (grey) for bmi and s5 in the boosted model. Both raise the predicted progression, with steps where the trees split. The ICE curves are roughly parallel: in this model, each predictor shifts everyone’s prediction by about the same amount.

6.8.3 Why this prediction? Shapley values

For one patient, the model predicts a progression of, say, 230, against an average prediction of 150. Which predictors account for the difference of 80? Shapley values, from cooperative game theory, split it fairly. Imagine adding the predictors one at a time, in random order, and record how much each one changes the prediction when it joins; a predictor’s Shapley value is its average contribution over all orders:

$$\phi_j = \sum_{S \subseteq F \setminus \{j\}} \frac{|S|!\,(|F| - |S| - 1)!}{|F|!} \Big( f(S \cup \{j\}) - f(S) \Big)$$

Here $f(S)$ is the prediction when only the predictors in S are known (the others averaged over the data). The values add up exactly: average prediction plus all $\phi_j$ equals the patient’s prediction. For a linear model, $\phi_j = \beta_j (x_j - \bar x_j)$.

background = X_train.sample(100, random_state=1)
patient = X_test.iloc[[5]]


def value_of(subset: tuple[int, ...]) -> float:
    """Average prediction when the predictors in `subset` take the patient's values."""
    data = background.copy()
    for j in subset:
        data.iloc[:, j] = patient.iloc[0, j]
    return float(boosted.predict(data).mean())


values = {s: value_of(s) for r in range(11) for s in itertools.combinations(range(10), r)}
shapley = np.zeros(10)
for j in range(10):
    for subset, v in values.items():
        if j in subset:
            continue
        weight = (
            math.factorial(len(subset))
            * math.factorial(10 - len(subset) - 1)
            / math.factorial(10)
        )
        shapley[j] += weight * (values[tuple(sorted(subset + (j,)))] - v)
base, final = values[()], values[tuple(range(10))]

fig, ax = plt.subplots(figsize=(8, 3.8))
order = np.argsort(np.abs(shapley))
running = base
for row, j in enumerate(order):
    color = style.CORAL if shapley[j] > 0 else style.BLUE
    ax.barh(row, shapley[j], left=running, color=color, alpha=0.8)
    ax.text(
        max(running, running + shapley[j]) + 1,
        row,
        f"{shapley[j]:+.0f}  ({NAMES[j]} = {patient.iloc[0, j]:.4g})",
        va="center",
        fontsize=8,
    )
    running += shapley[j]
ax.axvline(base, color=style.INK, ls="--", lw=1)
ax.axvline(final, color=style.INK, lw=1)
ax.set_yticks(range(10), np.array(NAMES)[order])
ax.set(
    xlabel=f"prediction: from the average {base:.0f} to this patient's {final:.0f}",
    xlim=(min(base, final) - 30, max(base, final) + 70),
)
plt.show()

Figure 6.8.4: Exact Shapley values for one held-out patient, computed over all 1,024 subsets of the ten predictors (other predictors averaged over 100 background patients). Starting from the average prediction, each bar adds one predictor’s contribution, red raising and blue lowering the prediction, ending at this patient’s prediction. High s5 and bmi push it up; the rest matter little.

Computing all subsets is only possible with few predictors. With many, Shapley values are estimated by sampling random orders, and for tree models there are fast exact algorithms. Plotting the Shapley values of many patients gives a global picture that also shows direction:

def sampled_shapley(row: np.ndarray, orders: int = 60) -> np.ndarray:
    """Average contribution of each predictor over random orders of adding them."""
    contributions = np.zeros(10)
    start = background.to_numpy()[:30]
    for _ in range(orders):
        data = start.copy()
        previous = boosted.predict(data).mean()
        for j in rng.permutation(10):
            data[:, j] = row[j]
            current = boosted.predict(data).mean()
            contributions[j] += current - previous
            previous = current
    return contributions / orders


explained = X_test.iloc[:150]
phis = np.array([sampled_shapley(r) for r in explained.to_numpy()])
order = np.argsort(np.abs(phis).mean(axis=0))

fig, ax = plt.subplots(figsize=(8, 4))
for row, j in enumerate(order):
    values_j = explained.iloc[:, j].to_numpy()
    scaled = (values_j - values_j.min()) / (np.ptp(values_j) + 1e-9)
    ax.scatter(
        phis[:, j], row + rng.uniform(-0.25, 0.25, len(phis)), c=scaled, cmap="coolwarm", s=8
    )
ax.axvline(0, color=style.INK, lw=0.8)
ax.set_yticks(range(10), np.array(NAMES)[order])
ax.set(xlabel="Shapley value (contribution to the prediction)")
plt.show()

Figure 6.8.5: Estimated Shapley values (60 random orders each) for 150 held-out patients, one dot per patient and predictor, coloured by the predictor’s value (blue low, red high). For s5 and bmi, high values push predictions up and low values down; for most other predictors, the dots stay near zero.

6.8.4 Fairness

When a model’s decisions affect people (loans, hiring, medical triage), it matters how its errors fall across groups. A model can be equally accurate overall and still make one kind of mistake far more often for one group. There are several definitions of fairness, and they conflict:

• Calibration within groups: a score of 30% means 30% in every group.
• Equal error rates: the same false positive rate and the same false negative rate in every group.
• Demographic parity: the same share of positive decisions in every group.

When the groups have different base rates of the outcome, a calibrated score cannot have equal error rates in both groups (unless it is perfect). Choosing among the definitions is a value judgement, not a statistical one; what statistics can do is make the trade-offs visible.

n = 20_000
group = rng.uniform(size=n) < 0.5
true_logit = np.where(group, -1.1, -2.2) + rng.normal(0, 1, n)
default = rng.uniform(size=n) < 1 / (1 + np.exp(-true_logit))
signal = true_logit + rng.normal(0, 0.8, n)
score_model = LogisticRegression().fit(np.column_stack([signal, group]), default)
risk = score_model.predict_proba(np.column_stack([signal, group]))[:, 1]
refused = risk > 0.2


def calibration(ax: Axes) -> None:
    """Observed default rate in bins of predicted risk, per group."""
    bins = np.linspace(0, 0.8, 9)
    for g, color, name in ((False, style.BLUE, "group A"), (True, style.CORAL, "group B")):
        mask = group == g
        index = np.digitize(risk[mask], bins)
        centres = [
            risk[mask][index == b].mean()
            for b in range(1, len(bins))
            if (index == b).sum() > 50
        ]
        observed = [
            default[mask][index == b].mean()
            for b in range(1, len(bins))
            if (index == b).sum() > 50
        ]
        ax.plot(centres, observed, "o-", color=color, lw=2, label=name)
    ax.plot([0, 0.8], [0, 0.8], color=style.INK, ls="--", lw=1)
    ax.xaxis.set_major_formatter(style.PERCENT)
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.legend(fontsize=8)
    ax.set(
        xlabel="predicted risk",
        ylabel="observed default rate",
        title="Calibrated in both groups",
    )


def error_rates(ax: Axes) -> None:
    """False positive rate, false negative rate and refusal share per group."""
    metrics = {}
    for g, name in ((False, "group A"), (True, "group B")):
        mask = group == g
        metrics[name] = [
            refused[mask & ~default].mean(),
            (~refused)[mask & default].mean(),
            refused[mask].mean(),
        ]
    labels = ["false positive rate", "false negative rate", "share refused"]
    x = np.arange(3)
    for offset, (name, values), color in zip(
        (-0.2, 0.2), metrics.items(), (style.BLUE, style.CORAL)
    ):
        ax.bar(x + offset, values, width=0.4, color=color, alpha=0.75, label=name)
    ax.set_xticks(x, labels, fontsize=8.5)
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.legend(fontsize=8)
    ax.set(ylabel="rate", title="Threshold 20% for everyone")


fig, axes = style.subplots(2, height=3.4, width=10.5)
calibration(axes[0])
error_rates(axes[1])
fig.tight_layout()
plt.show()

Figure 6.8.6: A simulated loan-default model applied to two groups with different default rates (10% and 25%). The score is calibrated in both (left: predicted risk matches the observed rate). With one threshold of 20% for both groups, the false positive rate (people refused who would have repaid) is much higher in group B, and the share refused differs (right).

6.8.5 What explanations do not say

• They describe the model, not the world. A predictor the model relies on is not thereby a cause; changing it in reality may not change the outcome (Chapter 5.4).
• Correlated predictors make any attribution partly arbitrary.
• Shuffling and moving one predictor at a time creates combinations that never occur (a tall person with a child’s weight), and the model’s behaviour there is meaningless.
• An explanation that looks sensible is not evidence that the model is right; check it on held-out data (Chapter 6.2).

6.8.6 In one picture

• Permutation importance: loss increase on held-out data when a predictor is shuffled; correlated predictors share credit.
• Partial dependence shows the average effect of moving a predictor; ICE curves show it per case.
• Shapley values split one prediction fairly among predictors and add up to it; for a linear model, $\beta_j(x_j - \bar x_j)$.
• Fairness criteria (calibration, equal error rates, parity) conflict when base rates differ; show the trade-offs.
• Explanations describe the model, not causes.

Next: Part 7, unsupervised learning, methods with no outcome to learn from.