Feature Importance and SHAP

The Production Scenario

A consumer lending bank has deployed a gradient-boosted tree model to automate loan decisions. The model performs well - AUC 0.91, default rate down 18% year-over-year. Then a regulator from the Consumer Financial Protection Bureau walks in with a list of 47 specific applicants who were denied and demands an explanation for each one. Not a general model summary. Each person, individually, with a reason that would hold up in a legal proceeding.

The ML engineer pulls up the XGBoost dashboard and screenshots the feature importance bar chart. "Income" is at the top with importance 0.34. "Debt-to-income ratio" is second at 0.28. The regulator shakes her head. "That tells me what your model generally pays attention to. It does not tell me why Maria Gonzalez, who earns $72,000 annually and has no late payments, was denied a$15,000 personal loan. You are in violation of ECOA adverse action notice requirements."

The engineer goes quiet. The global feature importance chart - the one they put in every model card and stakeholder presentation - is completely useless for this question. Worse, the model uses 23 features, and without knowing which ones pushed Maria's prediction toward "deny," the bank cannot issue a legally compliant adverse action letter. The bank's legal exposure is now in the millions.

Three days and a SHAP integration later, the engineer can generate a precise statement: "This application was declined primarily because the requested loan amount exceeds 2.1x your annual income (increases risk score by +0.31), your credit utilization is 78% (increases risk score by +0.18), and your employment tenure of 4 months is below our 12-month threshold (+0.09). Your clean payment history reduced risk by -0.14." That is an adverse action notice. That is SHAP.

Why Feature Importance Is Hard

Tree models produce multiple competing answers to "which features matter?" Each answer is statistically valid, measures something different, and can disagree substantially with the others. Understanding what each measures - and when each misleads - is the foundation of responsible ML deployment.

The three fundamental questions are different:

• Which features does the model structurally rely on? (MDI answers this)
• Which features are actually predictive on held-out data? (Permutation answers this)
• Which features drove this specific prediction? (SHAP answers this)

Conflating these questions produces wrong answers and, in regulated industries, legal exposure.

MDI: Mean Decrease in Impurity

MDI measures how much each feature contributes to impurity reduction across all splits in all trees:

$\text{MDI}(f_j) = \frac{1}{|\text{Trees}|}\sum_{t \in \text{Trees}}\sum_{\substack{n \in t \\ n\text{ splits on }f_j}} \frac{n_t}{n} \cdot \Delta I(n)$

where $n_t$ is the number of training samples reaching node $n$, $n$ is the total training size, and $\Delta I(n)$ is the impurity reduction at node $n$. The weighting $n_t/n$ gives more credit to splits near the root, which affect more examples.

MDI Bias with High-Cardinality Features: A Demonstration

The structural bias of MDI: continuous features with many unique values get more split opportunities, systematically inflating their importance scores.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.inspection import permutation_importance
from sklearn.model_selection import train_test_split
import xgboost as xgb

# ── Synthetic loan dataset ────────────────────────────────────────────────────
rng = np.random.default_rng(42)
n = 6000

# High-cardinality continuous features (many unique values)
income           = rng.lognormal(10.8, 0.5, n)
loan_amount      = rng.lognormal(9.6, 0.6, n)
credit_score     = np.clip(rng.normal(680, 80, n), 300, 850)
dti_ratio        = rng.beta(2, 5, n) * 0.8
employment_months = rng.exponential(36, n).clip(0, 240)
num_late_payments = rng.poisson(0.8, n)
credit_util      = rng.beta(3, 4, n)

# Low-cardinality categorical feature (only 4 distinct values - but actually useful!)
# Encoded as 0/1/2/3 = lease/used/own/contract
housing_type = rng.integers(0, 4, n).astype(float)

# Near-constant feature: should have near-zero importance
account_flag = rng.integers(0, 2, n).astype(float) * rng.normal(1, 0.01, n)

df = pd.DataFrame({
    "income":             income,
    "loan_amount":        loan_amount,
    "credit_score":       credit_score,
    "dti_ratio":          dti_ratio,
    "employment_months":  employment_months,
    "num_late_payments":  num_late_payments,
    "credit_util":        credit_util,
    "housing_type":       housing_type,       # low cardinality, real signal
    "account_flag":       account_flag,       # near-constant, minimal signal
})

# Target: default probability (housing_type is actually moderately predictive)
log_odds = (
    -3.0
    - 0.4 * (income / 50000)
    + 0.3 * (loan_amount / income)
    + 0.5 * (dti_ratio - 0.3)
    - 0.3 * ((credit_score - 680) / 100)
    + 0.4 * num_late_payments
    + 0.3 * credit_util
    - 0.2 * (employment_months / 12)
    + 0.15 * (housing_type == 1)  # renters slightly higher risk
)
df["default"] = (rng.logistic(log_odds) > 0).astype(int)
print(f"Default rate: {df['default'].mean():.1%}")

X = df.drop("default", axis=1)
y = df["default"]
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# ── Train XGBoost ─────────────────────────────────────────────────────────────
model = xgb.XGBClassifier(
    n_estimators=200, max_depth=5, learning_rate=0.05,
    subsample=0.8, colsample_bytree=0.8,
    random_state=42, eval_metric="auc", verbosity=0,
)
model.fit(
    X_train, y_train,
    eval_set=[(X_test, y_test)],
    verbose=False,
)
from sklearn.metrics import roc_auc_score
print(f"Test AUC: {roc_auc_score(y_test, model.predict_proba(X_test)[:,1]):.4f}")

