Explaining models to stakeholders

A manager asks what drives readmissions. A clinician asks why this patient scored 47%. Those are different questions. The first needs a global explanation, a summary of how the model behaves across everyone. The second needs a local one, about a single prediction. Both describe the model, not the patients. Most explanation mistakes come from mixing those up, or from forgetting that correlated features make any summary ambiguous.

This article explains a synthetic readmission risk model with the tools in sklearn.inspection, shows where each one misleads, and ends with the model card I’d hand over with it.

The example model

The data is 12,000 synthetic adult discharges, generated by explain_risk_model.py with a fixed seed, with 13.0% readmitted as an emergency within 30 days. The features are age, a frailty score from 1 to 9, emergency admissions and ED attendances in the last year, comorbidity and medication counts, length of stay, living alone and weekend discharge. I planted three things that make explanation hard:

• Correlated pairs. ED attendances mostly repeat emergency admissions (correlation 0.83 across all discharges), and the medication count tracks comorbidities (0.86). In the simulation only admissions and comorbidities affect the outcome. The partners have no effect of their own.
• Age and frailty move together (correlation 0.71).
• An interaction. Living alone raises risk mainly for patients with a frailty score of 6 or more. Frail patients who live alone are 2.1% of discharges.

I trained a gradient boosting model and a logistic regression on a random 70%, and tested on the other 3,600 discharges. Both reached an AUROC of 0.696. The true probabilities that generated the outcomes reach 0.703, so neither model is leaving much on the table. Boosting’s mean predicted risk was 13.1% against 13.0% observed. A real model needs a time split and a harder look at calibration, which my 30-day readmission risk model and calibration before AUC cover. I’ll explain the boosting model, because it has no coefficients to fall back on.

Global: permutation importance

Permutation importance shuffles one column in the test set, breaking its link with the outcome, and measures how much the score falls. scikit-learn does it in one call:

result = permutation_importance(
    model, X, y, scoring="roc_auc", n_repeats=N_REPEATS, random_state=SEED
)

Data table

Feature	Gradient boosting	Gradient boosting (low)	Gradient boosting (high)
Age	0.041	0.025	0.049
Frailty score	0.034	0.030	0.044
Emergency admissions, 12 months	0.025	0.019	0.034
Comorbidities	0.009	0.002	0.015
Length of stay	0.004	-0.001	0.009
ED attendances, 12 months	0.002	0.001	0.004
Lives alone	0.001	-0.001	0.003
Weekend discharge	0.000	-0.000	0.001
Medications	-0.001	-0.003	-0.000

Synthetic data. Source: projects/article-examples/round-two/explain_risk_model.py.

Age, frailty and prior admissions carry most of the model. Two readings are wrong, and both are tempting.

“Living alone doesn’t matter.” Its importance is 0.001. But for the frail patients it applies to, it matters a lot: for Patient A below, it’s the difference between 47% and 30%. Permutation importance averages over everyone, and a strong effect in 2% of patients averages out to almost nothing. Global importance tells you what the model leans on across the population, not what matters for a given patient.

“The model found the true drivers.” It measures how much this model uses a feature, not whether the feature causes anything. The scikit-learn guide says so directly: it “does not reflect the intrinsic predictive value of a feature by itself but how important this feature is for a particular model.”

Correlated features share the credit

Shuffle emergency admissions and the model can still recover some of that information from ED attendances, which mostly repeat it. The scikit-learn guide warns that this lowers the reported importance of both features, “though they might actually be important”. Shuffling a correlated pair together, with the same permutation for both columns, measures what the pair carries:

Together, admissions and ED attendances are worth 0.030, more than the two single numbers added up. The split within the pair is lopsided here because the model happened to lean on the column that drives the outcome. A different model, or a different seed, could split the credit more evenly, and the individual numbers would then look like two modest features instead of one important one.

In a stakeholder conversation, I report correlated features as a group: “prior hospital use”, not admissions and ED attendances separately. scikit-learn’s permutation_importance shuffles one column at a time, so grouped shuffling is a short loop of your own. The guide’s own suggestion is to cluster correlated features and keep one from each cluster.

Coefficients don’t escape the problem

Logistic regression looks easier to explain, because each feature has a coefficient. With correlated features, a coefficient is the effect of that feature holding all the others fixed, and “ED attendances rising while admissions stay the same” hardly happens in the data.

Data table

All nine features

Feature	Logistic regression	Logistic regression (low)	Logistic regression (high)
Age	0.19	0.11	0.27
Frailty score	0.35	0.25	0.42
Emergency admissions, 12 months	0.17	0.09	0.25
ED attendances, 12 months	-0.00	-0.07	0.09
Comorbidities	0.13	0.03	0.25
Medications	0.01	-0.10	0.10
Length of stay	0.14	0.09	0.18
Lives alone	0.11	0.05	0.15
Weekend discharge	0.02	-0.03	0.08

Without emergency admissions

