SHAP — explaining model predictions

You’ve shipped a credit-scoring model. A customer is denied. They call, and your support team has to tell them why. Or a regulator audits the model and asks how it decided. “The feature importances are highest for debt_ratio” is not a usable answer — that’s about the model overall, not about this customer.

You need a per-prediction explanation. That’s what SHAP gives you — it explains the model’s behavior, not necessarily the true cause of the outcome in the real world.

Why .feature_importances_ is misleading

Tree ensembles expose .feature_importances_. It measures how often each feature was used at a split, weighted by impurity reduction. That tells you about the structure of the trees, not about what the model predicted on a specific input.

Two well-known problems:

• Biased toward high-cardinality features. Features with many unique values (continuous numerics, IDs) win more splits simply because they’re more splittable. This inflates their importance even if they carry no real signal.
• No notion of per-row contribution. It’s a single number per feature for the whole model. Useless when a stakeholder asks “why did the model deny this applicant?”

For a single global ranking, permutation importance (shuffle a feature on held-out data, measure the score drop) is already a big upgrade. But for per-prediction explanations, SHAP is the standard.

What SHAP values actually are

SHAP — SHapley Additive exPlanations — comes from cooperative game theory. Treat each feature as a “player” contributing to the prediction. The Shapley value for feature i answers: across all possible subsets of features (coalitions) that could be presented to the model, what is the average marginal contribution of adding feature i to that subset?

The result is exactly the budget you want. For a single prediction:

prediction = base_value + Σ shap_value_i

where base_value is the model’s average prediction across the training set, and each shap_value_i is feature i’s push (positive or negative) on this specific row. Sum them up and you reconstruct the model’s prediction exactly. Per-feature. Per-row. Locally exact.

For tree models, the TreeExplainer computes SHAP values in polynomial time using the tree structure — fast enough for production use on millions of rows.

Explain a credit-decision model

import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(0)
n = 1500
df = pd.DataFrame({
    "credit_score":  rng.integers(300, 850, n),
    "income":        rng.gamma(2, 30000, n).round(0),
    "debt_ratio":    rng.beta(2, 5, n).round(3),
    "delinquencies": rng.poisson(0.4, n),
    "open_accounts": rng.integers(1, 15, n),
})
risk = (
    -0.006 * df["credit_score"]
    + 1.8  * df["debt_ratio"]
    + 0.5  * df["delinquencies"]
    + rng.normal(0, 0.4, n)
)
df["defaulted"] = (risk > 0.7).astype(int)

X, y = df.drop(columns=["defaulted"]), df["defaulted"]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, stratify=y, random_state=0)

model = HistGradientBoostingClassifier(
    max_iter=300, learning_rate=0.05, max_depth=5, random_state=0
).fit(X_tr, y_tr)

print(f"test accuracy: {model.score(X_te, y_te):.3f}")

# A denied applicant — pick one with high predicted default probability
proba = model.predict_proba(X_te)[:, 1]
i = int(np.argmax(proba))
applicant = X_te.iloc[i]
print(f"\napplicant predicted default probability: {proba[i]:.3f}")
print("\napplicant features:")
print(applicant)

Approximating SHAP — feature ablation

A simple way to decompose a prediction into per-feature contributions: hold each feature out and see how much the prediction moves. This isn’t exactly Shapley (it ignores feature ordering effects) but the intuition and the shape of the output match.

import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(0)
n = 1500
df = pd.DataFrame({
    "credit_score":  rng.integers(300, 850, n),
    "income":        rng.gamma(2, 30000, n).round(0),
    "debt_ratio":    rng.beta(2, 5, n).round(3),
    "delinquencies": rng.poisson(0.4, n),
    "open_accounts": rng.integers(1, 15, n),
})
risk = (
    -0.006 * df["credit_score"]
    + 1.8  * df["debt_ratio"]
    + 0.5  * df["delinquencies"]
    + rng.normal(0, 0.4, n)
)
df["defaulted"] = (risk > 0.7).astype(int)

X, y = df.drop(columns=["defaulted"]), df["defaulted"]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, stratify=y, random_state=0)

model = HistGradientBoostingClassifier(
    max_iter=300, learning_rate=0.05, max_depth=5, random_state=0
).fit(X_tr, y_tr)

# Pick the most-likely-default applicant
proba_test = model.predict_proba(X_te)[:, 1]
i = int(np.argmax(proba_test))
row = X_te.iloc[[i]]

# Baseline: prediction averaged over training set (the "base value")
base_value = model.predict_proba(X_tr)[:, 1].mean()
actual = model.predict_proba(row)[0, 1]

# Approximate per-feature contribution by replacing each feature with the train mean
contribs = {}
mean_row = X_tr.mean()
for col in X.columns:
    counterfactual = row.copy()
    counterfactual[col] = mean_row[col]
    cf_pred = model.predict_proba(counterfactual)[0, 1]
    contribs[col] = actual - cf_pred  # how much THIS feature pushed the prediction

print(f"base value         : {base_value:.3f}")
print(f"actual prediction  : {actual:.3f}")
print(f"\nper-feature contribution to the prediction:")
for f, v in sorted(contribs.items(), key=lambda x: -abs(x[1])):
    arrow = "↑ risk" if v > 0 else "↓ risk"
    print(f"  {f:<14} {v:+.3f}  {arrow}  (value: {row.iloc[0][f]})")

You now have an answer you can give the customer:

“Your application scored a default probability of 0.71 compared to the average of 0.30. The biggest drivers were your credit_score=412 contributing +0.18 to risk, delinquencies=3 contributing +0.14, and debt_ratio=0.62 contributing +0.09.”

This is exactly what a regulator wants to see, and it’s exactly what SHAP’s force_plot visualizes — a horizontal arrow chart from the base value to the final prediction, with each feature pushing it left or right.

The three SHAP plots you’ll actually use

When you use the real shap library, three visualizations cover 95% of the work:

The summary plot is what replaces feature_importances_. It tells you both which features matter and how high/low values push the prediction, which the impurity-based importance can never tell you.

Limits of SHAP

Two things to keep in mind:

• SHAP explains the model, not the world. If your model has learned a spurious correlation, SHAP will faithfully report that the spurious feature drove the prediction. SHAP doesn’t validate causality — it describes the model’s behavior.
• TreeExplainer is fast; KernelExplainer is not. For non-tree models (neural networks, SVMs), SHAP falls back to sampling-based estimation which can be slow on big datasets. Tree models are the sweet spot.

In one breath

• .feature_importances_ is misleading — it measures split structure (biased toward high-cardinality features) and gives one global number per feature, never a per-row explanation.
• SHAP = Shapley values from game theory: each feature’s average marginal contribution to a prediction, so prediction = base_value + Σ shap_i — locally exact, per-row, auditable.
• TreeExplainer computes them fast for tree models; the summary plot replaces feature_importances_, the force plot explains one decision, the dependence plot shows interactions.
• It’s often legally required in regulated domains (GDPR right to explanation, ECOA adverse-action reasons) — the standard because the contributions sum exactly and are consistent.
• SHAP explains the model, not the world — it faithfully reports a spurious/proxy feature as a driver, so it is not a causality or fairness check; and KernelExplainer (non-tree) is slow.

Quick check