# ── Compare all importance types ──────────────────────────────────────────────
print("\n=== MDI (split count, gain, coverage) ===")
for imp_type, label in [("weight", "Split Count"), ("gain", "Gain"), ("cover", "Coverage")]:
    scores = model.get_booster().get_score(importance_type=imp_type)
    top3 = sorted(scores.items(), key=lambda x: x[1], reverse=True)[:3]
    print(f"{label}: {[(f, f'{v:.3f}') for f, v in top3]}")

print("\n=== Permutation Importance (unbiased) ===")
perm = permutation_importance(
    model, X_test, y_test,
    n_repeats=15, scoring="roc_auc", random_state=42, n_jobs=-1
)
perm_df = pd.DataFrame({
    "feature": X.columns,
    "importance": perm.importances_mean,
    "std": perm.importances_std,
}).sort_values("importance", ascending=False)
print(perm_df.to_string(index=False))

MDA: Mean Decrease Accuracy (Permutation Importance)

Permutation importance breaks the relationship between a feature and the target by shuffling that feature's values in the validation set:

$\text{MDA}(f_j) = \text{Score}_{\text{baseline}} - \mathbb{E}_{\text{shuffles}}\left[\text{Score}_{\text{permuted}(f_j)}\right]$

Key properties:

• Model-agnostic: works with any black-box model
• Captures interaction effects: if removing feature $A$ breaks interactions with $B$, that registers
• Measured on actual predictive performance: not internal tree structure

The correlated feature problem: If income and income_monthly are correlated at 0.95, permuting income still leaves income_monthly to carry the signal. Both features appear unimportant even though together they are critical. The fix: grouped permutation importance - permute all correlated features simultaneously.

def grouped_permutation_importance(
    model,
    X: pd.DataFrame,
    y: np.ndarray,
    feature_groups: dict,
    n_repeats: int = 10,
    scoring: str = "roc_auc",
) -> pd.DataFrame:
    """
    Permute entire groups of correlated features simultaneously.
    feature_groups: {'group_name': ['feat1', 'feat2'], ...}
    """
    from sklearn.metrics import roc_auc_score
    baseline = roc_auc_score(y, model.predict_proba(X)[:, 1])

    results = []
    for group_name, features in feature_groups.items():
        drops = []
        for _ in range(n_repeats):
            X_perm = X.copy()
            # Permute all features in the group simultaneously
            perm_idx = np.random.permutation(len(X))
            for feat in features:
                if feat in X_perm.columns:
                    X_perm[feat] = X_perm[feat].values[perm_idx]
            score = roc_auc_score(y, model.predict_proba(X_perm)[:, 1])
            drops.append(baseline - score)

        results.append({
            "group": group_name,
            "features": features,
            "importance_mean": np.mean(drops),
            "importance_std": np.std(drops),
        })

    return pd.DataFrame(results).sort_values("importance_mean", ascending=False)


# Example: group highly correlated features
feature_groups = {
    "income_signals": ["income", "loan_amount", "dti_ratio"],
    "credit_history":  ["credit_score", "num_late_payments", "credit_util"],
    "employment":      ["employment_months"],
    "housing":         ["housing_type"],
    "noise":           ["account_flag"],
}

grouped_imp = grouped_permutation_importance(model, X_test, y_test, feature_groups)
print("\n=== Grouped Permutation Importance ===")
print(grouped_imp.to_string(index=False))

SHAP: SHapley Additive Explanations

SHAP is grounded in cooperative game theory. The central question: if a group of players (features) cooperate to produce an outcome (prediction), how do we fairly distribute credit among them?

The Shapley value for feature $i$, proven unique by Lloyd Shapley (1953), is:

$\phi_i = \sum_{S \subseteq F \setminus \{i\}} \frac{|S|!\,(|F|-|S|-1)!}{|F|!}\left[v(S \cup \{i\}) - v(S)\right]$

In plain terms: average the marginal contribution of feature $i$ across all possible orderings in which features could be introduced to the model. For each ordering, measure how much including feature $i$ changes the prediction compared to not including it.

The Three SHAP Axioms (Uniqueness Proof)

SHAP is the unique attribution method satisfying all three simultaneously:

Axiom	Mathematical form	Meaning
Efficiency	$\sum_{i=1}^{d} \phi_i = f(x) - \mathbb{E}[f(X)]$	Attributions sum exactly to prediction minus baseline
Symmetry	$v(S \cup \{i\}) = v(S \cup \{j\})$ for all $S \Rightarrow \phi_i = \phi_j$	Equal contributors get equal attribution
Dummy	$v(S \cup \{i\}) = v(S)$ for all $S \Rightarrow \phi_i = 0$	Useless features get zero attribution

The efficiency axiom is why SHAP is suitable for adverse action notices and audits: the attributions add up exactly to the model's output minus the baseline. You can show a customer precisely how their prediction was assembled from individual feature contributions, and the math is guaranteed to be consistent.

TreeSHAP: Exact and Polynomial-Time

Naively computing Shapley values requires evaluating the model on all $2^{|F|}$ feature subsets - exponential in the number of features. For a model with 30 features that is over a billion evaluations.

TreeSHAP (Lundberg et al., 2018) exploits the recursive structure of decision trees to compute exact Shapley values in $O(T L D^2)$ time, where $T$ is the number of trees, $L$ is the number of leaves, and $D$ is the maximum tree depth. For a typical XGBoost model with 200 trees, depth 6, and 1000 leaves per tree, SHAP values compute in milliseconds per sample.

The key insight: a decision tree partitions the feature space into rectangular regions. For any leaf, we know exactly which features determined the path. TreeSHAP uses dynamic programming over the tree structure to accumulate Shapley values without enumerating all subsets.

Complete SHAP Implementation

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

