Modern machine learning models often achieve strong predictive performance, but their internal logic can be difficult to interpret. In domains such as healthcare, finance, and policy, this lack of transparency creates problematic limitations. Decisions must be justified, validated, and aligned with domain knowledge. Interpretability becomes an important requirement.
SHAP, short for SHapley Additive exPlanations, provides a practical and theoretically grounded approach to understanding model predictions. The per-feature contributions it assigns, called SHAP values, support both detailed explanations of individual predictions and broader insights into overall model behavior.
SHAP values in practice
Many commonly used models such as gradient boosting, random forests, and neural networks are considered black boxes. While they capture complex relationships, they do not offer straightforward explanations for their outputs.
SHAP bridges this gap by translating model predictions into feature-level contributions. Each prediction is decomposed into additive components, making it possible to answer a simple but critical question. Which features drove this prediction, and in what direction?
That matters most in settings where:
• Model outputs influence high-stakes decisions
• Domain experts need to validate model behavior
• Regulatory or audit requirements demand transparency
Shapley values
SHAP is based on Shapley values from cooperative game theory. In that setting, each player contributes to a collective outcome, and the goal is to fairly assign credit.
In machine learning, the players are features and the outcome is the model prediction. SHAP assigns each feature a contribution value that reflects how much it shifts the prediction away from a base value, typically the average model output.
This leads to an additive explanation:
• Baseline prediction
• Plus contributions from each feature
• Equals the final prediction
The additive structure is more than a formality. It guarantees that the base value plus the feature contributions add up exactly to the model output for the prediction being explained.
The useful question is how each feature earns its share of that total. In the Shapley framework a feature is valued by its marginal contribution, the amount the model output moves when the feature joins a coalition of the features already present, that is, any subset of the others. Because that amount can depend on which features are already there, the contribution is averaged across every possible coalition, weighted so that each coalition size carries equal weight. No single ordering of the features decides the result. Figure 1 works through this averaging one feature at a time.
Figure 1
How a feature earns its SHAP value
A feature’s SHAP value is the average of its marginal contributions across the coalitions of the other features, weighted so each coalition size counts equally. Choose a feature and step through the coalitions.
| Coalition already present | Marginal | Weight |
|---|---|---|
| { } none | +0.59 | 1/4 |
| Waist | +0.93 | 1/12 |
| HDL | +0.59 | 1/12 |
| Family history | +0.59 | 1/12 |
| Waist, HDL | +0.93 | 1/12 |
| Waist, Family history | +0.93 | 1/12 |
| HDL, Family history | +0.59 | 1/12 |
| Waist, HDL, Family history | +0.93 | 1/4 |
| Weighted sum = SHAP value | +0.76 | 1 |
Numbers come from a transparent illustrative model of type 2 diabetes risk in log-odds (a fasting glucose by waist circumference interaction plus additive terms), with SHAP computed exactly over all coalitions against a fixed background sample. This is a separate example from the cardiovascular model worked below. It illustrates the mechanism only. All values are rounded to two decimals.
Figure 1. A feature’s SHAP value is the weighted average of its marginal contributions across every coalition of the other features. Select a feature and step through the coalitions, or press Play to watch the weighted sum build toward the final value. Fasting glucose contributes more when waist circumference is already present, because the two interact, while an additive feature such as high-density lipoprotein (HDL) cholesterol contributes the same amount to every coalition. Weighting every coalition by size is what defines the SHAP value.
Computational considerations
Exact computation of Shapley values requires evaluating all possible feature subsets, which is computationally infeasible for most real-world models.
In practice, more efficient adaptations are used:
• TreeSHAP for tree-based models
• KernelSHAP for model-agnostic settings
• Sampling-based approaches for large feature spaces
TreeSHAP is particularly important in applied work since it provides exact SHAP values for tree ensembles with polynomial complexity. That aspect makes it feasible to use SHAP routinely in production workflows involving gradient boosting models. A practical best practice is to use model-specific SHAP implementations whenever available. They are faster, and for tree ensembles TreeSHAP returns exact values where the generic methods approximate.
Local interpretability: explaining individual predictions
SHAP is well suited to explaining single predictions. Several visualization tools support this.
Waterfall plots
Waterfall plots show how each feature contributes to a prediction step by step. Starting from the base value, each feature pushes the prediction higher or lower.
This is especially useful in case-level analysis. For example, a clinical model may show that:
• Age and blood pressure increase predicted risk
• Favorable lab values reduce it
The result is a clear narrative that aligns with domain reasoning.
One detail matters before reading a waterfall, namely, the scale on which the contributions add up. For a classifier such as the gradient boosted model used later, SHAP explains the output in log-odds, so the base value and every contribution are expressed in log-odds rather than in probability. Figure 2 builds a single prediction on that scale and reports the probability separately.
Figure 2
Building one prediction from the base value
Starting from the base value, each feature adds its SHAP value in log-odds until the sum reaches the model output for that patient. Step through the contributions or switch patients.
Patient A has elevated fasting glucose and a large waist, but high HDL and no family history of diabetes. Illustrative values from the same transparent model as Figure 1.
Figure 2. Beginning at the base value E[f(X)], the average model output across the background patients, each feature adds its SHAP value in log-odds until the running total reaches the model output f(x) for that patient. Switch between the three patients or step through the contributions. The probability beside the chart is the logistic transform of the log-odds total, read off the two endpoints. The contributions sum to the log-odds, and the curve does the rest. Keeping the two scales separate avoids a common misreading.
Force plots
Force plots are among the most distinctive SHAP visualizations. They represent the additive nature of SHAP values through a dynamic balance of forces pushing and pulling the model’s output from the baseline toward the final prediction. Each force plot starts from the expected value, marked as the base value on the output axis, and the opposing forces meet at the model output f(x). Arrows extending to the right represent features that increase the model’s output, while those to the left reduce it. The width of each arrow corresponds to the magnitude of the SHAP value. The final prediction is the result of all these forces acting together.
Force plots are ideal for explaining individual predictions in settings such as clinical consults or patient education. For instance, in a model assessing the risk of postoperative complications, a clinician can use the force plot to show that while the patient’s age increased the predicted risk, favorable factors like normal renal function and low surgical complexity brought it back down. This clear visualization bridges the gap between statistical reasoning and clinical storytelling.
Beeswarm plots
The beeswarm plot is a dense summary of a SHAP analysis. It combines feature importance with the direction and spread of each feature’s contribution.
Each point represents a SHAP value for a single observation. The horizontal position shows the SHAP value, and color often encodes the original feature value.
That layout supports several readings at once:
• Which features are most important overall
• Whether higher values increase or decrease predictions
• Whether contributions are consistent or vary from patient to patient
In practice, beeswarm plots often reveal non-linear patterns or interactions that are not obvious from model coefficients alone.
Figure 3 connects one of these points to the whole population, and to the bar chart of mean absolute SHAP values that usually sits beside a beeswarm.
Figure 3
From single explanations to a global summary
Each patient has one SHAP value per feature (left). Averaging the absolute size of those values across patients gives the global importance bars (right). Select a feature or filter by value to see how local spread becomes one number.
A population of 60 illustrative patients from the same model as Figures 1 and 2, with SHAP computed exactly per patient. Mean |SHAP| aggregates magnitude. Direction sits in the per-patient values on the left, and causal effect comes from study design.
Figure 3. Each dot is one patient’s SHAP value for a feature. Averaging the absolute values across patients produces the global importance bars on the right. Select a feature, or filter the dots by feature value, to see how the local spread becomes a single number. The bars rank features by average magnitude, which discards direction. Waist circumference and HDL cholesterol have large average magnitude. They raise the prediction for some patients and lower it for others, and the bars cannot show that. Averaging the signed values does not recover it either, since the opposing contributions cancel. The same two features average −0.07 and −0.10 in signed terms, against magnitudes above 0.4. Reading the beeswarm alongside the bars keeps magnitude and direction distinct, and neither is a statement about cause.
Heatmaps
SHAP heatmaps provide a matrix-like overview of SHAP values across many instances and features. The rows correspond to features, ordered by mean absolute SHAP, and the columns represent individual patients. Each cell’s color reflects the SHAP value, often with red tones indicating positive contributions and blue tones indicating negative ones. This can uncover subgroups in the data where the model behaves differently. In healthcare applications, this often aligns with clinically meaningful subpopulations.
Example: Explaining a gradient boosted trees model
ExplainingPredictions
How SHAP turns an opaque gradient boosting model for cardiovascular risk into clear, additive explanations, from a global feature ranking down to a single patient’s prediction.
Explaining Predictions. How SHAP turns an opaque gradient boosting model for cardiovascular risk into clear, additive explanations, from a global feature ranking down to a single patient's prediction. Key topics covered: Cardiovascular dataset, Gradient boosting model, Beeswarm summary, Shapley values, Waterfall plot, Force plot, SHAP heatmap, Top features by mean |SHAP|.
This notebook provides a complete, reproducible workflow for training a gradient-boosted classifier on cardiovascular data and interpreting its predictions with SHAP (SHapley Additive exPlanations).
SHAP is grounded in cooperative game theory. Each feature is treated as a player in a coalition and assigned a SHAP value, which is its fair marginal contribution to a given prediction relative to the model’s baseline (expected) output. SHAP values satisfy three key properties that make them particularly suitable for clinical interpretation:
- Local accuracy: SHAP values sum to the model output:
base_value + Σφᵢ = f(x). - Consistency: if a model changes so that a feature contributes more no matter which other features are present, its SHAP value will not decrease.
- Missingness: features absent from a prediction receive a SHAP value of zero.
Workflow overview
| Section | Content |
|---|---|
| 1–2 | Environment setup and data loading |
| 3–5 | Feature selection, train-test split, and missing value imputation |
| 6 | Gradient Boosting Classifier training |
| 7 | TreeSHAP initialization and global SHAP value computation |
| 8–9 | Global summary: beeswarm and layered violin plots |
| 10 | Decision plot: cumulative feature effects per patient |
| 11 | Heatmap: SHAP values across the population |
| 12–13 | Local explanations: waterfall and force plots |
| 14 | SHAP vs impurity (MDI) feature importance comparison |
Dataset: Cardiovascular Disease dataset (Kaggle – Aidan, https://www.kaggle.com/datasets/colewelkins/cardiovascular-disease).
Contains anonymized records for approximately 70,000 patients. Features include blood pressure, cholesterol, body mass index (BMI), age, and physical activity.
Target: binary indicator of cardiovascular disease (1 = present, 0 = absent).
1. Setup: import dependencies
| Library | Purpose |
|---|---|
shap |
SHAP value computation and visualization (pip install shap) |
numpy / pandas |
Numerical arrays and tabular data manipulation |
matplotlib |
Plot rendering backend |
sklearn |
Model training and imputation utilities |
string |
Letter labels for anonymizing individual patients in plots |
Note:
enable_iterative_imputermust be imported beforeIterativeImputerbecause the class is still marked experimental in scikit-learn. The import registers it. The# noqacomment suppresses the “imported but unused” linter warning.
import string
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import shap
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingClassifier, RandomForestRegressor
from sklearn.experimental import enable_iterative_imputer # noqa: F401 - activates IterativeImputer
from sklearn.impute import IterativeImputer, SimpleImputer2. Load data and clean feature names
The raw CSV uses abbreviated column names (e.g., ap_hi for systolic arterial pressure, ap_lo for diastolic). We rename every column to a readable clinical label before anything else happens. Those names then carry through every downstream plot, which cuts the risk of misreading one.
Two features are immediately cast to pandas category dtype:
- Cholesterol, an ordinal three-level scale: 1 = normal, 2 = above normal, 3 = well above normal.
- Physical Activity, a binary flag: 1 = physically active, 0 = inactive.
units_map stores display labels that append measurement units (e.g., "Systolic BP (mmHg)"). These richer labels are reserved for individual-patient plots where the exact scale aids clinical interpretation.
df = pd.read_csv("cardio_data_processed.csv")
# Rename abbreviated raw column names to human-readable clinical labels
df.rename(
columns={
"ap_hi": "Systolic BP",
"ap_lo": "Diastolic BP",
"cholesterol": "Cholesterol",
"gluc": "Glucose",
"active": "Physical Activity",
"cardio": "Cardiovascular Disease",
"age_years": "Age",
"bmi": "Body Mass Index",
},
inplace=True,
)
# Encode ordinal / binary columns as categorical to signal their discrete nature
df["Physical Activity"] = df["Physical Activity"].astype("category")
df["Cholesterol"] = df["Cholesterol"].astype("category")
# Full labels with units, used in per-patient plots where scale context matters
units_map = {
"Age": "Age (years)",
"Physical Activity": "Physical Activity (0 | 1)",
"Systolic BP": "Systolic BP (mmHg)",
"Diastolic BP": "Diastolic BP (mmHg)",
"Cholesterol": "Cholesterol (1 | 2 | 3)",
"Body Mass Index": "Body Mass Index (kg/m\u00b2)",
}3. Feature selection
We predict cardiovascular disease from six clinically established risk factors:
| Feature | Type | Clinical Relevance |
|---|---|---|
| Age | Continuous | Cardiovascular risk increases substantially with age |
| Physical Activity | Binary | Regular exercise is a major protective factor |
| Systolic BP | Continuous | Primary hypertension marker; key cardiovascular disease (CVD) risk driver |
| Diastolic BP | Continuous | Complements systolic in assessing overall BP burden |
| Cholesterol | Ordinal (1-3) | Elevated cholesterol is a well-established CVD risk factor |
| Body Mass Index | Continuous | Obesity correlates with hypertension and metabolic syndrome |
Two parallel name lists are maintained:
features_short: used as model input column names and for concise axis labels in summary plots.features_full: annotated with measurement units; used in per-patient plots (waterfall, force) where the exact scale aids clinical reading.
# Short names: used as model column names and in summary plot axes
features_short = [
"Age",
"Physical Activity",
"Systolic BP",
"Diastolic BP",
"Cholesterol",
"Body Mass Index",
]
# Full names with units: used in per-patient plots for clinical readability
features_full = [units_map[c] for c in features_short]
X = df[features_short].copy()
y = df["Cardiovascular Disease"]4. Train–test split
We split the data 80% training / 20% test.
Key design choices:
stratify=y: preserves the class balance (CVD : non-CVD ratio) in both splits. Without stratification, random chance could produce a test set with a different prevalence than the training set, leading to misleading evaluation metrics.random_state=42: fixes the split for reproducibility.- BMI is rounded to two decimal places in the test set, reflecting realistic clinical measurement precision.
Data leakage principle: all preprocessing statistics (imputer parameters, scalers, encoders) must be estimated on the training set only and then applied to the test set. Fitting on the combined dataset would allow test-set information to influence training, inflating apparent model performance. We enforce this throughout.
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.20, random_state=42, stratify=y
)
# Make an explicit independent copy before modifying to avoid SettingWithCopyWarning
X_test = X_test.copy()
X_test["Body Mass Index"] = X_test["Body Mass Index"].round(2)5. Missing value imputation
Real-world clinical datasets routinely contain missing values. We apply different imputation strategies depending on the measurement scale of each feature.
Continuous features: Iterative Imputation (MICE)
IterativeImputer implements MICE (Multivariate Imputation by Chained Equations). It models each feature with missing values as a function of all other features, iterating until the imputed values converge. The scikit-learn class departs from the procedure the acronym names in one respect that matters here. It returns a single completed dataset rather than pooling several imputations, so it does not carry imputation uncertainty forward. A RandomForestRegressor is used as the base estimator, which handles nonlinear feature interactions without distributional assumptions.
initial_strategy='median': seeds the first iteration with the column median, which is robust to outliers.max_iter=20: maximum number of imputation rounds.n_estimators=500: sufficient trees for stable per-iteration estimates.
Categorical features: Mode imputation
SimpleImputer(strategy='most_frequent') replaces missing categories with the most frequent observed value. This is a conservative, non-parametric approach appropriate for ordinal and binary features.
Preventing data leakage
Both imputers are fit on the training set only (fit_transform), then applied to the test set with transform. The category dtype is restored after imputation because scikit-learn imputers return plain NumPy arrays and do not preserve pandas metadata.
print("Performing mixed-type imputation...")
# Separate features by measurement scale for targeted imputation strategies
num_cols = [c for c in features_short if c not in ("Physical Activity", "Cholesterol")]
cat_cols = ["Physical Activity", "Cholesterol"]
print(f" Continuous columns (MICE): {num_cols}")
print(f" Categorical columns (mode): {cat_cols}")
X_train_imp = X_train.copy()
X_test_imp = X_test.copy()
# --- Continuous features: iterative imputation (MICE) ---
if num_cols:
num_imputer = IterativeImputer(
estimator=RandomForestRegressor(n_estimators=500, random_state=42),
max_iter=20,
random_state=42,
initial_strategy="median", # robust first-pass estimate
)
# Fit on training data only, then transform both sets
X_train_imp[num_cols] = num_imputer.fit_transform(X_train_imp[num_cols])
X_test_imp[num_cols] = num_imputer.transform(X_test_imp[num_cols])
# --- Categorical features: mode imputation ---
if cat_cols:
cat_imputer = SimpleImputer(strategy="most_frequent")
# Fit on training data only, then transform both sets
X_train_imp[cat_cols] = cat_imputer.fit_transform(X_train_imp[cat_cols])
X_test_imp[cat_cols] = cat_imputer.transform(X_test_imp[cat_cols])
X_train = X_train_imp.copy()
X_test = X_test_imp.copy()
# Restore category dtype (scikit-learn imputers return plain numpy arrays)
for col in cat_cols:
X_train[col] = X_train[col].astype("category")
X_test[col] = X_test[col].astype("category")
# Ensure continuous columns are stored as float64
if num_cols:
X_train[num_cols] = X_train[num_cols].astype(np.float64)
X_test[num_cols] = X_test[num_cols].astype(np.float64)
print("Imputation complete.")
# Validate: report any remaining missing values
for label, frame in (("Training", X_train), ("Test", X_test)):
n_missing = frame.isnull().sum().sum()
if n_missing > 0:
print(f" Warning: {n_missing} missing values remain in {label} set")
print(frame.isnull().sum()[frame.isnull().sum() > 0])
else:
print(f" {label} set: no missing values.")Performing mixed-type imputation...
Continuous columns (MICE): ['Age', 'Systolic BP', 'Diastolic BP', 'Body Mass Index']
Categorical columns (mode): ['Physical Activity', 'Cholesterol']
Imputation complete.
Training set: no missing values.
Test set: no missing values.If no warnings are printed above, imputation completed without residual missing values. Note that category dtypes are explicitly restored post-imputation. Scikit-learn’s imputers operate on NumPy arrays internally and return float arrays, discarding pandas metadata such as categorical encoding.
6. Training the gradient boosting classifier
GradientBoostingClassifier builds an ensemble of decision trees sequentially. Each successive tree is trained to correct the residual errors of the current ensemble, a strategy known as gradient boosting (Friedman, 2001).
Hyperparameter rationale:
| Parameter | Value | Rationale |
|---|---|---|
n_estimators |
1000 | Large ensemble for high expressive power |
learning_rate |
0.05 | Low shrinkage: each tree contributes modestly, reducing overfitting |
subsample |
0.8 | Stochastic GBM: each tree sees 80% of training samples, adding variance reduction analogous to bagging |
random_state |
42 | Reproducibility |
The pairing of many trees with a low learning rate is a well-established regularization strategy. The model learns slowly but thoroughly. subsample < 1.0 introduces randomness that further guards against overfitting and typically improves generalization.
gb_clf = GradientBoostingClassifier(
n_estimators=1000,
learning_rate=0.05,
subsample=0.8,
random_state=42,
)
gb_clf.fit(X_train, y_train)GradientBoostingClassifier(learning_rate=0.05, n_estimators=1000,
random_state=42, subsample=0.8)In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
GradientBoostingClassifier(learning_rate=0.05, n_estimators=1000,
random_state=42, subsample=0.8)7. Computing SHAP values with TreeSHAP
Why TreeSHAP?
shap.TreeExplainer uses the exact TreeSHAP algorithm (Lundberg et al., 2020), which exploits the tree structure to compute SHAP values in polynomial time, orders of magnitude faster than evaluating all 2ⁿ feature subsets. That speed makes exact SHAP computation practical for large ensembles.
Output scale: raw log-odds
model_output="raw" returns SHAP values on the model’s internal log-odds (decision function) scale:
f(x) = base_value + φ₁ + φ₂ + ... + φₙ
where base_value is the average log-odds prediction across the training set and φᵢ is the SHAP value for feature i. To convert to a probability: p = 1 / (1 + exp(-f(x))).
Keeping values on the log-odds scale is preferred for gradient boosting models because it preserves the additive structure of the model’s internal representation.
Sampling for efficiency
We draw up to 10,000 test instances for global SHAP summaries. This subset is large enough to represent the full distribution while keeping computation time manageable.
API version note: older SHAP versions return
shap_values()as alistwith one array per class for binary classifiers. Modern versions return a single array for the positive class. Theisinstance(raw_shap, list)check handles both behaviors transparently.
# Build the SHAP explainer for the trained gradient boosting ensemble.
# model_output="raw" returns SHAP values on the log-odds scale.
explainer = shap.TreeExplainer(gb_clf, model_output="raw")
# Sample up to 10,000 test instances for efficient global visualisation
sample_n = min(10000, len(X_test))
X_sample = X_test.sample(sample_n, random_state=42)
# Compute SHAP values for the sampled test set
raw_shap = explainer.shap_values(X_sample)
# Handle both old SHAP API (list per class) and new API (single array)
if isinstance(raw_shap, list):
shap_vals = raw_shap[1] # index 1 = positive class (CVD = 1)
base_val = explainer.expected_value[1]
else:
shap_vals = raw_shap
base_val = explainer.expected_value
# Guarantee base_val is a plain Python float regardless of SHAP version.
# Some versions return a 0-d numpy array, which cannot be used in f-string
# format specs (:.4f) or passed to np.full() as a scalar reliably.
base_val = float(np.squeeze(base_val))
# Wrap into a SHAP Explanation object, required by modern plot functions
expl_global = shap.Explanation(
values=shap_vals,
base_values=np.full(sample_n, base_val),
data=X_sample.values,
feature_names=features_short,
)
print(f"Base value (expected log-odds): {base_val:.4f}")
print(f"Equivalent baseline probability: {1 / (1 + np.exp(-base_val)):.4f}")
print(f"SHAP value array shape: {shap_vals.shape} (samples x features)")Base value (expected log-odds): -0.0219
Equivalent baseline probability: 0.4945
SHAP value array shape: (10000, 6) (samples x features)8. Beeswarm plot: global feature importance and effect direction
The beeswarm plot is the primary SHAP summary visualization. Each point represents one patient–feature pair:
- Horizontal position (x-axis): SHAP value: how much this feature moved the prediction (in log-odds) away from the baseline. Positive = increased predicted CVD risk; negative = decreased predicted risk.
- Color: raw feature value (red = high, blue = low).
- Vertical ordering: features sorted top-to-bottom by mean |SHAP|, i.e., overall importance.
In a single plot, the beeswarm reveals three dimensions of information simultaneously:
- Feature importance: the vertical rank.
- Effect direction: whether high or low feature values push the prediction up or down.
- Effect heterogeneity: the spread of the point cloud reflects how consistently a feature acts across patients.
shap.plots.beeswarm(expl_global, show=False)
plt.title("SHAP Beeswarm: Cardiovascular Disease Model")
plt.tight_layout()
plt.savefig("shap_beeswarm.png", dpi=200)
plt.show()SHAP beeswarm across 10,000 test patients
Interpreting the beeswarm plot:
- Systolic BP (top row) is the dominant predictor. High values (red) produce large positive SHAP values, strongly increasing predicted CVD risk. The effect is approximately monotone. Progressively higher blood pressure maps to progressively higher log-odds.
- Age shows a consistent positive gradient. Older patients (red) carry higher SHAP values, reflecting the well-established age–CVD relationship.
- Cholesterol follows the same direction. A higher ordinal level (red = level 3) is associated with increased predicted risk.
- Body Mass Index contributes positively at high values, consistent with clinical evidence linking obesity to cardiovascular outcomes.
- Diastolic BP has a smaller but visible positive contribution. The beeswarm cannot show whether that contribution is independent of systolic BP. The two are strongly correlated, and how credit is split between correlated features depends on how the explainer treats the features left out of a coalition.
- Physical Activity is the only feature with a predominantly negative SHAP contribution for high values. Active patients (red, value = 1) receive lower predicted risk. That describes how the model uses the feature, and stops short of any claim about exercise and the disease.
The wide spread of points along the x-axis shows that a feature’s contribution varies from patient to patient. A linear model produces spread too, whenever the feature itself varies, so width alone is not evidence of nonlinearity.
9. Layered violin plot: SHAP value distribution by feature
The layered violin plot conveys the same information as the beeswarm but replaces individual points with kernel density estimates. For large datasets (here 10,000 instances), this produces a less cluttered view and makes the density of SHAP values at each position more visible.
- Width of each layer: density of SHAP values at that point on the x-axis.
- Color encodes feature value (red = high, blue = low), as in the beeswarm.
- Vertical ordering is identical to the beeswarm, with features ranked by mean |SHAP|.
Use the violin when you want to emphasize the shape of the SHAP distribution (e.g., bimodal effects for binary features), and the beeswarm when you want to reveal individual outliers.
shap.plots.violin(expl_global, plot_type="layered_violin", show=False)
plt.title("SHAP Layered Violin Plot: Cardiovascular Disease Model")
plt.tight_layout()
plt.savefig("shap_layered_violin.png", dpi=200)
plt.show()Layered violin of SHAP distributions
Interpreting the violin plot:
The violin plot shows the same SHAP values as the beeswarm in a different form, restating those findings without adding to them. For Systolic BP, the wide right-shifted (positive) density mass for high values and the left-shifted density for low values show a clear, wide separation, consistent with its top rank by mean absolute SHAP. The Physical Activity violin shows a characteristically bimodal shape. The two density layers (active vs inactive) are separated and point in opposite directions, reflecting the binary nature of this feature and the consistent, directionally opposing contributions the model assigns to its two levels.
10. Decision plot: cumulative feature effects across patients
The decision plot visualizes how each feature cumulatively shifts a prediction from the model baseline to its final output, with one line per patient. Reading the plot from the bottom up, each row adds one feature’s SHAP value to the running total.
- x-axis: predicted value. We use
link="logit"to convert from log-odds to probability (0–1), which is more clinically intuitive. The decomposition is additive in log-odds only. The link warps the axis, so equal SHAP values are not equal horizontal distances and the contributions do not sum on the plotted scale. - y-axis: features, listed in the order their SHAP values are accumulated (bottom to top by default).
- Each line: one patient. Lines that end at the same point on the top row share similar final predictions. Diverging paths reveal that different features drive the prediction for different patients.
- Line color: reflects the final predicted probability (cool = low risk, warm = high risk).
To keep the plot readable, 26 patients are drawn at random. The letter map assigns each row a letter from its position in the full dataframe modulo 26, so a random draw of 26 repeats some letters and misses others rather than covering A to Z once each, and decision_plot does not render the letters on the figure.
# Map each original DataFrame row index to an anonymised letter label (A–Z, cycling)
LETTERS = string.ascii_uppercase
def num_to_letter(n: int) -> str:
"""Convert a 0-based integer to a letter: 0->A, 1->B, ..., 25->Z, 26->A, ..."""
return LETTERS[n % 26]
letter_map = {idx: num_to_letter(i) for i, idx in enumerate(df.index)}
# Select 26 random patients for a readable decision plot
np.random.seed(42)
subset_idx = np.random.choice(
X_sample.index,
size=min(26, sample_n),
replace=False,
)
X_decision = X_sample.loc[subset_idx].copy()
X_decision.index = [letter_map[i] for i in X_decision.index]
# Compute SHAP values for the decision plot subset
dec_shap = explainer.shap_values(X_decision)
dec_shap = dec_shap[1] if isinstance(dec_shap, list) else dec_shap
plt.figure(figsize=(11, 7))
shap.decision_plot(
base_val,
dec_shap,
X_decision,
feature_names=features_short,
link="logit", # convert log-odds to probability for a clinically meaningful x-axis
highlight=None,
show=False,
)
plt.title("SHAP Decision Plot: 26 Randomly Selected Patients")
plt.tight_layout()
plt.savefig("shap_decision_plot.png", dpi=200)
plt.show()Decision plot for 26 random patients
Interpreting the decision plot:
Lines clustered near the right of the plot represent high-risk patients. Those near the left are low-risk. What decision plots show is where divergence occurs. Lines begin to separate most strongly at the Systolic BP band, the same feature the beeswarm ranks first. Both readings come from the same SHAP values, so the plot shows that ranking again and adds no independent support. For some patients, lines cross the midline, where an initially below-average predicted risk is pushed up by the contributions from blood pressure or age, making each feature’s share of the total transparent and quantified. The bottom-to-top order is a plotting convention, and the model follows no such sequence. Patients whose lines stay near-vertical across multiple features (no large SHAP contributions) have a predicted risk close to the population baseline.
11. SHAP heatmap: population-level feature attribution
The SHAP heatmap provides a matrix view of SHAP values across many patients simultaneously:
- Rows: features, ordered by mean absolute SHAP with the most important at the top.
- Columns: individual patients (up to 5,000), ordered by hierarchical clustering of their SHAP profiles.
- Cell color: SHAP value: red = positive contribution to the predicted risk, blue = negative.
- Side panels: an f(x) trace above the matrix tracking the prediction for each patient, and a mean absolute SHAP bar panel at the right. The function draws no dendrogram.
This visualization is good for spotting patient clusters with distinct explanatory structures. For example, a cluster whose predicted risk is attributed mostly to blood pressure versus one attributed mostly to age and metabolic factors. Such patterns may guide population-level clinical strategies or suggest model validation across subgroups.
# Select up to 5,000 patients for the heatmap; more samples improve the cluster resolution
np.random.seed(42)
heat_idx = np.random.choice(
X_sample.index,
size=min(5000, sample_n),
replace=False,
)
X_heat = X_sample.loc[heat_idx].copy()
# Compute SHAP values for the heatmap subset
heat_sh = explainer.shap_values(X_heat)
heat_sh = heat_sh[1] if isinstance(heat_sh, list) else heat_sh
# Note: shap.plots.heatmap creates its own figure internally.
# Do NOT call plt.figure() beforehand, or an empty blank figure will be rendered.
shap.plots.heatmap(
shap.Explanation(
values=heat_sh,
base_values=np.full(len(X_heat), base_val),
data=X_heat.values,
feature_names=features_short,
),
show=False,
)
plt.title("SHAP Heatmap: 5,000 Test Patients")
plt.tight_layout()
plt.savefig("shap_heatmap.png", dpi=200)
plt.show()SHAP heatmap of 5,000 test patients
Interpreting the heatmap:
The hierarchical clustering typically reveals two broad patient cohorts, those with predominantly red cells along the Systolic BP row (high positive SHAP → high predicted risk) and those with predominantly blue cells (low or negative BP contribution → low predicted risk). Within the high-risk cluster, secondary separation often reflects age and cholesterol contributions. Bands of consistent color across many patients along a single row show that feature’s globally consistent direction of contribution, which is also visible in the beeswarm. Diagonal or patchy patterns along a feature row suggest heterogeneous, nonlinear, or interaction-dependent contributions.
12. Waterfall plots: dissecting individual predictions
The waterfall plot provides a detailed, per-patient decomposition of a single prediction:
- The plot starts at
E[f(X)]: the model’s baseline (expected) log-odds output. - Each bar adds or subtracts that feature’s SHAP value, building step-by-step toward the final prediction
f(x)shown at the top. - Red bars (rightward): features that push the log-odds prediction up.
- Blue bars (leftward): features that push the log-odds prediction down.
- Bar width reflects SHAP magnitude. The feature value for this specific patient is shown alongside each label.
We generate waterfall plots for four randomly selected patients. The title shows the model’s predicted class and probability, giving immediate clinical context for the decomposition.
# Fix the random seed so the same four patients are selected each run
np.random.seed(42)
plot_width = 9 # inches
plot_height = 5 # inches
for idx in np.random.choice(X_sample.index, size=4, replace=False):
letter = letter_map[idx]
x_i = X_sample.loc[[idx]]
# Compute SHAP values for this single patient
sv_i = explainer.shap_values(x_i)
sv_i = sv_i[1] if isinstance(sv_i, list) else sv_i
# Use features_full (with units) so axis labels are clinically informative
shap.plots.waterfall(
shap.Explanation(
values=sv_i[0],
base_values=base_val,
data=x_i.values[0],
feature_names=features_full,
),
show=False,
)
proba = gb_clf.predict_proba(x_i)[0, 1]
pred = gb_clf.predict(x_i)[0]
plt.gcf().set_size_inches(plot_width, plot_height)
plt.title(f"Patient {letter} \u2014 {'CVD' if pred else 'Non-CVD'} (P = {proba:.3f})")
plt.tight_layout()
plt.savefig(f"waterfall_plot_patient_{letter}.png", dpi=200)
plt.show()Waterfall decomposition for patient M
Waterfall decomposition for patient L
Waterfall decomposition for patient T
Waterfall decomposition for patient O
Interpreting the waterfall plots:
Each patient’s waterfall tells a unique clinical story. For a high-risk patient (CVD, high probability), expect large red bars for Systolic BP and/or Age driving the prediction well above the baseline. For a low-risk patient (Non-CVD, low probability), wide blue bars, most commonly reflecting normal blood pressure, younger age, and physical activity, pull the prediction below the baseline. Comparing two patients with similar final probabilities but different bar patterns reveals that the model’s internal reasoning differs. The predicted risk rests mostly on blood pressure for one patient and mostly on age and cholesterol for another. This patient-level narrative is the primary value of waterfall plots in clinical decision support.
13. Force plot: visualizing feature forces for a single patient
The force plot is an alternative single-instance visualization that emphasizes the balance of opposing forces acting on a prediction:
- The prediction starts at
base value: the model’s population-average output. - Red arrows (rightward): features pushing toward higher CVD risk.
- Blue arrows (leftward): features pushing toward lower risk.
- Arrow width encodes SHAP magnitude.
- The final prediction is where all forces come to rest.
Compared to the waterfall, the force plot sacrifices the step-by-step accumulation view in favor of a more immediate visual impression of the dominant factors. It is particularly effective for patient communication and clinical case presentations. The tug-of-war metaphor is intuitive for non-technical audiences.
Selecting a patient: adjust
PATIENT_POSITION(0-indexed) to examine a different test-set patient.
# Select a patient by their 0-based position in the test set.
# Change PATIENT_POSITION to examine a different patient.
PATIENT_POSITION = 22
USER_PATIENT_IDX = X_test.index[PATIENT_POSITION]
x_sel = X_test.loc[[USER_PATIENT_IDX]]
# Compute SHAP values for the selected patient
sv_sel = explainer.shap_values(x_sel)
sv_sel = sv_sel[1] if isinstance(sv_sel, list) else sv_sel
letter_sel = letter_map[USER_PATIENT_IDX]
shap.force_plot(
base_val,
sv_sel[0],
x_sel.values,
feature_names=features_full,
matplotlib=True,
figsize=(24, 4),
show=False,
)
proba_sel = gb_clf.predict_proba(x_sel)[0, 1]
pred_sel = gb_clf.predict(x_sel)[0]
plt.title(
f"Patient {letter_sel} \u2014 {'CVD' if pred_sel else 'Non-CVD'} (P = {proba_sel:.3f})"
)
plt.tight_layout()
plt.savefig(f"force_plot_patient_{letter_sel}.png", dpi=200)
plt.show()Force plot for patient K
Interpreting the force plot:
The width and direction of each arrow immediately communicate which features moved this patient’s prediction most. A patient classified as CVD will show dominant wide red arrows, typically Systolic BP and/or Age, that overpower any blue (negative) contributions. A Non-CVD patient will show wide blue arrows, most commonly from normal blood pressure and an active lifestyle, with modest rightward contributions from age. This intuitive balance makes force plots effective for explaining individual model decisions to clinicians or patients unfamiliar with statistical model internals, bridging the gap between data-driven inference and clinical storytelling.
14. SHAP importance vs impurity (MDI) feature importance
To contextualize SHAP’s ranking, we compare it against the model’s built-in impurity importance (mean decrease in impurity, MDI). scikit-learn’s feature_importances_ is conventionally called Gini importance, and the code below labels the column that way, but for a gradient boosting classifier the base learners are regression trees fitted to pseudo-residuals, so the criterion is squared error (friedman_mse by default) rather than the Gini index of a classification forest:
| Metric | Definition | Computed from |
|---|---|---|
| Mean abs(SHAP) | Average absolute SHAP value across test instances: average impact on model output magnitude | Model output (log-odds) |
| Impurity importance (MDI) | Proportion of total node impurity reduction attributable to each feature across all trees | Tree split structure |
Known limitations of impurity importance (MDI):
- Cardinality bias: continuous features with many possible split points (e.g., Systolic BP) have more opportunities to appear in splits and tend to be overweighted.
- Disconnected from output: MDI is computed from the tree split structure, with the model output playing no part. It adds up the impurity reduction achieved at each split, weighted by the number of samples reaching the node, and never measures how much a feature changes predictions on held-out data.
- No directionality: MDI cannot indicate whether a feature increases or decreases the prediction.
SHAP importance is generally preferred for interpretation because it is derived directly from the model’s output function, measured on held-out data, and decomposes each prediction exactly.
# Mean absolute SHAP value per feature (computed on the 10,000-sample test subset)
mean_abs_shap = np.abs(shap_vals).mean(axis=0)
summary_df = (
pd.DataFrame(
{
"Feature": features_short,
"Mean |SHAP|": mean_abs_shap,
"Gini Importance": gb_clf.feature_importances_,
}
)
.sort_values("Mean |SHAP|", ascending=False)
)
print("Mean absolute SHAP vs Gini importance (sorted by SHAP):")
print(summary_df.to_string(index=False))Mean absolute SHAP vs Gini importance (sorted by SHAP):
Feature Mean |SHAP| Gini Importance
Systolic BP 0.828423 0.686283
Age 0.286570 0.131229
Cholesterol 0.195615 0.073295
Body Mass Index 0.147117 0.078425
Diastolic BP 0.106140 0.021533
Physical Activity 0.076715 0.009235Interpreting the results:
Both methods agree that Systolic BP is overwhelmingly the feature this model leans on most, but the rank ordering is not identical across both metrics. The two agree on the top two and the bottom two ranks and disagree on the middle pair, where mean absolute SHAP places Cholesterol above Body Mass Index and MDI reverses them. Agreement where it does occur is a consistency check on the model’s feature importance structure, and falls well short of proving it robust.
However, MDI assigns even greater relative dominance to Systolic BP (68.6% of total, a larger share than it takes in a normalized SHAP ranking), partly due to its cardinality bias toward continuous features with many split points.
Notably, Diastolic BP and Physical Activity gain the most relative importance in the move from MDI to SHAP, though every feature except Systolic BP gains some. This is consistent with MDI’s tendency to credit features that reduce impurity at many splits, whether or not they move the output. Two candidate explanations sit close to hand, that Physical Activity is binary and that Diastolic BP is strongly correlated with Systolic BP, which takes the splits first, but the notebook tests neither.
Both metrics are computed from the same fitted model, so their agreement is a consistency check and carries no independent evidence for the overall importance ranking. What the SHAP ranking does bear on, separately, is the model’s clinical plausibility. The top-ranked features, systolic blood pressure, age, and cholesterol, are precisely the risk factors most strongly supported by epidemiological evidence for cardiovascular disease. This alignment between data-driven SHAP values and established clinical knowledge is a plausibility check. External validation would need an external dataset. The notebook’s subject is the explanation itself, and held-out accuracy, AUC and calibration are the measures that speak to the quality of the model being explained.
Key Takeaways
A SHAP value records how far a feature moved this model’s output from the base value. What would follow if the feature were changed lies outside its scope, because the object being explained is the model rather than the disease.
For the gradient boosted classifier here the explainer returns the raw decision function, so the base value of −0.022 and every contribution are log-odds. Probability is the logistic transform of the total, and probabilities themselves do not add.
Ranking features by average attribution magnitude drops the sign, so a feature that raises the prediction for some patients and lowers it for others still ranks high. The beeswarm beside the bars is what restores direction and spread.
The definition averages over every coalition of features, one model evaluation per subset, so cost grows exponentially with the feature count. TreeSHAP exploits the tree structure to return exact values for tree ensembles in polynomial time.
Systolic and diastolic pressure carry overlapping information, and every SHAP value rests on a convention for the features left out of a coalition. No background sample is passed here, so that convention comes from the fitted trees themselves.
Since both are computed from the same fitted model, their agreement is not independent evidence. They do not fully agree here in any case: SHAP places cholesterol above body mass index, and MDI reverses the pair.
Data & License
Cardiovascular Disease dataset (about 70,000 patient records) published on Kaggle by Aidan (user “colewelkins”). Source: kaggle.com/datasets/colewelkins/cardiovascular-disease.License: Open Database License (ODbL) v1.0 for the database and Database Contents License (DbCL) v1.0 for its contents; used with attribution. ODbL 1.0 · DbCL 1.0 This dataset appears to be a re-upload; the license shown is the one declared on Kaggle, and its upstream chain of rights is unverified.
© 2025 Philip Sarajlic. All rights reserved for the article’s original text and figures; the third-party data remains under the license shown above.
Code examples in this article are licensed under the Common Public Attribution License Version 1.0 (CPAL-1.0), an OSI-approved copyleft license based on the Mozilla Public License 1.1. Initial Developer: Philip Sarajlic.
Attribution required by CPAL Exhibit B: © 2025 Philip Sarajlic · “Based on code by Philip Sarajlic” · philipsarajlic.com · no graphic image. This attribution must be displayed in Larger Works.
Modifications must be released in Source Code form under CPAL-1.0. Making the code usable by anyone other than you over a network is External Deployment under the license and is treated as distribution, so the Source Code must be made available to those users.
Full text: opensource.org/license/cpal-1-0 (SPDX identifier CPAL-1.0)
