SHAP Values: Why Your Model Predicted That

The Black Box Problem: When Accuracy Isn’t Enough

You’ve spent weeks tuning your model. You cross-validated it, tweaked the hyperparameters, and finally hit 94% accuracy on the test set. Your team is thrilled. Your boss wants it in production.

Then someone asks the question that catches you off guard. Mr. Owens, the lending team’s manager, leans across the table: “Why did it deny Maria’s loan application?”

You look at your model. It says no. But why it says no — that’s where things get fuzzy. You can see the features that went in: age, income, credit score, employment history. You just can’t point to any single one and say, “This is the reason.” The model is a black box. Maria is still waiting for an answer, and Mr. Owens needs one before he faces her — or, worse, a regulator.

This is where real-world machine learning gets hard. Accuracy is where it starts. What matters is trust, and trust comes from understanding.

Why This Matters More Than You Think

Picture a doctor with a perfect diagnosis rate who can’t explain why she thinks you have a disease. She just says, “You do.” Would you trust her? Most people would get a second opinion.

The same logic applies to your model. Here’s what’s actually going on:

Regulators now demand explanations. The European Union’s General Data Protection Regulation (GDPR) requires that automated decisions affecting people — like loan denials or job rejections — be explainable. If your model rejects someone, you need to tell them why. “The algorithm said so” isn’t good enough anymore.

Your model might be cheating. High accuracy can hide a multitude of sins. A model might learn to predict based on proxy variables — features that correlate with what you want to predict but for the wrong reasons. Say a hiring model trained on historical data learns that certain zip codes predict “good employee.” What it’s really doing is replicating past hiring bias. The model works, but it’s unfair.

You can’t debug what you can’t see. When your model fails on a new customer or a new market, you need to understand why it failed. Could be a data quality issue, a feature that no longer predicts, or a hidden bias. Without explanations, you’re flying blind.

Your stakeholders need to believe in it. A loan officer, a doctor, a hiring manager — these people won’t use a model they don’t understand. They’ll override it, ignore it, or worse, use it as cover for their own biases (“The algorithm said so”). Real adoption requires real understanding.

What We’re Going to Learn

We’re going to solve this problem using SHAP values — a tool that tells you, for any single prediction, which features pushed the model toward that decision and by how much.

What makes SHAP worth using is its game-theory foundation, which makes it mathematically principled. We won’t need to understand the math, though. We’ll build intuition first, then look at the code, then interpret what the numbers actually mean.

By the end, you’ll be able to:

• Explain to your boss exactly why your model made a specific prediction
• Spot when your model is learning the wrong patterns
• Debug model failures by seeing which features drove each decision
• Build stakeholder trust by showing, not just telling

Let’s start with a concrete example.

Understanding the Intuition: A Game Theory Thought Experiment

The Poker Game Analogy

Picture a poker game with three friends. You’re all betting on who has the best hand. The round ends. You win — the pot is yours. Your friends ask: “Why did you win?”

“I had a good hand” works as an answer, but it’s vague. What if you broke it down?

• Without your ace of spades, you would’ve lost. That card was worth, say, 40% of your win.
• Your pair of kings added another 35%.
• Your position at the table — betting last — contributed 25%.

Add those up and you get 100% of your win. Each share is fair — it accounts for what that feature contributed in the context of all the other features.

That’s the core intuition behind SHAP values. Rather than asking, “Which features matter?” we ask, “For this specific prediction, how much did each feature contribute?”

Why This Is Hard (And Why SHAP Solves It)

The tricky part is splitting credit fairly. Say your model has three features: age, income, and credit score. They all matter. But how much did each one contribute to this particular prediction?

You might think: just check the feature importance. But feature importance tells you which features matter on average across all predictions. It says nothing about this specific prediction.

Or you might turn to coefficients from a linear regression. But your model isn’t linear — it’s a random forest, a neural network, or a gradient boosting model. Coefficients don’t apply.

Or you might remove each feature one at a time and watch how the prediction shifts. But the order matters. Remove income first, and credit score might look more important. Remove credit score first, and income might too. It’s ambiguous.

Game theorists already tackled this. In poker, there’s a fair way to split credit called the Shapley value. It works by imagining all possible orderings of features joining the game, then averaging each feature’s contribution across those orderings. The math is complex, but the intuition is simple: it’s the only way to split credit that’s fair to everyone.

SHAP values apply this idea to machine learning. For each prediction, SHAP calculates how much each feature contributed — in a way that’s mathematically fair and consistent.

The Shapley value for feature $i$ is defined as:

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

In plain English: for every possible subset $S$ of features that doesn’t include feature $i$, measure how much the model’s prediction changes when you add feature $i$ to that subset. Weight each measurement by how many orderings produce that subset, then average everything together.