# ── Create TreeExplainer ──────────────────────────────────────────────────────
# Use a real background dataset for honest "absent feature" marginalization
background = shap.sample(X_train, 500, random_state=42)
explainer = shap.TreeExplainer(model, background)

# ── Compute SHAP values for test set ─────────────────────────────────────────
shap_values = explainer(X_test)  # Returns Explanation object (n_test x d)

# ── Verify the efficiency axiom ───────────────────────────────────────────────
sample_idx = 0
shap_sum      = shap_values[sample_idx].values.sum()
model_logit   = model.predict(xgb.DMatrix(X_test.iloc[[sample_idx]]))[0]
baseline      = explainer.expected_value
print(f"SHAP sum:              {shap_sum:.6f}")
print(f"prediction - baseline: {model_logit - baseline:.6f}")
print(f"Efficiency check passed: {abs(shap_sum - (model_logit - baseline)) < 1e-5}")

# ── Global: mean absolute SHAP importance ─────────────────────────────────────
mean_abs_shap = np.abs(shap_values.values).mean(axis=0)
shap_importance = pd.DataFrame({
    "feature": X.columns,
    "mean_abs_shap": mean_abs_shap,
}).sort_values("mean_abs_shap", ascending=False)
print("\n=== SHAP Global Importance (mean|φ|) ===")
print(shap_importance.to_string(index=False))

SHAP Visualization Plots

# ── 1. Summary plot: replaces global bar chart ────────────────────────────────
# Each dot = one prediction; x-axis = SHAP value; color = feature value
plt.figure(figsize=(10, 6))
shap.summary_plot(shap_values, X_test, plot_type="dot", show=False)
plt.tight_layout()
plt.savefig("shap_summary.png", dpi=150)
plt.close()

# ── 2. Waterfall: single prediction explanation (Maria's application) ─────────
# Find a denied applicant (predicted high risk)
pred_proba = model.predict_proba(X_test)[:, 1]
denied_idx = np.where(pred_proba > 0.7)[0][0]  # first high-risk prediction

shap.plots.waterfall(shap_values[denied_idx], show=False, max_display=10)
plt.savefig("shap_waterfall.png", dpi=150, bbox_inches="tight")
plt.close()

# ── 3. Force plot: compact single-prediction visualization ───────────────────
shap.force_plot(
    explainer.expected_value,
    shap_values[denied_idx].values,
    X_test.iloc[denied_idx],
    matplotlib=True,
    show=False,
)
plt.savefig("shap_force.png", dpi=150, bbox_inches="tight")
plt.close()

# ── 4. Dependence plot: feature relationship + interactions ──────────────────
# Shows how a feature's SHAP value varies with its raw value
# Color by the interacting feature (auto-detected)
shap.dependence_plot(
    "credit_util",
    shap_values.values,
    X_test,
    interaction_index="dti_ratio",   # color by interacting feature
    show=False,
)
plt.savefig("shap_dependence.png", dpi=150, bbox_inches="tight")
plt.close()

# ── 5. Bar plot: mean absolute SHAP by feature ───────────────────────────────
shap.summary_plot(shap_values, X_test, plot_type="bar", show=False)
plt.savefig("shap_bar.png", dpi=150, bbox_inches="tight")
plt.close()

print("All SHAP plots saved.")

SHAP Interaction Values

SHAP interaction values $\Phi_{ij}$ decompose each SHAP value into the feature's main effect and its interaction with every other feature:

$\phi_i(x) = \Phi_{ii}(x) + \sum_{j \neq i} \Phi_{ij}(x)$

where $\Phi_{ii}(x)$ is the main effect of feature $i$ and $\Phi_{ij}(x) = \Phi_{ji}(x)$ is the interaction between features $i$ and $j$.

# Compute SHAP interaction values (more expensive: O(TLD^3))
# Use a subset for speed
X_sample = X_test.iloc[:200]
shap_interaction = explainer.shap_interaction_values(X_sample)
# shape: (n_samples, n_features, n_features)
# shap_interaction[i, j, k] = interaction effect of features j and k on sample i

# Main effects (diagonal)
main_effects = np.array([
    shap_interaction[i].diagonal() for i in range(len(X_sample))
])  # shape: (n_samples, n_features)

# Interaction matrix: mean absolute interaction between all feature pairs
mean_interactions = np.abs(shap_interaction).mean(axis=0)  # (n_features, n_features)

# Visualize top interactions
fig, ax = plt.subplots(figsize=(10, 8))
im = ax.imshow(mean_interactions, cmap="Blues", aspect="auto")
ax.set_xticks(range(len(X.columns)))
ax.set_yticks(range(len(X.columns)))
ax.set_xticklabels(X.columns, rotation=45, ha="right", fontsize=9)
ax.set_yticklabels(X.columns, fontsize=9)
plt.colorbar(im, ax=ax, label="Mean |SHAP Interaction|")
ax.set_title("SHAP Feature Interaction Heatmap\n(diagonal = main effects, off-diagonal = interactions)")
plt.tight_layout()
plt.savefig("shap_interactions.png", dpi=150)
plt.close()

# Find the strongest interaction pair
np.fill_diagonal(mean_interactions, 0)  # remove main effects
flat_idx = np.unravel_index(np.argmax(mean_interactions), mean_interactions.shape)
feat_a, feat_b = X.columns[flat_idx[0]], X.columns[flat_idx[1]]
print(f"Strongest interaction: {feat_a} × {feat_b} = {mean_interactions[flat_idx]:.4f}")

SHAP for Feature Selection

SHAP values enable principled recursive feature elimination - remove features that contribute least to any individual prediction, not just globally.

from sklearn.model_selection import cross_val_score