Feature	Logistic regression	Logistic regression (low)	Logistic regression (high)
Age	0.19	0.11	0.27
Frailty score	0.35	0.25	0.42
ED attendances, 12 months	0.14	0.10	0.19
Comorbidities	0.13	0.03	0.25
Medications	0.01	-0.09	0.10
Length of stay	0.14	0.09	0.19
Lives alone	0.11	0.05	0.16
Weekend discharge	0.02	-0.03	0.08

Synthetic data. Source: projects/article-examples/round-two/explain_risk_model.py.

With all nine features, the coefficient for ED attendances is −0.003 per standard deviation, with a 90% bootstrap interval from −0.073 to 0.091: no effect. Remove emergency admissions and it becomes 0.143 (0.095 to 0.192). Nothing about ED attendances changed. Only the company it keeps did. The medication count sits near zero both times, with an interval as wide as the comorbidity coefficient’s, because the two are hard to tell apart.

So “the model says ED attendances don’t matter” and “the model says they do” can both be true statements about two reasonable models. Say which model, and group the features.

Partial dependence: the shape of an effect

Partial dependence shows how the average prediction changes as one feature moves, with every other feature left as it was. Individual conditional expectation (ICE) curves show the same thing for single patients. sklearn.inspection.partial_dependence returns both with kind="both":

res = partial_dependence(
    model,
    X,
    [feature],
    method="brute",
    kind="both",
    response_method="predict_proba",
    custom_values={feature: grid},
)
return res["grid_values"][0], res["average"][0], res["individual"][0][ice_rows]

method="brute" is needed for the individual curves, and response_method="predict_proba" keeps the output on the probability scale that people understand.

Data table

Age

Feature value	Average (partial dependence)	Patient 1	Patient 2	Patient 3	Patient 4	Patient 5
20	9%	13%	5%	7%	4%	4%
25	9%	13%	5%	7%	4%	4%
30	9%	13%	5%	7%	4%	4%
35	8%	11%	5%	6%	4%	4%
40	9%	12%	8%	6%	4%	6%
45	9%	13%	8%	7%	4%	6%
50	9%	12%	7%	7%	4%	6%
55	10%	14%	8%	8%	5%	7%
60	10%	13%	8%	8%	5%	6%
65	10%	13%	8%	7%	5%	6%
70	11%	14%	8%	8%	6%	8%
75	11%	14%	9%	8%	6%	8%
80	14%	20%	12%	13%	7%	10%
85	15%	22%	14%	13%	8%	11%
90	15%	22%	14%	13%	7%	10%
95	16%	22%	14%	14%	7%	10%
100	15%	22%	13%	14%	7%	10%

Length of stay

Feature value	Average (partial dependence)	Patient 1	Patient 2	Patient 3	Patient 4	Patient 5
0	13%	13%	17%	14%	5%	7%
1	13%	13%	17%	14%	5%	7%
2	13%	12%	14%	13%	4%	7%
3	12%	11%	14%	13%	4%	6%
4	12%	11%	13%	13%	4%	6%
5	13%	12%	14%	13%	4%	7%
6	13%	12%	14%	13%	4%	7%
7	13%	11%	14%	13%	4%	7%
8	13%	11%	14%	13%	4%	7%
9	13%	12%	14%	14%	5%	7%
10	13%	12%	14%	14%	5%	7%
11	14%	12%	15%	14%	5%	7%
12	14%	12%	15%	14%	5%	7%
13	14%	12%	15%	14%	5%	7%
14	13%	11%	15%	14%	4%	7%
15	15%	16%	15%	14%	6%	9%
16	15%	16%	15%	14%	6%	9%
17	17%	17%	16%	17%	8%	11%
18	19%	20%	19%	19%	9%	13%
19	22%	24%	25%	22%	10%	15%
20	22%	24%	25%	22%	10%	15%
21	22%	24%	25%	22%	10%	15%
22	26%	29%	30%	26%	14%	19%
23	26%	29%	30%	26%	14%	19%
24	26%	29%	30%	26%	14%	19%
25	26%	29%	30%	26%	14%	19%
26	26%	29%	30%	26%	14%	19%
27	26%	29%	30%	26%	14%	19%
28	26%	29%	30%	26%	14%	19%
29	26%	29%	30%	26%	14%	19%
30	26%	29%	30%	26%	14%	19%

Frailty score

Feature value	Average (partial dependence)	Patient 1	Patient 2	Patient 3	Patient 4	Patient 5
1	9%	7%	10%	13%	4%	6%
2	9%	7%	10%	13%	5%	6%
3	11%	10%	12%	18%	5%	6%
4	14%	13%	13%	23%	6%	7%
5	15%	13%	20%	30%	6%	8%
6	19%	15%	31%	40%	7%	9%
7	19%	15%	31%	40%	7%	9%
8	19%	15%	31%	40%	7%	9%
9	19%	15%	31%	40%	7%	9%

Synthetic data. Source: projects/article-examples/round-two/explain_risk_model.py.

For frailty, the average curve rises from 9.2% at a score of 1 to 19.1% at 6 and above. The individual lines tell a different story at a frailty score of 6: they run from 7.4% to 40.1%, and the two steepest belong to patients in their 80s who live alone. That’s the planted interaction, invisible in the average and in the importance chart.

