Model Interpretability

Real-World Scenario

A bank’s loan model is challenged by a regulator who wants to know why a customer was denied. The data scientist can’t say “the black box said no.” With SHAP values, they produce a per-applicant breakdown: “Credit history contributed −0.42 to the score, income contributed +0.28, and debt ratio contributed −0.31.” The regulator is satisfied.

Feature Importances (Tree Models)

from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier, HistGradientBoostingClassifier
from sklearn.inspection import permutation_importance
from sklearn.model_selection import train_test_split
import numpy as np
import pandas as pd

FEATURE_NAMES = [
    "income", "age", "credit_score", "debt_ratio", "employment_years",
    "num_accounts", "late_payments", "loan_amount", "assets", "expenses",
]
X, y = make_classification(n_samples=3000, n_features=10, n_informative=7,
                            random_state=42)
X_df = pd.DataFrame(X, columns=FEATURE_NAMES)
X_train, X_test, y_train, y_test = train_test_split(X_df, y, test_size=0.2, random_state=42)

model = RandomForestClassifier(n_estimators=200, random_state=42, n_jobs=-1)
model.fit(X_train, y_train)

# Built-in impurity-based importance (fast but biased toward high-cardinality features)
impurity_imp = pd.Series(model.feature_importances_, index=FEATURE_NAMES)
print("Impurity-based importance (biased — use as starting point only):")
print(impurity_imp.sort_values(ascending=False).round(4).to_string())

# Permutation importance (unbiased — measures actual predictive contribution)
result = permutation_importance(
    model, X_test, y_test,
    n_repeats=20, random_state=42, scoring="roc_auc", n_jobs=-1
)
perm_imp = pd.DataFrame({
    "mean":  result.importances_mean,
    "std":   result.importances_std,
}, index=FEATURE_NAMES).sort_values("mean", ascending=False)

print("\nPermutation importance (test set AUC drop):")
for feat, row in perm_imp.iterrows():
    print(f"  {feat:<20} {row['mean']:>8.4f} ± {row['std']:.4f}")

$ python mi_importance.py
Impurity-based importance (biased — use as starting point only):
debt_ratio          0.2858
num_accounts        0.1499
assets              0.0982
income              0.0868
age                 0.0826
credit_score        0.0765
late_payments       0.0713
employment_years    0.0677
expenses            0.0550
loan_amount         0.0260

Permutation importance (test set AUC drop):
  debt_ratio             0.2070 ± 0.0132
  num_accounts           0.0337 ± 0.0055
  income                 0.0235 ± 0.0036
  assets                 0.0226 ± 0.0036
  late_payments          0.0174 ± 0.0032
  credit_score           0.0165 ± 0.0033
  age                    0.0123 ± 0.0020
  employment_years       0.0035 ± 0.0017
  expenses               0.0020 ± 0.0011
  loan_amount           -0.0002 ± 0.0005

Both methods agree debt_ratio dominates, but they disagree on the runner-up: impurity-based importance ranks num_accounts second (0.1499) while permutation importance ranks it second too but at a much smaller relative gap to income/assets — and loan_amount’s permutation importance is slightly negative, meaning shuffling it occasionally improves AUC by chance. That’s the tell that loan_amount contributes essentially nothing real to this model, something the impurity score (0.0260, ranked dead last but still positive) doesn’t make as obvious.

SHAP Values

# pip install shap
import shap
import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import train_test_split

FEATURE_NAMES = [
    "income", "age", "credit_score", "debt_ratio", "employment_years",
    "num_accounts", "late_payments", "loan_amount", "assets", "expenses",
]

from sklearn.datasets import make_classification
X, y = make_classification(n_samples=3000, n_features=10, n_informative=7, random_state=42)
X_df = pd.DataFrame(X, columns=FEATURE_NAMES)
X_train, X_test, y_train, y_test = train_test_split(X_df, y, test_size=0.2, random_state=42)

model = HistGradientBoostingClassifier(max_iter=200, random_state=42)
model.fit(X_train, y_train)

# TreeExplainer: optimized for gradient boosting and random forests
explainer   = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)   # shape: (n_samples, n_features)