def shap_recursive_feature_elimination(
    model_class,
    model_params: dict,
    X_train: pd.DataFrame,
    y_train: np.ndarray,
    X_test: pd.DataFrame,
    y_test: np.ndarray,
    min_features: int = 3,
) -> dict:
    """
    Iteratively remove the lowest-SHAP feature until performance drops.
    Returns: dict mapping n_features -> test AUC
    """
    remaining_features = list(X_train.columns)
    results = {}

    while len(remaining_features) >= min_features:
        # Train model on remaining features
        model_temp = model_class(**model_params, random_state=42, verbosity=0)
        model_temp.fit(X_train[remaining_features], y_train)

        # Compute mean absolute SHAP per feature
        explainer_temp = shap.TreeExplainer(model_temp)
        shap_vals = explainer_temp.shap_values(X_train[remaining_features])
        if isinstance(shap_vals, list):
            shap_vals = shap_vals[1]  # binary classification: use class-1 SHAP
        mean_abs = np.abs(shap_vals).mean(axis=0)

        # Test performance
        auc = roc_auc_score(y_test, model_temp.predict_proba(X_test[remaining_features])[:, 1])
        results[len(remaining_features)] = {
            "auc": auc,
            "features": remaining_features.copy(),
            "removed": None if len(remaining_features) == len(X_train.columns)
                       else remaining_features[np.argmin(mean_abs)],
        }
        print(f"Features: {len(remaining_features):2d}  AUC: {auc:.4f}  "
              f"Features: {remaining_features}")

        if len(remaining_features) <= min_features:
            break

        # Remove feature with lowest SHAP importance
        worst_idx = np.argmin(mean_abs)
        worst_feat = remaining_features[worst_idx]
        remaining_features.remove(worst_feat)

    return results

# Run RFE
rfe_results = shap_recursive_feature_elimination(
    xgb.XGBClassifier,
    {"n_estimators": 100, "max_depth": 4, "learning_rate": 0.1},
    X_train, y_train, X_test, y_test,
    min_features=3,
)

SHAP for Data Leakage Detection

A SHAP dependence plot that shows a nearly perfect linear relationship between a feature and its SHAP value is a major data leakage signal - the model has found a direct shortcut to the label through a feature that should not be available at inference time.

def detect_leakage_via_shap(
    shap_values,
    X: pd.DataFrame,
    threshold_r2: float = 0.92,
) -> pd.DataFrame:
    """
    Flag features where SHAP value is nearly perfectly predictable
    from the raw feature value. Perfect linearity often signals leakage.

    A legitimate strong predictor has non-linear SHAP patterns;
    a leaked feature often has a perfectly linear monotone SHAP curve.
    """
    from sklearn.linear_model import LinearRegression
    from sklearn.metrics import r2_score

    suspects = []
    shap_matrix = shap_values.values if hasattr(shap_values, "values") else shap_values

    for i, col in enumerate(X.columns):
        x_vals = X[col].values.reshape(-1, 1)
        y_vals = shap_matrix[:, i]

        # Check linear fit quality
        lr = LinearRegression().fit(x_vals, y_vals)
        r2 = r2_score(y_vals, lr.predict(x_vals))

        # Check rank correlation (catches monotone but non-linear leakage)
        from scipy.stats import spearmanr
        corr, _ = spearmanr(x_vals.ravel(), y_vals)

        suspects.append({
            "feature": col,
            "linear_r2": round(r2, 4),
            "rank_correlation": round(abs(corr), 4),
            "suspicious": r2 > threshold_r2 or abs(corr) > 0.98,
        })

    df_suspects = pd.DataFrame(suspects).sort_values("linear_r2", ascending=False)
    if df_suspects["suspicious"].any():
        print("WARNING: Potential data leakage detected:")
        print(df_suspects[df_suspects["suspicious"]].to_string(index=False))
    else:
        print("No obvious leakage detected by SHAP linearity test.")

    return df_suspects


leakage_report = detect_leakage_via_shap(shap_values, X_test)

SHAP for Fairness Analysis

Group-level SHAP attribution reveals whether the model is relying on protected attributes or proxies for them. This is essential for fair lending compliance.

def fairness_audit_shap(
    shap_values,
    X: pd.DataFrame,
    protected_col: str,
    group_labels: dict,
    n_top: int = 5,
) -> pd.DataFrame:
    """
    Compare mean SHAP attributions across demographic groups.
    Large differences in attribution may indicate disparate impact.
    """
    shap_matrix = shap_values.values if hasattr(shap_values, "values") else shap_values
    results = {}

    for group_val, group_label in group_labels.items():
        mask = X[protected_col] == group_val
        group_shap = shap_matrix[mask]
        results[group_label] = {
            col: group_shap[:, i].mean()
            for i, col in enumerate(X.columns)
        }

    attribution_df = pd.DataFrame(results)

    # Flag features with large attribution differences between groups
    group_names = list(group_labels.values())
    if len(group_names) == 2:
        attribution_df["difference"] = (
            attribution_df[group_names[0]] - attribution_df[group_names[1]]
        ).abs()
        attribution_df = attribution_df.sort_values("difference", ascending=False)

    return attribution_df.head(n_top)


# Example: compare SHAP attributions for different housing types
fairness_report = fairness_audit_shap(
    shap_values,
    X_test,
    protected_col="housing_type",
    group_labels={0: "lease", 1: "used", 2: "own", 3: "contract"},
    n_top=5,
)
print("\n=== Fairness Audit: SHAP Attribution by Housing Type ===")
print(fairness_report.to_string())

Adverse Action Notice Generator