Plain English	Statistical symbol	Python equivalent
The Shapley value (credit assigned to feature $i$)	$\phi_i$	shap_values[customer_idx, i, 1]
The set of all features	$F$	X_test.columns
A coalition (subset of features, not including $i$)	$S \subseteq F \setminus \{i\}$	All feature subsets iterated internally by shap.TreeExplainer
Model prediction using only features in coalition $S$	$f_S(x_S)$	Approximated via background dataset passed to TreeExplainer
Marginal contribution of adding feature $i$ to coalition $S$	$f_{S \cup \{i\}}(x_{S \cup \{i\}}) - f_S(x_S)$	The difference SHAP averages across all coalitions
Weighting factor (how many orderings produce this coalition)	$\frac{|S|!,(	F
Baseline prediction (average model output)	$\mathbb{E}[f(x)]$	explainer.expected_value[1]

The additive property: SHAP values satisfy $f(x) = \mathbb{E}[f(x)] + \sum_{i=1}^{|F|} \phi_i$. This means the base value plus all SHAP values always sums exactly to the model’s prediction — which is why the force plot is a valid decomposition, not an approximation.

What This Means in Practice

Say your model predicts a customer has a 75% chance of buying a product. The baseline — what you’d predict with no information — is 50%. So the model lifted the prediction by 25 percentage points.

SHAP values break that 25-point increase apart: 10 points from the customer’s high income, 8 from their young age, 5 from their high engagement score, 2 from their account tenure. The rest is noise.

Now you can explain the prediction. You can also catch problems this way. If you expected engagement to matter more than tenure but SHAP says the opposite, something may be off with how you engineered tenure. Or engagement score might be a proxy for something else.

Here it is in code.

Building Your First SHAP Explanation

Setting Up the Data and Model

We’ll use a real dataset: customer churn. We’ll build a model to predict whether a customer will leave, then use SHAP to explain a specific prediction.

First, let’s load the data and build a simple model:

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.preprocessing import LabelEncoder
import shap
import matplotlib.pyplot as plt

# Load a sample dataset (we'll create a synthetic one for reproducibility)
np.random.seed(42)
n_samples = 1000

# Create synthetic customer data
data = pd.DataFrame({
    'age': np.random.randint(18, 80, n_samples),
    'monthly_charges': np.random.uniform(20, 150, n_samples),
    'tenure_months': np.random.randint(0, 72, n_samples),
    'total_charges': np.random.uniform(100, 10000, n_samples),
    'num_services': np.random.randint(1, 8, n_samples),
})

# Create a target variable (churn) with some realistic patterns
data['churn'] = (
    (data['tenure_months'] < 12).astype(int) * 0.5 +  # New customers churn more
    (data['monthly_charges'] > 100).astype(int) * 0.3 +  # High charges increase churn
    (data['num_services'] < 2).astype(int) * 0.2 +  # Few services = higher churn
    np.random.uniform(0, 0.3, n_samples)  # Add some randomness
) > 0.5

data['churn'] = data['churn'].astype(int)

print(f"Dataset shape: {data.shape}")
print(f"Churn rate: {data['churn'].mean():.1%}")
print(f"\nFirst few rows:")
print(data.head())

• np.random.seed(42) — Sets the random number generator to a fixed starting point so the synthetic data comes out the same every time you run the notebook. Change this number and you get a different (but equally valid) dataset.
• pd.DataFrame({...}) — Constructs a table directly from a Python dictionary; each key becomes a column name, each value becomes the column’s data.
• np.random.randint(18, 80, n_samples) — Draws n_samples integers uniformly from the half-open interval [18, 80). The upper bound is excluded, so the oldest simulated customer is 79.
• np.random.uniform(20, 150, n_samples) — Draws floats uniformly between 20 and 150. Unlike randint, these are continuous values.
• (data['tenure_months'] < 12).astype(int) — The comparison produces a boolean Series (True/False); .astype(int) converts it to 1s and 0s so it can participate in arithmetic.
• * 0.5 — Weights that binary indicator: new customers (tenure < 12 months) contribute 0.5 to the churn score, making them the strongest churn signal.
• > 0.5 — Thresholds the weighted sum into a binary label. Customers whose combined score exceeds 0.5 are labelled as churned.
• data['churn'].mean() — Because churn is 0/1, the mean equals the proportion of churned customers — a quick sanity check on class balance. This dataset comes out to 30.5%.

This gives us a dataset of 1,000 customers. Each row has features like age, how much they pay per month, how long they’ve been a customer, and how many services they use. The target is whether they churned (left) or not.

Now let’s train a model:

# Split into train and test
X = data.drop('churn', axis=1)
y = data['churn']

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# Train a random forest
model = RandomForestClassifier(n_estimators=100, random_state=42, max_depth=10)
model.fit(X_train, y_train)

# Check accuracy
train_score = model.score(X_train, y_train)
test_score = model.score(X_test, y_test)

print(f"Train accuracy: {train_score:.1%}")
print(f"Test accuracy: {test_score:.1%}")

• data.drop('churn', axis=1) — Returns a new DataFrame with the churn column removed. axis=1 means “drop a column” (as opposed to axis=0, which drops a row).
• train_test_split(X, y, test_size=0.2, random_state=42) — Shuffles the data and reserves 20% for testing. random_state=42 makes the split reproducible. Returns four objects: training features, test features, training labels, test labels — in that order.
• RandomForestClassifier(n_estimators=100, random_state=42, max_depth=10) — Creates (but does not yet train) a random forest with 100 decision trees, each capped at depth 10 to prevent overfitting. random_state seeds the internal bootstrapping so results are reproducible.
• model.fit(X_train, y_train) — Actually trains the forest: each tree is grown on a bootstrap sample of the training data, and features are randomly subsampled at each split.
• model.score(X_train, y_train) — Computes classification accuracy (fraction of correct predictions) on the training set. Comparing this to test_score is a quick overfitting check: a large gap signals the model has memorised training data.

The model hits 98.8% accuracy on the training set and 89.5% on the test set. That gap is worth noticing — it’s a sign the forest is fitting some training-set noise it can’t fully generalize — but 89.5% is still a solid working model, and it’s what the rest of this article uses. Either way, accuracy alone tells us nothing about why it makes each prediction.

Creating SHAP Explanations

Here’s where SHAP comes in. We’ll build a SHAP explainer to calculate how much each feature contributes to each prediction:

# Create a SHAP explainer
explainer = shap.TreeExplainer(model)

# Calculate SHAP values for the test set
shap_values = explainer.shap_values(X_test)

print(f"SHAP values shape: {np.array(shap_values).shape}")
print(f"Number of features: {len(X_test.columns)}")
print(f"Number of test samples: {len(X_test)}")

• shap.TreeExplainer(model) — Creates a SHAP explainer specialised for tree-based models (random forests, gradient boosting, XGBoost, LightGBM, etc.). It uses an exact, polynomial-time algorithm that exploits the tree structure — much faster than the general KernelExplainer, which uses sampling. You pass the trained model object directly; no background data is required for TreeExplainer — it derives its baseline from the tree structure itself (the “tree-path” baseline), not from a supplied sample.
• explainer.shap_values(X_test) — Runs the TreeSHAP algorithm over every row in X_test. On a current shap install this returns a single NumPy array shaped (n_test_samples, n_features, n_classes). Older versions returned a list of two arrays, one per class — code written against that older shape (shap_values[1], shap_values[0]) will silently pick the wrong axis on the current shape rather than the class it means to.
• np.array(shap_values).shape — Prints the shape directly. You’ll see (200, 5, 2) — 200 test rows, 5 features, 2 classes.
• len(X_test.columns) — The number of feature columns; useful as a quick sanity check that the SHAP array’s middle dimension matches.

So what just happened? The TreeExplainer inspected your random forest model and computed, for every prediction in the test set, how much each feature contributed. The output is a set of SHAP values — one per feature, per prediction, per class.

For a classification model like this, SHAP gives you contributions toward both classes at once. Since this is binary classification (churn or not churn), the last axis of the array has length 2. We’ll focus on the positive class (churn = 1), which is the last axis, not the first:

# For binary classification, the last axis of shap_values indexes the class
# shap_values[:, :, 1] is for the positive class (churn = 1)
shap_values_churn = shap_values[:, :, 1]

print(f"SHAP values for churn class shape: {shap_values_churn.shape}")
print(f"This means: {shap_values_churn.shape[0]} predictions, {shap_values_churn.shape[1]} features")

• shap_values[:, :, 1] — Selects every row, every feature, and class index 1 (churn) from the last axis. Index [:, :, 0] would give contributions toward the negative class (no churn). For binary classification these are mirror images on the class axis: shap_values[:, :, 0] + shap_values[:, :, 1] equals zero for every row-feature cell, because each class’s contribution is measured as a push away from the other.
• shap_values_churn.shape — A tuple (n_rows, n_features). Printing .shape[0] and .shape[1] separately makes the two dimensions explicit for readers who are new to NumPy array indexing.

Explaining a Single Prediction

Pick one customer and see why the model made its prediction:

# Pick the first customer in the test set
customer_idx = 0
customer = X_test.iloc[customer_idx]

print(f"Customer features:")
print(customer)
print(f"\nModel prediction (probability of churn): {model.predict_proba(X_test.iloc[[customer_idx]])[0][1]:.1%}")
print(f"Actual churn: {y_test.iloc[customer_idx]}")

• X_test.iloc[customer_idx] — Uses integer-position indexing (.iloc) to grab the row at position 0. This returns a pandas Series with feature names as the index — convenient for printing.
• X_test.iloc[[customer_idx]] — The double brackets [[...]] return a single-row DataFrame rather than a Series. predict_proba requires a 2-D input, so this distinction matters.
• model.predict_proba(...)[0][1] — predict_proba returns an array of shape (n_rows, n_classes). [0] selects the first (and only) row; [1] selects the probability of class 1 (churn). [0] would give the probability of no-churn.
• y_test.iloc[customer_idx] — Retrieves the actual label for this customer from the test set, so we can compare the model’s prediction to ground truth.

This customer is 42 years old, pays $89.99/month, has been a customer for 69 months, has spent $640.88 in total, and uses 3 services. The model gives them a 4.1% chance of churning — and they didn’t churn, so that’s a correct call on a genuinely loyal-looking customer.

Now the SHAP values for that same customer:

# Get SHAP values for this customer
customer_shap = shap_values_churn[customer_idx]

print(f"\nSHAP values for this customer:")
for feature, shap_val in zip(X_test.columns, customer_shap):
    print(f"  {feature}: {shap_val:.4f}")

# The base value is the model's average prediction
base_value = explainer.expected_value[1]
print(f"\nBase value (average prediction): {base_value:.1%}")
print(f"Sum of SHAP values: {customer_shap.sum():.4f}")
print(f"Final prediction: {base_value + customer_shap.sum():.1%}")

# Don't just print the check - enforce it
predicted_proba = model.predict_proba(X_test.iloc[[customer_idx]])[0][1]
assert abs(base_value + customer_shap.sum() - predicted_proba) < 1e-6

• shap_values_churn[customer_idx] — Selects the row corresponding to our customer from the SHAP value array. The result is a 1-D NumPy array of length n_features, one value per feature.
• zip(X_test.columns, customer_shap) — Pairs each feature name with its corresponding SHAP value so we can print them together in a loop.
• explainer.expected_value[1] — The model’s average predicted probability of churn across the training data. This is the “baseline” from which all SHAP values are measured. Index [1] selects the positive class, matching shap_values[:, :, 1].
• customer_shap.sum() — Sums all five SHAP values. By the additive property of Shapley values, base_value + customer_shap.sum() must equal model.predict_proba(...)[0][1] exactly (up to floating-point precision).
• assert abs(base_value + customer_shap.sum() - predicted_proba) < 1e-6 — Turns the additive property from something printed and eyeballed into something the code actually enforces. If the indexing above were wrong, this line would fail loudly instead of quietly printing a plausible-looking number.
• :.4f / :.1% — Format specifiers: .4f prints four decimal places; .1% multiplies by 100 and appends a % sign with one decimal place.

Running this prints:

age: 0.0076
monthly_charges: -0.1250
tenure_months: -0.1506
total_charges: 0.0272
num_services: -0.0326

Base value (average prediction): 31.5%
Sum of SHAP values: -0.2735
Final prediction: 4.1%

Here’s how to read this. The model starts from a baseline — the average churn rate across all customers, just over 30% for this dataset. Each feature then pushes that baseline up or down for this particular customer. This customer’s long tenure (69 months) is the single biggest factor pushing the prediction down, followed closely by their below-average monthly bill. Their total spend and age nudge slightly upward, and their service count pulls slightly down again. The final prediction — 4.1% — and the additivity check both land exactly where the assertion says they must.

Let’s visualize it:

# Create a force plot for this customer
shap.force_plot(
    explainer.expected_value[1],
    customer_shap,
    customer,
    matplotlib=True,
    show=False
)
plt.tight_layout()
plt.savefig('shap_force_plot.png', dpi=100, bbox_inches='tight')
plt.show()

print("Force plot saved!")

• shap.force_plot(...) — Draws the classic SHAP force plot: a horizontal bar starting at the base value, with red arrows pushing right (toward higher churn probability) and blue arrows pushing left (toward lower probability). The final position of the arrow tip is the model’s prediction.
• explainer.expected_value[1] — Passed as the first argument; this anchors the left edge of the plot at the baseline probability.
• customer_shap — The array of SHAP values for this customer; determines the length and direction of each colored bar.
• customer — The feature values (a pandas Series); these are printed as labels on each bar so you can see what value of each feature drove the contribution.
• matplotlib=True — Forces the plot to render as a static Matplotlib figure rather than the default interactive JavaScript widget. Required when saving to a file or running outside a Jupyter notebook.
• show=False — Prevents SHAP from calling plt.show() internally, giving us control over when the figure is displayed (useful when you want to call plt.tight_layout() or plt.savefig() first).
• plt.savefig('shap_force_plot.png', dpi=100, bbox_inches='tight') — Saves the figure to disk. dpi=100 sets resolution; bbox_inches='tight' trims excess whitespace around the plot.

The force plot is a horizontal bar chart with the baseline on the left. Red bars push the prediction to the right, toward churn, while blue bars push it left, away from churn. Each bar’s length shows the magnitude of that feature’s contribution.

Interpreting the Numbers

The pattern above is real data from a real customer, but it’s worth walking through a cleaner, rounder illustration to build intuition before moving on. Suppose a different customer’s output looked like this:

• Base value: 0.28 (28% churn rate on average)
• tenure_months: -0.12 (this customer has been with us for 60 months, which reduces churn risk by 12 percentage points)
• monthly_charges: +0.08 (they pay $120/month, which increases churn risk by 8 percentage points)
• num_services: -0.05 (they use 5 services, which reduces churn risk by 5 percentage points)
• age: +0.02 (they’re 35 years old, which slightly increases churn risk)
• total_charges: -0.03 (they’ve spent a lot total, which slightly reduces churn risk)

Final prediction: 0.28 - 0.12 + 0.08 - 0.05 + 0.02 - 0.03 = 0.18 (18% churn risk)

From there you can explain the prediction: “This customer has a low churn risk because they’ve been loyal for 5 years and use multiple services. Their high monthly bill is a concern, though. Overall, we predict an 18% chance they’ll leave.”

You can hand that to your boss. A business user can follow it.

Understanding Feature Importance at Scale

From One Prediction to Many

Explaining one customer is useful, but what about the whole model — which features matter most across all predictions?

SHAP lets us aggregate values across every prediction to see what drives the model’s decisions overall.

# Create a summary plot showing which features matter most
shap.summary_plot(shap_values_churn, X_test, plot_type="bar", show=False)
plt.tight_layout()
plt.savefig('shap_summary_bar.png', dpi=100, bbox_inches='tight')
plt.show()

print("Summary bar plot saved!")

• shap.summary_plot(shap_values_churn, X_test, plot_type="bar") — Aggregates SHAP values across all test-set predictions and draws a horizontal bar chart. Each bar’s length equals the mean absolute SHAP value for that feature: $\frac{1}{n}\sum_{j=1}^{n}|\phi_i^{(j)}|$. Features are sorted from most to least important.
• plot_type="bar" — Selects the aggregated bar chart. The default (plot_type="dot") shows a beeswarm plot where each dot is one prediction, colored by feature value — more information-dense but harder to read at a glance.
• shap_values_churn — The (n_test, n_features) array for the positive class. Passing the full array (not a single row) is what triggers the aggregation.
• X_test — Passed alongside SHAP values so the plot can label the x-axis with actual feature names from the DataFrame columns.
• show=False / plt.savefig(...) — Same pattern as the force plot: suppress auto-display so we can save cleanly.

This bar plot shows the average absolute SHAP value for each feature. Think of it as the average amount each feature moves the needle on predictions.

If tenure_months has the highest bar, tenure is the most important driver of churn predictions. If num_services has the lowest, that feature barely moves the needle.

The key distinction: this isn’t the same as traditional feature importance. The kind you get from model.feature_importances_ in scikit-learn tells you which features the model relies on most. SHAP tells you which features change the prediction most. Related ideas, but not identical.

SHAP vs. Other Feature-Importance Methods

When someone asks “which features matter?”, there are four common answers. Here’s when each one wins — and when it misleads.

Method	What it measures	Scope	Key strength	Key weakness
SHAP values	Marginal contribution of each feature, averaged over all coalitions	Per-prediction and global	Mathematically fair; local + global; works on any model	Slower than alternatives for large datasets; KernelExplainer can be very slow on non-tree models
Tree / MDI feature importance (model.feature_importances_)	How often a feature is used in splits, weighted by impurity reduction	Global only	Instant (computed during training); no extra library needed	Biased toward high-cardinality features; gives no per-prediction insight; can be misleading when features are correlated
Permutation importance (sklearn.inspection.permutation_importance)	Drop in model score when a feature’s values are randomly shuffled	Global only	Model-agnostic; measures actual impact on performance metric	Unstable when features are correlated (shuffling one correlated feature can be “rescued” by another); slow on large datasets
LIME (Local Interpretable Model-agnostic Explanations)	Coefficients of a locally-fitted linear model around one prediction	Per-prediction only	Fast; model-agnostic; intuitive linear explanation	Explanations can be unstable (different runs → different results); no global view; the local linear approximation may be poor in high-curvature regions

When to reach for SHAP

• ✅ You need to explain a specific decision to a regulator, customer, or auditor (e.g., why Maria was denied).
• ✅ You want a consistent, globally-coherent story: the same framework that explains one prediction also ranks features across all predictions.
• ✅ Your model is a tree-based ensemble (random forest, XGBoost, LightGBM) — TreeExplainer is exact and fast.
• ✅ You suspect proxy variables or bias and need to audit individual predictions.

When to be cautious with SHAP

• ⚠️ Non-tree models at scale: KernelExplainer (for neural nets, SVMs, etc.) uses sampling and can take minutes per batch. Consider LIME for a quick first pass.
• ⚠️ Highly correlated features: SHAP distributes credit across correlated features in a way that can make each one look less important than it “really” is — which is actually the correct mathematical answer, but can surprise stakeholders expecting one clear winner.
• ⚠️ You only need a global ranking, fast: If you just want to know which features to drop before retraining, permutation importance is cheaper and often good enough.

The short version

Use SHAP when you need to explain this prediction. Use permutation importance when you need a quick global audit. Use MDI importance only as a rough first look. Save LIME for when you need a fast, model-agnostic local explanation — which is exactly what the next article covers.

Seeing the Full Picture: Dependence Plots

Here’s a more useful question: for each feature, how does its SHAP value shift as the feature value itself changes?

# Create a dependence plot for tenure
shap.dependence_plot(
    "tenure_months",
    shap_values_churn,
    X_test,
    show=False
)
plt.tight_layout()
plt.savefig('shap_dependence_tenure.png', dpi=100, bbox_inches='tight')
plt.show()

print("Dependence plot for tenure saved!")

• shap.dependence_plot("tenure_months", shap_values_churn, X_test) — Draws a scatter plot with tenure_months on the x-axis and each customer’s SHAP value for tenure_months on the y-axis. Each point is one test-set customer.
• "tenure_months" — The feature whose relationship with its own SHAP value you want to visualise. You can also pass an integer column index instead of a string name.
• Automatic interaction coloring — SHAP automatically selects a second feature (the one most correlated with tenure_months in the SHAP value space) and colors each dot by that feature’s value. A color gradient that splits the scatter cloud reveals an interaction effect: the impact of tenure depends on the value of that second feature.
• show=False / plt.savefig(...) — Same save-before-show pattern used throughout.

The x-axis shows customer tenure; the y-axis shows the SHAP value — how much tenure pushed the prediction toward churn. The pattern isn’t a smooth gradient — it’s a cliff. Averaged by tenure band, the mean SHAP value is +0.659 below 12 months, then drops to -0.147 from 12–30 months, -0.145 from 30–50 months, and -0.137 above 50 months. Nearly all of the effect happens at the twelve-month mark; past that cliff, the line is close to flat, not steadily descending. That’s not a coincidence — the data-generating process above uses a hard threshold, (tenure_months < 12), not a continuous function of tenure, so there’s no gradient in the underlying data for a gradient in the plot to pick up.

Each dot’s color represents another feature, automatically chosen as the one most correlated with tenure. A visible color pattern signals an interaction, meaning the effect of tenure depends on that second feature.

If you want to see a genuine gradient in a dependence plot rather than a cliff-then-flat pattern, you’d need to generate tenure’s effect on the label as a continuous function of tenure instead of a threshold — otherwise the honest reading of this plot is “there’s a sharp cutoff at 12 months, and not much changes with tenure on either side of it.”

Spotting Problems: When SHAP Reveals Bias

The Hidden Proxy Variable

Here’s where SHAP becomes a debugging tool. Say you build a hiring model to predict candidate success. It hits 85% accuracy. Then someone asks, “Is the model biased by gender?”

You check the feature list. Gender isn’t in the model. You think you’re safe.

Let’s use SHAP to look closer:

# Hypothetical: suppose we had a hiring dataset with a hidden bias
# Let's create a synthetic example

np.random.seed(42)
n_candidates = 500

hiring_data = pd.DataFrame({
    'years_experience': np.random.uniform(0, 20, n_candidates),
    'test_score': np.random.uniform(40, 100, n_candidates),
    'years_at_last_job': np.random.uniform(0, 15, n_candidates),
    'has_advanced_degree': np.random.randint(0, 2, n_candidates),
})

# Create a target: success (but with hidden bias)
hidden_gender = np.random.randint(0, 2, n_candidates)

# Give years_at_last_job a real, modest correlation with hidden_gender.
# In the historical data this simulates, tenure records for women show more
# career gaps -- so years_at_last_job becomes a genuine (if imperfect) proxy
# for the protected attribute, the way the article's own prose suggests below.
hiring_data['years_at_last_job'] = np.clip(
    hiring_data['years_at_last_job'] - hidden_gender * 7 + np.random.normal(0, 1, n_candidates),
    0, 15
)

hiring_data['success'] = (
    (hiring_data['years_experience'] > 5).astype(int) * 0.4 +
    (hiring_data['test_score'] > 70).astype(int) * 0.4 +
    (hiring_data['has_advanced_degree'] == 1).astype(int) * 0.2 +
    (hidden_gender == 0).astype(int) * 0.3 +  # Men get a boost
    np.random.uniform(-0.2, 0.2, n_candidates)
) > 0.5

hiring_data['success'] = hiring_data['success'].astype(int)

print(f"Success rate overall: {hiring_data['success'].mean():.1%}")
print(f"Success rate for men (hidden_gender=0): {hiring_data[hidden_gender == 0]['success'].mean():.1%}")
print(f"Success rate for women (hidden_gender=1): {hiring_data[hidden_gender == 1]['success'].mean():.1%}")
print(f"corr(hidden_gender, years_at_last_job): {np.corrcoef(hidden_gender, hiring_data['years_at_last_job'])[0, 1]:.3f}")

• hidden_gender = np.random.randint(0, 2, n_candidates) — Creates a binary gender variable (0 = man, 1 = woman) that is never added to hiring_data. It exists only in this Python variable, simulating a protected attribute that was present when historical labels were created but is absent from the feature matrix the model sees.
• hiring_data['years_at_last_job'] = np.clip(... - hidden_gender * 7 + noise, 0, 15) — This is the proxy mechanism. years_at_last_job starts out independent of gender, then gets shifted down for women (hidden_gender == 1) and clipped back into the original [0, 15] range. The result: years_at_last_job is now a visible feature that correlates with the hidden attribute, exactly the kind of proxy the article’s own hiring-discrimination story depends on. Without this step, years_at_last_job would be pure noise with respect to gender, and there would be nothing for SHAP — or anything else — to detect.
• (hidden_gender == 0).astype(int) * 0.3 — Adds 0.3 to the success score for men only. This encodes the historical bias directly in the label: men were rated more likely to succeed even when credentials were equal.
• hiring_data[hidden_gender == 0]['success'].mean() — Boolean indexing: selects only the rows where hidden_gender is 0 (men), then computes the mean success rate. Comparing this to the women’s rate reveals the bias baked into the labels.
• The key point: because hidden_gender never enters hiring_data, the model cannot use it directly — but it can learn from years_at_last_job, a visible feature that happens to correlate with gender in this historical data.

Running this prints an overall success rate of 72.2%, a 90.0% success rate for men, and a 55.8% success rate for women — even though the model never sees gender directly. The historical data carries the bias, and years_at_last_job correlates with it at -0.622.

Now let’s train a model and actually use SHAP to spot this — not just eyeball success rates by group:

X_hiring = hiring_data.drop('success', axis=1)
y_hiring = hiring_data['success']

X_train_h, X_test_h, y_train_h, y_test_h = train_test_split(
    X_hiring, y_hiring, test_size=0.2, random_state=42
)

model_hiring = RandomForestClassifier(n_estimators=100, random_state=42, max_depth=10)
model_hiring.fit(X_train_h, y_train_h)

print(f"Model accuracy: {model_hiring.score(X_test_h, y_test_h):.1%}")

# Now use SHAP
explainer_hiring = shap.TreeExplainer(model_hiring)
shap_values_hiring = explainer_hiring.shap_values(X_test_h)
shap_values_success = shap_values_hiring[:, :, 1]  # Positive class

# Find matched pairs: candidates with similar core credentials
# (years_experience, test_score, has_advanced_degree) but different gender.
# We deliberately leave years_at_last_job OUT of the "similar" check --
# it's the proxy feature we're trying to catch in the act.
n = len(X_test_h)
pairs = []
for i in range(n):
    for j in range(i + 1, n):
        if hidden_gender[X_test_h.index[i]] != hidden_gender[X_test_h.index[j]]:
            ci, cj = X_test_h.iloc[i], X_test_h.iloc[j]
            if (abs(ci['years_experience'] - cj['years_experience']) < 2
                    and abs(ci['test_score'] - cj['test_score']) < 2
                    and ci['has_advanced_degree'] == cj['has_advanced_degree']):
                pairs.append((i, j))

print(f"Matched pairs found: {len(pairs)}")

# Don't just print whichever pair matches first -- check the direction
# of the gap across every matched pair, and assert it.
gaps = []
for i, j in pairs:
    gi = hidden_gender[X_test_h.index[i]]
    pi = model_hiring.predict_proba(X_test_h.iloc[[i]])[0][1]
    pj = model_hiring.predict_proba(X_test_h.iloc[[j]])[0][1]
    men_pred, women_pred = (pi, pj) if gi == 0 else (pj, pi)
    gaps.append(men_pred - women_pred)
gaps = np.array(gaps)

print(f"Mean prediction gap (men - women) across matched pairs: {gaps.mean():.4f}")
print(f"Fraction of pairs favoring men: {(gaps > 0).mean():.1%}")
assert gaps.mean() > 0, "expected the matched-pair gap to favor men on average"

# Show one representative pair with its SHAP breakdown
best_idx = int(np.argmax(gaps))
i, j = pairs[best_idx]
gi, gj = hidden_gender[X_test_h.index[i]], hidden_gender[X_test_h.index[j]]

print(f"\nCandidate {i} (gender={gi}) vs Candidate {j} (gender={gj})")
print(X_test_h.iloc[[i, j]])
print(f"Prediction i: {model_hiring.predict_proba(X_test_h.iloc[[i]])[0][1]:.1%}")
print(f"Prediction j: {model_hiring.predict_proba(X_test_h.iloc[[j]])[0][1]:.1%}")
print(f"SHAP contributions, candidate {i}: {dict(zip(X_test_h.columns, np.round(shap_values_success[i], 4)))}")
print(f"SHAP contributions, candidate {j}: {dict(zip(X_test_h.columns, np.round(shap_values_success[j], 4)))}")

• shap_values_hiring[:, :, 1] — Selects contributions for the positive class (success = 1) from the last axis, consistent with the churn example earlier.
• for i in range(len(X_test_h)): for j in range(i+1, ...) — A nested loop that considers every unique pair of test candidates exactly once (since j always starts above i, we avoid duplicate pairs).
• Matching on years_experience, test_score, has_advanced_degree only — This is the important methodological choice: “similar credentials” means similar on the features that are supposed to determine hiring outcomes. years_at_last_job is deliberately excluded from the similarity check, because it’s the feature under suspicion — if we required it to match too, we’d match away the exact signal we’re trying to catch.
• gaps.mean() and the assert — Rather than printing whichever single pair happens to match first (which could go either direction by chance), this aggregates the prediction gap across every matched pair and checks that the average direction favors men — the systematic pattern, not a one-off.
• SHAP contributions printed per candidate — This is where SHAP actually gets used: for the pair with the largest gap, the two candidates’ years_at_last_job contributions point in opposite directions even though their other three features are nearly identical — that’s the proxy variable showing up directly in the SHAP breakdown, not just in the final prediction gap.

Running this on the test set finds 11 matched pairs, a mean prediction gap of +0.2007 favoring men, and men favored in 10 of the 11 pairs. The most extreme pair: a man and a woman with nearly identical years_experience (12.13 vs. 12.55), test_score (43.80 vs. 44.51), and has_advanced_degree (both 0) — but the man’s years_at_last_job is 7.19 while the woman’s is 0.00. Their predictions: 95.3% for the man, 51.9% for the woman. And the SHAP breakdown shows exactly why — years_at_last_job contributes +0.142 for the man and -0.147 for the woman, an almost perfect mirror image, while the three genuinely credential-based features contribute almost the same amount for both. That’s the proxy variable caught in the act: same core credentials, opposite prediction, and SHAP points straight at the one feature responsible.

The issue could be bias in the training data. Or it might be a proxy variable, like “years at last job,” which correlates with gender due to historical discrimination. Either way, SHAP helped you find it — and this time, it’s SHAP that found it, not just a comparison of group success rates.

Putting It All Together: A Practical Workflow

The Three Questions SHAP Answers

SHAP boils down to three questions you’ll keep returning to.

Question 1: Why did the model make this prediction? A force plot or waterfall plot handles a single prediction. This is the one you bring to your boss, your regulator, or your customer.

Question 2: Which features matter most overall? Use a summary bar plot. It shows which features drive the model’s decisions on average.

Question 3: Is there a hidden pattern or bias? Reach for dependence plots and matched-pairs comparisons like the one above, and watch for unexpected patterns. If two candidates with near-identical core credentials get very different predictions, something is wrong — and SHAP can tell you which feature is responsible.

A Complete Example

Here’s a final, complete example:

# Let's create a realistic scenario: predicting customer lifetime value

np.random.seed(42)
n_customers = 500

ltv_data = pd.DataFrame({
    'customer_age_years': np.random.uniform(1, 20, n_customers),
    'monthly_spend': np.random.uniform(10, 500, n_customers),
    'engagement_score': np.random.uniform(0, 100, n_customers),
    'num_purchases': np.random.randint(1, 50, n_customers),
    'days_since_last_purchase': np.random.randint(0, 365, n_customers),
})

# Target: high lifetime value, split at the *median* of a continuous score --
# not a collection of independent thresholds that happen to land wherever
# they land. This guarantees an actual ~50/50 split instead of an accidental
# 80/20 one.
ltv_score = (
    ltv_data['customer_age_years'] * 0.3 +
    ltv_data['monthly_spend'] * 0.4 +
    ltv_data['engagement_score'] * 0.2 +
    ltv_data['num_purchases'] * 0.1 +
    np.random.uniform(-5, 5, n_customers)
)
ltv_data['high_ltv'] = (ltv_score > ltv_score.median()).astype(int)

print(f"High LTV rate: {ltv_data['high_ltv'].mean():.1%}")

# Train model
X_ltv = ltv_data.drop('high_ltv', axis=1)
y_ltv = ltv_data['high_ltv']

X_train_l, X_test_l, y_train_l, y_test_l = train_test_split(
    X_ltv, y_ltv, test_size=0.2, random_state=42
)

model_ltv = RandomForestClassifier(n_estimators=100, random_state=42, max_depth=10)
model_ltv.fit(X_train_l, y_train_l)

print(f"Model accuracy: {model_ltv.score(X_test_l, y_test_l):.1%}")

# Create SHAP explainer
explainer_ltv = shap.TreeExplainer(model_ltv)
shap_values_ltv = explainer_ltv.shap_values(X_test_l)
shap_values_ltv_positive = shap_values_ltv[:, :, 1]

# 1. Explain a single prediction
print("\n=== EXPLAINING A SINGLE PREDICTION ===")
customer_idx = 0
customer_features = X_test_l.iloc[customer_idx]
prediction = model_ltv.predict_proba(X_test_l.iloc[[customer_idx]])[0][1]
customer_shap_vals = shap_values_ltv_positive[customer_idx]
base_val = explainer_ltv.expected_value[1]

print(f"Customer features:")
for feat, val in customer_features.items():
    print(f"  {feat}: {val:.2f}")

print(f"\nBase prediction (average): {base_val:.1%}")
print(f"\nFeature contributions:")
for feat, shap_val in zip(X_test_l.columns, customer_shap_vals):
    direction = "↑" if shap_val > 0 else "↓"
    print(f"  {feat}: {shap_val:+.4f} {direction}")

print(f"\nFinal prediction: {prediction:.1%}")
print(f"Verification: {base_val + customer_shap_vals.sum():.1%}")
assert abs(base_val + customer_shap_vals.sum() - prediction) < 1e-6

# 2. Feature importance across all predictions
print("\n=== FEATURE IMPORTANCE (AVERAGE SHAP VALUES) ===")
feature_importance = np.abs(shap_values_ltv_positive).mean(axis=0)
for feat, imp in sorted(zip(X_test_l.columns, feature_importance), key=lambda x: x[1], reverse=True):
    print(f"  {feat}: {imp:.4f}")

# 3. Check for unexpected patterns
print("\n=== CHECKING FOR PATTERNS ===")
print("\nCorrelation between monthly_spend and its SHAP values:")
corr = np.corrcoef(
    X_test_l['monthly_spend'].values,
    shap_values_ltv_positive[:, X_test_l.columns.get_loc('monthly_spend')]
)[0, 1]
print(f"  Correlation: {corr:.3f}")
print(f"  Interpretation: As monthly_spend increases, SHAP values {'increase' if corr > 0 else 'decrease'} (as expected)")

• ltv_score and (ltv_score > ltv_score.median()) — Builds a single continuous score first, then thresholds it at its own median. This is what actually guarantees a 50/50 class split, unlike summing several independent > threshold indicators (each with its own base rate), which can easily land at a lopsided split despite a comment claiming otherwise.
• customer_features.items() — Iterates over the (feature name, value) pairs of the pandas Series. Equivalent to zip(customer_features.index, customer_features.values) but more idiomatic.
• direction = "↑" if shap_val > 0 else "↓" — A one-line conditional expression (Python ternary operator) that picks an arrow character based on the sign of the SHAP value, making the printed output easier to scan.
• {shap_val:+.4f} — The + flag forces a sign character to always be printed (e.g., +0.0812 or -0.1234), making positive and negative contributions visually distinct.
• assert abs(base_val + customer_shap_vals.sum() - prediction) < 1e-6 — Same enforced-additivity pattern as the churn example.
• np.abs(shap_values_ltv_positive).mean(axis=0) — Takes the absolute value of every SHAP value (so positive and negative contributions don’t cancel), then averages across rows (axis=0 = across customers). The result is a 1-D array of length n_features: the mean absolute SHAP importance per feature.
• sorted(..., key=lambda x: x[1], reverse=True) — Sorts the list of (feature_name, importance) tuples by the second element (importance) in descending order.
• X_test_l.columns.get_loc('monthly_spend') — Returns the integer position of monthly_spend in the column list, which is needed to index into the NumPy SHAP array (which has no column names, only integer indices).
• np.corrcoef(...)[0, 1] — np.corrcoef returns a 2×2 correlation matrix; [0, 1] extracts the off-diagonal element, which is the Pearson correlation between the two input arrays.

This gives a 50.0% high-LTV rate (a genuine median split), a 96.0% test accuracy, and for the first test customer: age 13.22 years, monthly spend $232.11, engagement score 21.07, 28 purchases, and a purchase 67 days ago. The model gives them only a 16.0% chance of high lifetime value, even though a few of those numbers sound promising in isolation. SHAP explains why: monthly_spend is the single biggest drag on the prediction (-0.2927) — at $232, this customer’s spend sits close to the sample average rather than standing out, so it doesn’t push the prediction up the way a much higher spend would. Their engagement score, well below the sample average, adds a smaller further pull down (-0.0423). Their relatively recent purchase and above-average purchase count only nudge the prediction back up slightly (+0.0159 and effectively flat), not nearly enough to close a 34-point gap from the 50% baseline. And feature importance across the whole test set confirms monthly_spend dominates globally too, at 0.4662 — far ahead of the next feature.

That workflow gives you three things:

1. For a specific prediction: You can explain it to anyone — your boss, a regulator, a customer. “This customer has only a 16% chance of high lifetime value. Their monthly spend and engagement score are both working against them — neither stands out enough from the average to push the prediction up, and that’s most of the story here.”

2. For model debugging: You can see which features matter most. If a feature you expected to matter doesn’t, investigate. If a feature you didn’t expect to matter does, investigate.

3. For bias detection: You can spot when the model makes different predictions for candidates with similar core credentials, and use SHAP itself — not just the prediction gap — to see which feature is carrying the bias.

---

Key Takeaways

You now understand why your model made a specific prediction. Here’s what we covered:

• The black box problem: High accuracy isn’t enough. Regulators, stakeholders, and customers need explanations.
• The game theory intuition: SHAP values fairly split credit among features, just like the Shapley value in poker.
• The practical workflow: Use force plots for single predictions, summary plots for overall importance, and dependence plots to spot patterns — including patterns that turn out to be a cliff rather than a gradient.
• The debugging power: SHAP helps you find hidden biases, proxy variables, and unexpected patterns that accuracy alone would miss — as long as you actually call SHAP on the suspect comparison, not just eyeball group averages.

The next step? Start using SHAP on your own models. Pick a prediction that matters — a loan denial, a hiring decision, a medical diagnosis — and explain it. Show it to your stakeholders. You’ll be surprised how much trust it builds.

And when someone asks, “Why did your model predict that?” you’ll have a clear, data-driven answer.

Mr. Owens nods when he sees the force plot for Maria’s application — he can finally point to specific numbers and say why the model said no. But then he leans back and frowns: “This is thorough, but is there a simpler, faster tool I could use when a customer is sitting right in front of me?” That question is exactly what the next article tackles, pitting LIME against SHAP to find the right translator for every situation.

---

Check Your Understanding

Questions below move from simple recall up to open-ended design, roughly following Bloom’s Taxonomy.

Remember What are the three questions the article says SHAP can answer, and which plot type addresses each one?

Understand In your own words, explain why “just remove each feature one at a time” is an unreliable way to measure feature contribution, using the article’s point about ordering (remove income first vs. credit score first).

Apply Using the article’s worked example (base value 0.28, tenure -0.12, monthly_charges +0.08, num_services -0.05, age +0.02), calculate the final prediction if you also add a total_charges SHAP value of -0.03 — and verify it matches the article’s stated 0.18.

Analyze The hiring-bias example finds that gender isn’t a feature in the model, yet a matched-pairs SHAP comparison reveals two candidates with nearly identical core credentials getting very different predictions — and SHAP shows years_at_last_job pulling in opposite directions for the two of them. Walk through why this specific pattern (same core inputs, opposite SHAP contribution on one feature, different outputs) is the actual signal of bias here, rather than looking at a “gender” SHAP value directly (which doesn’t exist since gender was never a model input).

Evaluate The article distinguishes SHAP importance (“which features change the prediction most”) from traditional feature importance (“which features the model uses most”), calling them “related but not identical.” Critique relying on SHAP summary plots alone to decide which features to keep in a simplified model: what could a feature with high average |SHAP value| but low traditional importance actually indicate, and would dropping the “unimportant” traditional features be safe?

Create Design a SHAP-based debugging workflow for a new scenario: a resume-screening model where you suspect it’s implicitly penalizing candidates from a specific set of universities (a proxy for something like socioeconomic background) even though “university” isn’t a direct feature. Using the article’s hiring-bias pattern (matched-pairs comparison backed by SHAP, not just group averages), describe what you’d compare and what result would confirm your suspicion.

---

Related articles

• LIME vs SHAP: Choosing the Right Translator for You)
• Partial Dependence Plots and ICE Curves: See How Your Model Really Thinks)
• Why Did the Model Say No? Explaining Black-Box Decisions to the Business)

---

References & Further reading

• Lundberg, S. M., & Lee, S.-I. (2017). “A Unified Approach to Interpreting Model Predictions.” Advances in Neural Information Processing Systems 30 (NeurIPS 2017). https://arxiv.org/abs/1705.07874 — The foundational paper introducing SHAP values and the TreeSHAP algorithm.
• shap library documentation. https://shap.readthedocs.io — Full API reference, tutorials for tree models, deep learning models, and linear models, plus interactive plot examples.
• Kaggle: “Home Credit Default Risk” competition. https://www.kaggle.com/c/home-credit-default-risk — A real-world loan-default dataset where SHAP explanations are both practically meaningful and widely discussed in top-scoring notebooks; an excellent sandbox for applying the techniques in this article to Maria’s exact scenario.

---