# SHAP values for a single prediction
sample_idx = 0
sample     = X_test.iloc[sample_idx]
sv         = shap_values[sample_idx]

print(f"Sample prediction: {model.predict_proba(sample.values.reshape(1, -1))[0, 1]:.4f}")
base_value = np.asarray(explainer.expected_value).item()
print(f"Base value (global mean): {base_value:.4f}")
print(f"\nFeature contributions:")
contributions = pd.Series(sv, index=FEATURE_NAMES).sort_values(key=abs, ascending=False)
for feat, val in contributions.items():
    direction = "↑" if val > 0 else "↓"
    print(f"  {feat:<20} {direction} {val:+.4f}")

# Global feature importance (mean |SHAP|)
mean_abs_shap = np.abs(shap_values).mean(axis=0)
global_imp    = pd.Series(mean_abs_shap, index=FEATURE_NAMES).sort_values(ascending=False)
print(f"\nGlobal SHAP importance (mean |SHAP|):")
print(global_imp.round(4).to_string())

$ python mi_shap_tree.py
Sample prediction: 0.9999
Base value (global mean): -0.2553

Feature contributions:
  num_accounts         ↑ +3.3514
  debt_ratio           ↑ +1.5280
  credit_score         ↑ +1.3816
  assets               ↑ +1.2705
  income               ↑ +0.9017
  late_payments        ↑ +0.7781
  employment_years     ↑ +0.2291
  age                  ↑ +0.2233
  expenses             ↓ -0.2139
  loan_amount          ↓ -0.1301

Global SHAP importance (mean |SHAP|):
debt_ratio          4.6288
late_payments       1.7771
num_accounts        1.6644
assets              1.1772
credit_score        1.1200
income              0.9961
age                 0.8900
expenses            0.6131
employment_years    0.3105
loan_amount         0.2445

In this shap release, TreeExplainer.expected_value comes back as a 1-element array rather than a scalar, so it needs np.asarray(...).item() before formatting — without that fix f"{explainer.expected_value:.4f}" raises a TypeError. The base value (-0.2553) is in log-odds space, not a probability, which is why it’s negative even though this particular sample’s predicted probability (0.9999) is near-certain: the ten feature contributions summed with the base value push the log-odds strongly positive before the sigmoid converts it back to a probability near 1. Note that this sample’s single largest driver by raw contribution (num_accounts) isn’t the same feature as the model’s largest global driver (debt_ratio) — local and global importance answer different questions.

SHAP for Any Model (KernelExplainer)

import shap
import numpy as np
import pandas as pd
from sklearn.neural_network import MLPClassifier
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split

FEATURE_NAMES = ["f1","f2","f3","f4","f5","f6","f7","f8","f9","f10"]
X, y = make_classification(n_samples=1000, n_features=10, random_state=42)
X_df = pd.DataFrame(X, columns=FEATURE_NAMES)
X_train, X_test, y_train, y_test = train_test_split(X_df, y, test_size=0.2, random_state=42)

# Neural network (tree SHAP doesn't apply)
pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("mlp",    MLPClassifier(hidden_layer_sizes=(64, 32), max_iter=200, random_state=42)),
])
pipe.fit(X_train, y_train)

# KernelExplainer works on any model but is slower
# Use a small background set (50-100 samples is usually enough)
background  = shap.sample(X_train, 50, random_state=42)
explainer   = shap.KernelExplainer(
    lambda x: pipe.predict_proba(pd.DataFrame(x, columns=FEATURE_NAMES))[:, 1],
    background,
)

# Explain a small subset (KernelSHAP is slow)
X_explain   = X_test.iloc[:20]
shap_values = explainer.shap_values(X_explain, nsamples=100)

mean_abs = np.abs(shap_values).mean(axis=0)
imp = pd.Series(mean_abs, index=FEATURE_NAMES).sort_values(ascending=False)
print("Neural network SHAP importance:")
print(imp.round(4).to_string())