def generate_adverse_action_notice(
    explainer,
    applicant: pd.Series,
    top_n: int = 3,
) -> str:
    """
    Generate a legally-informed adverse action explanation using SHAP.
    Satisfies ECOA/Reg-B requirements for adverse action letters.
    """
    shap_vals     = explainer(applicant.to_frame().T)
    contributions = pd.Series(shap_vals[0].values, index=applicant.index)

    # Positive SHAP = increased default risk = denial reason
    denial_reasons = (
        contributions[contributions > 0]
        .sort_values(ascending=False)
        .head(top_n)
    )
    mitigating = (
        contributions[contributions < 0]
        .sort_values()
        .head(top_n)
    )

    lines = [
        "ADVERSE ACTION NOTICE",
        "=" * 44,
        f"Applicant feature values: {dict(applicant.round(2))}",
        "",
        "Principal reasons for this credit decision:",
    ]
    for i, (feat, val) in enumerate(denial_reasons.items(), 1):
        lines.append(
            f"  {i}. {feat.replace('_', ' ').title()}: {applicant[feat]:.2f}  "
            f"(risk contribution: +{val:.3f})"
        )

    if not mitigating.empty:
        lines.append("")
        lines.append("Factors that partially offset this decision:")
        for feat, val in mitigating.items():
            lines.append(
                f"  - {feat.replace('_', ' ').title()}: {applicant[feat]:.2f}  "
                f"(risk reduction: {val:.3f})"
            )

    lines.append("")
    lines.append("Note: Contact us within 60 days to review this decision.")
    return "\n".join(lines)


# Generate notice for the denied applicant
notice = generate_adverse_action_notice(explainer, X_test.iloc[denied_idx])
print(notice)

Production SHAP Drift Monitoring

SHAP values have a second killer application: detecting feature-level distribution shift weeks before AUC degradation becomes measurable.

def compute_mean_abs_shap(explainer, X: pd.DataFrame) -> pd.Series:
    """Return mean |SHAP| per feature for a dataset batch."""
    sv = explainer(X)
    return pd.Series(np.abs(sv.values).mean(axis=0), index=X.columns)


def monitor_shap_drift(
    baseline_shap: pd.Series,
    current_shap: pd.Series,
    alert_threshold: float = 0.20,
) -> pd.DataFrame:
    """
    Compare SHAP importance distributions between baseline and current window.
    A relative change > alert_threshold triggers an investigation alert.

    SHAP drift fires weeks before AUC degradation becomes measurable.
    """
    report = pd.DataFrame({
        "baseline": baseline_shap,
        "current":  current_shap,
    })
    report["relative_change"] = (
        (report["current"] - report["baseline"])
        / (report["baseline"].abs() + 1e-9)
    )
    report["alert"] = report["relative_change"].abs() > alert_threshold
    return report.sort_values("relative_change", key=abs, ascending=False)


# Simulate baseline (month 1) vs current (month 6) scoring batches
baseline_shap = compute_mean_abs_shap(explainer, X_test.iloc[:300])
current_shap  = compute_mean_abs_shap(explainer, X_test.iloc[300:600])

drift_report = monitor_shap_drift(baseline_shap, current_shap, alert_threshold=0.20)
print("\n=== SHAP Drift Monitoring Report ===")
print(drift_report.to_string())

alerted_features = drift_report[drift_report["alert"]]
if len(alerted_features) > 0:
    print(f"\nALERT: {len(alerted_features)} feature(s) show >20% importance shift:")
    for feat, row in alerted_features.iterrows():
        direction = "increased" if row["relative_change"] > 0 else "decreased"
        print(f"  {feat}: {direction} by {abs(row['relative_change']):.1%}")

:::tip SHAP drift beats PSI for root cause analysis Population Stability Index (PSI) tells you something in the model's input distribution changed. SHAP drift tells you which specific feature changed importance and by how much, pointing directly at the root cause. If employment_months drops from 12% of SHAP importance to 2%, either the feature pipeline is broken or the relationship between employment tenure and default has fundamentally changed. Either way, you know exactly where to look. :::

Feature Importance Methods Comparison

Method	Scope	Speed	Handles Correlated Features	Individual Explanations	Regulatory Use
MDI split count	Global	Instant	Biased (high-cardinality)	No	No
MDI gain	Global	Instant	Biased (depth-dependent)	No	No
Permutation (individual)	Global	Slow	Splits credit (problematic)	No	No
Permutation (grouped)	Global	Moderate	Correct for groups	No	No
SHAP global (mean	φ	)	Global	Fast (TreeSHAP)	Correctly distributes
SHAP local (per-sample)	Local	Fast (TreeSHAP)	Correctly distributes	Yes	Yes
SHAP interactions	Local+Global	Slower (O(TLD³))	Yes	Yes	Advanced

:::danger Never use MDI alone for high-stakes decisions MDI feature importance is systematically biased toward high-cardinality continuous features and depth-0 root splits. Using MDI to claim "feature X is the most important predictor" in a regulatory filing or model card without cross-referencing permutation or SHAP importance can lead to materially incorrect characterizations of model behavior. In any regulated context, SHAP values are the minimum standard for defensible explanation. :::

YouTube Resources

Video	Channel	Why Watch It
SHAP Values - Main Ideas	StatQuest with Josh Starmer	Best intuitive introduction to Shapley values and SHAP
SHAP for Machine Learning - Full Tutorial	Aleksander Molnar	Deep dive into all SHAP plot types with real code
Feature Importance - MDI vs Permutation	Towards Data Science	MDI bias demonstration with real dataset
TreeSHAP - Fast Exact SHAP for Trees	Christoph Molnar	TreeSHAP algorithm explanation at the right depth
SHAP for Model Monitoring	Evidently AI	Production SHAP drift monitoring patterns