Partial dependence has a catch, and it’s the correlated-features problem again. The scikit-learn guide states that PDPs and ICE curves “assume that the input features of interest are independent from the complement features”, and that with correlated features they “create absurd data points”. The age curve is an example. At age 40 it shows 8.9%. To get that number, the model scored every test patient as a 40-year-old while keeping their own frailty, including the frail ones. But none of the 328 test patients under 45 had a frailty score of 6 or more, and their mean score was 1.4 against 3.1 overall. Their observed readmission rate was 5.5%. The curve’s 8.9% describes a population that doesn’t exist.

Two more things to point out when you show these curves. Tree models go flat where the data runs thin: the length-of-stay curve stops rising at 22 days, while the simulated risk keeps climbing. And a curve is the model’s response to a change, not the effect of making that change.

Local: one patient, and what-ifs

Clinicians usually want the local question answered. A what-if changes one input and reports the new prediction. Patient A is 84, frailty score 7, lives alone, with two emergency admissions and three ED attendances in the last year, and a nine-day stay:

This is useful for showing what the model is paying attention to for this patient: frailty and living alone. It’s also easy to misread, in two ways.

First, the ED attendance scenario is impossible. In this data every emergency admission comes with an ED attendance, so a patient can’t have two admissions and no attendances. A what-if can create a patient who couldn’t exist, and the model will still give an answer.

Second, a what-if says what the model would predict, not what would happen. In the simulation, ED attendances have no effect on readmission at all, yet removing them lowers the prediction by five points, because the model uses them as a stand-in for admissions. Say “the model would score her lower”, never “her risk would fall”. In real data this matters most for things you could change. A long stay is partly a marker of how ill someone was, and the model can’t tell you what shortening it would do.

Uncertainty and limits

A single number like 47.3% sounds more precise than it is. One way to show how much it depends on the training sample is to refit the model on resampled data and look at the spread of predictions for the same patient:

Data table

Patient	Gradient boosting	Gradient boosting (low)	Gradient boosting (high)
Patient A	47%	35%	64%
Patient B	8%	7%	11%
Patient C	19%	16%	34%
Patient D	58%	47%	72%

Synthetic data. Source: projects/article-examples/round-two/explain_risk_model.py.

Across 30 refits, Patient A’s risk ranged from 35.4% to 63.6% (5th to 95th percentile). Patient B, a 58-year-old with no prior admissions, stayed between 6.7% and 11.4%. The spread is much wider for Patients A, C and D, whose combinations of age, frailty and history are rarer in the training data.

When I present a risk, I use these habits:

• Say it as a frequency. “About five in ten patients like her are readmitted within 30 days” is clearer than 47.3%, and it’s honest that the outcome is uncertain even when the risk is known.
• Give a range for unusual patients, and say the range comes from the model, not from the patient.
• Say what the model can’t see: social support beyond living alone, what happens after discharge, anything not in the record.
• Say what it wasn’t built for. A score built to rank patients for a phone call isn’t a basis for treatment decisions.

Model cards

A model card is a short document that travels with a model. Mitchell and colleagues proposed the format in Model Cards for Model Reporting (2019), with nine sections: model details, intended use, factors, metrics, evaluation data, training data, quantitative analyses, ethical considerations, and caveats and recommendations. To my mind, the intended-use section, especially the out-of-scope uses, does the most to prevent misuse. Here’s the skeleton I’d start from for this model:

model_details: Gradient boosted trees (scikit-learn HistGradientBoostingClassifier),
  nine features, trained on synthetic discharges. Owner, version and date go here.
intended_use:
  primary: Rank adult discharges for a follow-up phone call.
  users: Discharge team, who review the list before calling.
  out_of_scope: Treatment decisions; patients under 18; any use without the list review.
factors: Age band, sex, deprivation quintile, frailty.
metrics: AUROC, calibration (mean predicted against observed, and by risk decile),
  readmissions reached by the top 10% of the list.
evaluation_data: 3,600 held-out discharges; test AUROC 0.696, mean predicted 13.1%
  against 13.0% observed.
training_data: 8,400 synthetic discharges.
quantitative_analyses: The metrics above for each factor, with intervals.
ethical_considerations: Ranking by risk can send most calls to the oldest patients;
  say who is left out.
caveats_and_recommendations: Correlated features, so report importance by group.
  Check calibration monthly. Retrain if it drifts.

The 30-day readmission risk model shows what the quantitative analyses look like, including who ends up on the call list.

Before the meeting

1. Decide whether the question is global or local, and pick the tool to match.
2. Group correlated features before reporting importance.
3. Show ICE curves next to any partial dependence plot, and point out where the grid runs into combinations that don’t exist.
4. Phrase what-ifs as statements about the model.
5. Give a range and a frequency, not a bare percentage.
6. Bring the model card, and read out the out-of-scope uses.

Reproduce

The data, both models, every table and the four charts come from one script, which runs in under a minute:

uv run python projects/article-examples/round-two/explain_risk_model.py