$ python mi_shap_kernel.py
Neural network SHAP importance:
f7     0.1769
f3     0.1313
f9     0.0856
f1     0.0730
f4     0.0268
f2     0.0244
f6     0.0233
f5     0.0206
f10    0.0169
f8     0.0166

KernelExplainer treats the model as a black box — it perturbs each of the 20 explained samples against the 50-sample background and fits a local linear model, which is why it works on an MLPClassifier where TreeExplainer doesn’t apply. The trade-off shows up as a ConvergenceWarning from the MLP itself (200 iterations wasn’t quite enough for this random split to fully converge) plus a visibly slower run than TreeExplainer even at just 20 samples and 100 perturbations — KernelExplainer doesn’t scale the way tree-specific SHAP does.

Partial Dependence Plots (PDP)

from sklearn.inspection import partial_dependence, PartialDependenceDisplay
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
import numpy as np
import pandas as pd

FEATURE_NAMES = [
    "income", "age", "credit_score", "debt_ratio", "employment_years",
    "num_accounts", "late_payments", "loan_amount", "assets", "expenses",
]
X, y = make_classification(n_samples=3000, n_features=10, n_informative=7, random_state=42)
X_df = pd.DataFrame(X, columns=FEATURE_NAMES)
X_train, X_test, y_train, y_test = train_test_split(X_df, y, test_size=0.2, random_state=42)

model = HistGradientBoostingClassifier(max_iter=200, random_state=42)
model.fit(X_train, y_train)

# Compute partial dependence for the top feature
pd_result = partial_dependence(
    model, X_train,
    features=[0],                     # feature index (income)
    kind="average",                   # "average" = PDP, "individual" = ICE curves
    method="brute",                   # required for response_method="predict_proba"
    response_method="predict_proba",  # probabilities, not the default log-odds
    percentiles=(0.05, 0.95),
    grid_resolution=50,
)
values   = pd_result["grid_values"][0]
avg_pred = pd_result["average"][0]

print("Partial dependence of 'income' on prediction:")
print(f"{'Income (std)':<15} {'Predicted prob':>15}")
print("-" * 32)
for v, p in zip(values[::5], avg_pred[::5]):
    bar = "█" * int(p * 30)
    print(f"{v:15.2f} {p:15.4f}  {bar}")

# 2D interaction PDP
pd_2d = partial_dependence(
    model, X_train,
    features=[(0, 2)],     # income × credit_score interaction
    kind="average",
    grid_resolution=20,
)
print(f"\n2D PDP shape: {pd_2d['average'][0].shape}")  # (20, 20) grid

$ python mi_pdp.py
Partial dependence of 'income' on prediction:
Income (std)     Predicted prob
--------------------------------
          -3.11          0.3864  ███████████
          -2.52          0.3929  ███████████
          -1.93          0.4138  ████████████
          -1.33          0.4641  █████████████
          -0.74          0.5131  ███████████████
          -0.15          0.5316  ███████████████
           0.44          0.5498  ████████████████
           1.03          0.5524  ████████████████
           1.63          0.5332  ███████████████
           2.22          0.5439  ████████████████

2D PDP shape: (20, 20)

partial_dependence’s default method="auto" picks the "recursion" method for gradient-boosted trees, and "recursion" only supports response_method="decision_function" — it silently returns log-odds, not probabilities, even though nothing in the call suggested that. Left at the defaults, the “Predicted prob” column would print negative numbers and the int(p * 30) bar-width calculation would silently produce empty bars for every negative value instead of erroring. Passing method="brute" with response_method="predict_proba" explicitly, as done above, is what actually gets a 0-to-1 probability curve — which here rises from 0.39 to a plateau around 0.55 as income increases, flattening out rather than climbing indefinitely.

LIME for Individual Predictions

# pip install lime
from lime.lime_tabular import LimeTabularExplainer
import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split