Interview Questions and Answers

Q1: Explain the MDI formula and its bias toward high-cardinality features with a concrete example.

MDI computes $\text{MDI}(f_j) = \frac{1}{T}\sum_t\sum_{n:\text{split on }j} \frac{n_t}{n} \cdot \Delta I(n)$ - the weighted average impurity reduction across all nodes that split on feature $j$. The bias: a continuous feature with 1000 unique values has 999 candidate thresholds; a binary feature has 1. The continuous feature gets far more split opportunities, accumulating MDI score simply by appearing in more nodes - even if each individual split is no better than the binary feature's single split. Concrete example: a customer ID feature (unique per row) would achieve maximum MDI because it can perfectly split the training set at every node, even though it is completely useless for generalization. In practice, continuous features like income or credit_score systematically appear more important than binary flags even when their actual predictive contribution is similar. The fix: always compare MDI against permutation importance computed on a held-out validation set.

Q2: What is TreeSHAP and why is it O(TLD²) instead of O(2^d) for naive Shapley computation?

Naive Shapley values require evaluating the model on all $2^d$ feature subsets to measure marginal contributions - exponential in $d$. TreeSHAP (Lundberg et al., 2018) exploits the tree structure: in a decision tree, the path from root to any leaf is determined by exactly which features' values are used as split criteria along the path. For any leaf, we can analytically compute how much each feature on the path contributed to reaching that leaf, and how much each feature would change the prediction if it were absent (replaced by its background distribution). The dynamic programming algorithm accumulates these contributions bottom-up in $O(LD^2)$ per tree (L leaves, D depth), and across T trees gives $O(TLD^2)$. This is polynomial, not exponential, making exact SHAP computation feasible for any tree model in production.

Q3: How do SHAP interaction values decompose feature effects, and why are they useful?

SHAP interaction values $\Phi_{ij}(x)$ decompose each feature's SHAP value into a main effect ($\Phi_{ii}$, the feature's independent contribution) and interaction effects ($\Phi_{ij}$ for $j \neq i$, how the feature's contribution changes depending on another feature's value). The decomposition satisfies $\phi_i(x) = \sum_j \Phi_{ij}(x)$. They are useful for: (1) discovering interactions: if $\Phi_{ij}$ is large, features $i$ and $j$ interact in a non-additive way - the effect of one depends on the value of the other; (2) explaining complex predictions: a prediction might be driven by an interaction between income and loan amount (the ratio matters, not either alone); (3) model debugging: unexpected strong interactions between an input feature and a protected attribute (like employment status and housing type) may indicate proxy discrimination. The cost is $O(TLD^3)$ computation - more expensive than standard SHAP but feasible for small-medium sample sizes.

Q4: How would you use SHAP to detect data leakage in a model?

Train the model and compute SHAP dependence plots for every feature. A feature with data leakage shows a nearly perfectly monotone linear relationship between its raw value and its SHAP value - because the model has found a direct shortcut to the label through a feature that should not be available at inference time. Legitimate strong predictors have more complex, non-linear SHAP patterns (the feature's effect depends on other features, on different segments of the data, etc.). Implement this formally by fitting a linear regression of SHAP values on raw feature values and checking $R^2$ - an $R^2 > 0.95$ is a red flag requiring investigation. Also check Spearman rank correlation between feature value and SHAP value - near $\pm 1$ suggests the model has learned a near-perfect monotone mapping that should not exist.

Q5: Design a SHAP-based feature selection pipeline for a production model. What are the tradeoffs vs permutation importance?

SHAP-based RFE: (1) train model on all features; (2) compute mean absolute SHAP per feature; (3) remove the feature with the lowest mean absolute SHAP; (4) repeat until performance drops below an acceptable threshold. Advantages over permutation importance RFE: (a) SHAP correctly handles correlated features by distributing credit appropriately - permutation importance falsely undervalues all features in a correlated group; (b) SHAP is faster to compute than full permutation importance (TreeSHAP is polynomial, permutation requires many model evaluations); (c) SHAP-based importance tracks which features drive individual predictions, so removing a low-SHAP feature genuinely removes local predictive contribution, not just global average contribution. Tradeoff: SHAP assumes the trained model's feature importance generalizes - if there are very many features and the model itself is poorly regularized, SHAP importance on the training model may not reflect importance on new data. Use cross-validation to validate the selected feature set.

Q6: What does the SHAP efficiency axiom mean for regulatory compliance?

The efficiency axiom states: $\sum_{i=1}^{d} \phi_i = f(x) - \mathbb{E}[f(X)]$. In plain terms, the SHAP values for all features sum exactly to the model's prediction minus the baseline (expected) prediction. This has a direct regulatory consequence: when a loan application is denied, you can present the applicant with a statement like "your predicted default risk is 0.73, compared to the baseline of 0.42. Feature A contributed +0.18, Feature B contributed +0.12, Feature C reduced risk by -0.08, and these three factors sum to exactly +0.31 = 0.73 - 0.42." This additive consistency is what makes SHAP values legally defensible under ECOA adverse action notice requirements - the explanation is mathematically guaranteed to be internally consistent, unlike qualitative explanations that can be cherry-picked or arbitrary. Other explanation methods (LIME, gradient-based saliency) do not satisfy this exact consistency guarantee.

Complete SHAP Analysis Pipeline

import shap
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_classification

def shap_full_analysis(model, X_train: np.ndarray, X_test: np.ndarray,
                         feature_names: list, max_display: int = 20):
    """
    Complete SHAP analysis pipeline:
    1. Compute SHAP values (TreeSHAP for tree models)
    2. Global importance bar plot
    3. Beeswarm summary plot
    4. SHAP dependence plots for top features
    5. Waterfall plot for individual predictions
    """
    # ── Step 1: Compute SHAP values ───────────────────────────────────────────
    explainer = shap.TreeExplainer(model)
    shap_values = explainer.shap_values(X_test)

    # For binary classification: shap_values is (n_samples, n_features)
    # For multi-class: shap_values is a list of length n_classes
    if isinstance(shap_values, list):
        sv = shap_values[1]  # positive class SHAP values
    else:
        sv = shap_values

    print(f"SHAP values computed: shape={sv.shape}")
    print(f"Expected value (baseline): {explainer.expected_value:.4f}")

    # Verify efficiency axiom: SHAP values sum to prediction - baseline
    from sklearn.base import is_classifier
    if is_classifier(model) and hasattr(model, 'predict_proba'):
        preds = model.predict_proba(X_test[:5])[:, 1]
        shap_sums = sv[:5].sum(axis=1) + explainer.expected_value
        print("\nEfficiency check (first 5 samples):")
        for i in range(5):
            print(f"  pred={preds[i]:.4f}, SHAP_sum={shap_sums[i]:.4f}, "
                  f"diff={abs(preds[i]-shap_sums[i]):.6f}")

    # ── Step 2: Global importance ──────────────────────────────────────────────
    mean_abs_shap = np.abs(sv).mean(axis=0)
    feature_importance_df = pd.DataFrame({
        'feature': feature_names,
        'mean_abs_shap': mean_abs_shap,
    }).sort_values('mean_abs_shap', ascending=False)

    fig, axes = plt.subplots(1, 2, figsize=(16, 6))

    top_feats = feature_importance_df.head(max_display)
    axes[0].barh(range(len(top_feats)), top_feats['mean_abs_shap'].values[::-1],
                  color='#3b82f6', alpha=0.8)
    axes[0].set_yticks(range(len(top_feats)))
    axes[0].set_yticklabels(top_feats['feature'].values[::-1], fontsize=9)
    axes[0].set_xlabel("Mean |SHAP value|")
    axes[0].set_title(f"Global Feature Importance (Top {max_display})", fontsize=12)

    # Beeswarm-style plot: scatter SHAP values by feature
    top_k = min(10, len(feature_names))
    top_feature_idx = np.argsort(mean_abs_shap)[-top_k:]

    for i, feat_idx in enumerate(top_feature_idx):
        feat_shap = sv[:, feat_idx]
        axes[1].scatter(feat_shap, np.full_like(feat_shap, i),
                         c=X_test[:, feat_idx], cmap='coolwarm', s=5, alpha=0.5,
                         vmin=np.percentile(X_test[:, feat_idx], 5),
                         vmax=np.percentile(X_test[:, feat_idx], 95))

    axes[1].set_yticks(range(top_k))
    axes[1].set_yticklabels([feature_names[i] for i in top_feature_idx], fontsize=9)
    axes[1].axvline(0, color='black', linewidth=0.8)
    axes[1].set_xlabel("SHAP value (impact on model output)")
    axes[1].set_title("SHAP Values Distribution\n(color = feature value: red=high, blue=low)")

    plt.tight_layout()
    plt.show()

    return sv, explainer, feature_importance_df


def shap_waterfall_single_prediction(sv_single: np.ndarray, expected_value: float,
                                      feature_names: list, feature_values: np.ndarray,
                                      prediction: float, title: str = "SHAP Explanation"):
    """
    Waterfall plot for a single prediction: shows how each feature
    pushes the model output from the baseline to the final prediction.
    """
    # Sort by absolute SHAP value
    idx = np.argsort(np.abs(sv_single))[::-1][:10]   # top 10
    sorted_shap  = sv_single[idx]
    sorted_names = [f"{feature_names[i]}={feature_values[i]:.3g}" for i in idx]

    colors = ['#ef4444' if v > 0 else '#3b82f6' for v in sorted_shap]

    fig, ax = plt.subplots(figsize=(10, 6))

    # Cumulative SHAP contributions
    cumulative = expected_value
    for i, (shap_val, name, color) in enumerate(zip(sorted_shap[::-1],
                                                      sorted_names[::-1], colors[::-1])):
        ax.barh(i, shap_val, left=cumulative, color=color, alpha=0.8, edgecolor='white')
        ax.text(cumulative + shap_val/2, i,
                f"{shap_val:+.3f}", ha='center', va='center', fontsize=8, color='white',
                fontweight='bold')
        cumulative += shap_val

    ax.set_yticks(range(len(sorted_names)))
    ax.set_yticklabels(sorted_names[::-1], fontsize=9)
    ax.axvline(expected_value, color='gray', linestyle='--', linewidth=1,
                label=f'Baseline={expected_value:.3f}')
    ax.axvline(prediction, color='black', linestyle='-', linewidth=2,
                label=f'Prediction={prediction:.3f}')
    ax.set_xlabel("Model Output")
    ax.set_title(title)
    ax.legend()
    plt.tight_layout()
    plt.show()


# Run full analysis
X, y = make_classification(n_samples=5000, n_features=15,
                             n_informative=8, n_redundant=3, random_state=42)
feature_names = [f"feature_{i}" for i in range(15)]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

rf = RandomForestClassifier(n_estimators=200, max_depth=8, random_state=42)
rf.fit(X_train, y_train)

sv, explainer, importance_df = shap_full_analysis(rf, X_train, X_test, feature_names)

# Explain one specific prediction
sample_idx = 42
shap_waterfall_single_prediction(
    sv[sample_idx],
    explainer.expected_value,
    feature_names,
    X_test[sample_idx],
    rf.predict_proba(X_test[[sample_idx]])[0, 1],
    title=f"Explanation for Sample {sample_idx} - Prediction: "
          f"{rf.predict_proba(X_test[[sample_idx]])[0, 1]:.3f}"
)

SHAP Dependence Plots: Discovering Interactions

SHAP dependence plots show the relationship between a feature's value and its SHAP contribution, colored by the value of the most interacting feature:

import shap
import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import make_classification

X, y = make_classification(n_samples=3000, n_features=10,
                             n_informative=6, n_redundant=2, random_state=42)
feature_names = [f"feature_{i}" for i in range(10)]

rf = RandomForestClassifier(n_estimators=200, random_state=42, n_jobs=-1)
rf.fit(X, y)

explainer   = shap.TreeExplainer(rf)
shap_values = explainer.shap_values(X)
sv = shap_values[1] if isinstance(shap_values, list) else shap_values

def shap_dependence_analysis(sv: np.ndarray, X: np.ndarray,
                               feature_names: list, top_k: int = 4):
    """
    Plot SHAP dependence plots for top-k features.
    Color by the feature that most interacts with each target feature.

    The interaction feature is identified by the strongest color gradient
    (shap library does this automatically via 'auto' interaction index).
    """
    mean_abs = np.abs(sv).mean(axis=0)
    top_features = np.argsort(mean_abs)[-top_k:][::-1]

    fig, axes = plt.subplots(1, top_k, figsize=(5*top_k, 5))

    for ax, feat_idx in zip(axes, top_features):
        feat_name = feature_names[feat_idx]

        # Find most interacting feature (maximize |covariance of SHAP with other feature|)
        shap_col = sv[:, feat_idx]
        interact_scores = [
            abs(np.corrcoef(shap_col, X[:, j])[0, 1])
            for j in range(X.shape[1]) if j != feat_idx
        ]
        interact_idx = [j for j in range(X.shape[1]) if j != feat_idx][np.argmax(interact_scores)]
        interact_name = feature_names[interact_idx]

        scatter = ax.scatter(X[:, feat_idx], shap_col,
                              c=X[:, interact_idx], cmap='coolwarm', s=8, alpha=0.5)
        plt.colorbar(scatter, ax=ax, label=interact_name, fraction=0.046)

        # Smooth trend line
        order = np.argsort(X[:, feat_idx])
        from scipy.ndimage import uniform_filter1d
        x_sorted = X[order, feat_idx]
        y_sorted = shap_col[order]
        y_smooth = uniform_filter1d(y_sorted, size=len(y_sorted)//20 + 1)
        ax.plot(x_sorted, y_smooth, color='black', linewidth=2, alpha=0.8)

        ax.axhline(0, color='gray', linewidth=0.8, linestyle='--')
        ax.set_xlabel(feat_name)
        ax.set_ylabel(f"SHAP({feat_name})")
        ax.set_title(f"SHAP Dependence\n({feat_name} × {interact_name})", fontsize=9)

    plt.suptitle("SHAP Dependence Plots (color = top interacting feature)", fontsize=11)
    plt.tight_layout()
    plt.show()

shap_dependence_analysis(sv, X, feature_names, top_k=4)

SHAP for Model Comparison and Selection

SHAP can be used to compare two models side-by-side - not just their accuracy but the explanations they produce:

import shap
import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=5000, n_features=15,
                             n_informative=8, n_redundant=3, random_state=42)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=42)