FEATURE_NAMES = [
    "income", "age", "credit_score", "debt_ratio", "employment_years",
    "num_accounts", "late_payments", "loan_amount", "assets", "expenses",
]
X, y = make_classification(n_samples=3000, n_features=10, n_informative=7, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

model = HistGradientBoostingClassifier(max_iter=200, random_state=42)
model.fit(X_train, y_train)

explainer = LimeTabularExplainer(
    training_data=X_train,
    feature_names=FEATURE_NAMES,
    class_names=["approved", "denied"],
    mode="classification",
    random_state=42,
)

# Explain a single prediction
sample_idx  = 0
sample      = X_test[sample_idx]
explanation = explainer.explain_instance(
    sample,
    model.predict_proba,
    num_features=5,    # top 5 contributing features
    num_samples=1000,
)

pred_prob = model.predict_proba(sample.reshape(1, -1))[0, 1]
print(f"Prediction: denied with probability {pred_prob:.3f}")
print(f"\nTop 5 contributing features (LIME):")
for feat, weight in explanation.as_list():
    direction = "increases denial" if weight > 0 else "decreases denial"
    print(f"  {feat:<35}  {weight:+.4f}  ({direction})")

$ python mi_lime.py
Prediction: denied with probability 1.000

Top 5 contributing features (LIME):
  num_accounts <= -0.46                +0.3051  (increases denial)
  0.01 < debt_ratio <= 1.28            +0.1132  (increases denial)
  -0.86 < credit_score <= 0.42         +0.0486  (increases denial)
  assets <= -0.70                      +0.0479  (increases denial)
  0.07 < age <= 1.26                   +0.0476  (increases denial)

This is the same test sample SHAP explained earlier, and the two methods agree on direction: num_accounts was SHAP’s single largest per-sample driver (+3.3514, pushing toward the positive class) and LIME independently ranks it first here too. LIME reports its explanation as feature ranges (num_accounts <= -0.46) rather than SHAP’s exact per-feature values, because it fits a local linear model over binned/perturbed neighbors around this one prediction rather than computing an exact game-theoretic attribution — faster per call, but the explanation is only valid in the neighborhood LIME sampled around this specific point.

Practice

1. In the feature importance block, change permutation_importance’s n_repeats from 20 to 5. Do the std values grow, and why would that make loan_amount’s near-zero importance harder to trust?

   Solution

   Fewer repeats means each feature’s importance estimate averages over less randomness, so std grows — for a feature like loan_amount whose true importance is already near zero, a noisier estimate makes it harder to tell “genuinely unimportant” apart from “unlucky shuffle.”

2. In the SHAP TreeExplainer block, change the model from HistGradientBoostingClassifier to RandomForestClassifier(n_estimators=200, random_state=42) and compare the global SHAP importance ranking to the lesson’s.

   Solution

   The top feature (debt_ratio) should stay consistent since it’s genuinely the strongest signal in this synthetic dataset regardless of which tree ensemble learns it, but the exact ranking of the weaker features can shift because Random Forest and HistGradientBoosting build trees differently.

3. In the KernelExplainer block, reduce background from 50 to 10 samples and increase nsamples from 100 to 500. How does the importance ranking compare, and which change would you expect to matter more for accuracy — background size or nsamples?

   Solution

   nsamples (how many perturbations KernelExplainer evaluates per explained point) usually matters more for stability than background size, since a too-small background set mainly widens the “reference” distribution against which contributions are measured, while too few nsamples makes the local linear fit itself noisy.

4. In the PDP block, compute partial dependence for debt_ratio (feature index 3) instead of income, using the same method="brute", response_method="predict_proba" fix. Does the curve rise or fall as debt_ratio increases?

   Solution

   Expect the curve to fall as debt_ratio increases — higher debt-to-income ratios should push the predicted probability down, the opposite direction from income’s rising curve, matching debt_ratio’s large negative-leaning SHAP contributions on the sample explained earlier in the lesson (once you account for label direction).

5. In the LIME block, increase num_samples from 1000 to 100 (fewer perturbations). Run the explanation twice with the same random_state — does it still produce identical output, and what does that imply about num_samples and stability?

   Solution

   With the same random_state, output stays reproducible run-to-run, but with fewer perturbations (100 instead of 1000) the local linear fit has less data to work with, so a different random_state would be more likely to produce a noticeably different feature ranking than it would at 1000 samples — LIME’s explanations are only as stable as the sampling budget you give them.