feature_names = [f"f{i}" for i in range(15)]

# Train two models to compare
rf_model  = RandomForestClassifier(n_estimators=200, random_state=42).fit(X_tr, y_tr)
gbm_model = GradientBoostingClassifier(n_estimators=200, learning_rate=0.05,
                                        random_state=42).fit(X_tr, y_tr)

# SHAP for both
exp_rf  = shap.TreeExplainer(rf_model)
exp_gbm = shap.TreeExplainer(gbm_model)

sv_rf  = exp_rf.shap_values(X_te)
sv_rf  = sv_rf[1] if isinstance(sv_rf, list) else sv_rf
sv_gbm = exp_gbm.shap_values(X_te)
sv_gbm = sv_gbm[1] if isinstance(sv_gbm, list) else sv_gbm

# Feature importance comparison
importance_rf  = np.abs(sv_rf).mean(axis=0)
importance_gbm = np.abs(sv_gbm).mean(axis=0)

# Rank correlation - do models agree on what matters?
from scipy.stats import spearmanr
rank_corr, _ = spearmanr(importance_rf, importance_gbm)
print(f"RF vs GBM SHAP importance rank correlation: {rank_corr:.4f}")
print("  > 0.8: models agree - either can be used with confidence")
print("  < 0.6: models disagree - investigate discrepant features")

fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(importance_rf, importance_gbm, s=40, alpha=0.8, color='#7c3aed')

for i, name in enumerate(feature_names):
    if importance_rf[i] > np.percentile(importance_rf, 75) or \
       importance_gbm[i] > np.percentile(importance_gbm, 75):
        ax.annotate(name, (importance_rf[i], importance_gbm[i]), fontsize=8)

ax.set_xlabel("RF mean |SHAP|")
ax.set_ylabel("GBM mean |SHAP|")
ax.set_title(f"RF vs GBM Feature Importance Agreement\n(Spearman r={rank_corr:.3f})")
ax.axline((0,0), slope=1, color='gray', linestyle='--', alpha=0.5, label='Perfect agreement')
ax.legend()
plt.tight_layout()
plt.show()

:::tip 🎮 Interactive Playground

Visualize this concept: Try the SHAP Values & Feature Attribution demo on the EngineersOfAI Playground - no code required.

:::