In this tutorial, we implement SHAP workflows as a practical framework for interpreting machine learning models beyond basic feature-importance plots. We start by training tree-based models and then compare different SHAP explainers, including Tree, Exact, Permutation, and Kernel methods, to understand how accuracy and runtime change across model-aware and model-agnostic approaches. We also examine how maskers affect explanations when features are correlated, how interaction values reveal pairwise feature effects, and how link functions alter interpretation between the log-odds and probability spaces. Also, we use Owen values, cohort testing, SHAP-based feature selection, drift monitoring, and custom black-box explanations to build a complete interpretability workflow that can run directly in Google Colab. Copy CodeCopiedUse a different Browser !pip install -q –upgrade shap xgboost transformers import warnings, time, numpy as np, pandas as pd, matplotlib.pyplot as plt from scipy import stats from scipy.cluster import hierarchy warnings.filterwarnings(“ignore”) import shap, xgboost as xgb from sklearn.datasets import fetch_california_housing, load_breast_cancer from sklearn.model_selection import train_test_split from sklearn.metrics import roc_auc_score, r2_score shap.initjs() np.random.seed(42) print(f”SHAP: {shap.__version__}n”) housing = fetch_california_housing() X = pd.DataFrame(housing.data, columns=housing.feature_names) y = pd.Series(housing.target, name=”MedHouseVal”) reg = xgb.XGBRegressor(n_estimators=300, max_depth=5, learning_rate=0.05, subsample=0.9, random_state=42, n_jobs=-1).fit(X_tr, y_tr) print(f”Housing regressor R² = {reg.score(X_te, y_te):.3f}”) def reg_predict(X): return reg.predict(np.asarray(X)) We install the required libraries and import the core tools for SHAP, XGBoost, statistics, visualization, and model evaluation. We load the California housing dataset and train an XGBoost regression model. We also define a clean prediction wrapper so that SHAP can explain the model without running into compatibility issues with bound model methods. Copy CodeCopiedUse a different Browser print(“n” + “=”*72) print(“PART 1: Explainer comparison — correctness & speed”) print(“=”*72) X_sample = X_te.iloc[:25] bg_small = shap.sample(X_tr, 50, random_state=42) def _wrap_kernel(expl, X, bg_mean): vals = expl.shap_values(X, nsamples=200, silent=True) return shap.Explanation(values=vals, base_values=np.full(len(X), bg_mean), data=X.values, feature_names=X.columns.tolist()) runs = {} t0 = time.time(); tree_expl = shap.TreeExplainer(reg); sv_tree = tree_expl(X_sample) runs[“Tree (exact, model-aware)”] = (sv_tree, time.time() – t0) t0 = time.time() sv_exact = shap.Explainer(reg_predict, bg_small, algorithm=”exact”)(X_sample) runs[“Exact (model-agnostic)”] = (sv_exact, time.time() – t0) t0 = time.time() sv_perm = shap.Explainer(reg_predict, bg_small, algorithm=”permutation”)(X_sample) runs[“Permutation”] = (sv_perm, time.time() – t0) t0 = time.time() ke = shap.KernelExplainer(reg_predict, shap.sample(X_tr, 50, random_state=42).values) sv_kern = _wrap_kernel(ke, X_sample, ke.expected_value) runs[“Kernel”] = (sv_kern, time.time() – t0) ref = sv_tree.values.flatten() print(f”n{‘Method’:30s} {‘time(s)’:>8s} {‘ρ vs Tree’:>10s} {‘max|Δ|’:>8s}”) for name, (sv, dt) in runs.items(): flat = sv.values.flatten() rho = np.corrcoef(ref, flat)[0, 1] err = np.abs(ref – flat).max() print(f”{name:30s} {dt:8.2f} {rho:10.4f} {err:8.4f}”) print(“nTakeaway: Tree is the only exact + fast option for tree ensembles.”) print(“Exact ≈ Permutation when permutation has enough samples; Kernel is noisier and slowest.”) print(“n” + “=”*72) print(“PART 2: Maskers — Independent vs Partition under correlation”) print(“=”*72) corr = X_tr.corr().abs() top_pair = corr.where(np.triu(np.ones_like(corr, dtype=bool), k=1)) .stack().sort_values(ascending=False).head(3) print(“Top correlated pairs (|ρ|):”) for (a, b), v in top_pair.items(): print(f” {a:10s} {b:10s} |ρ| = {v:.3f}”) masker_ind = shap.maskers.Independent(X_tr, max_samples=100) masker_part = shap.maskers.Partition(X_tr, max_samples=100) sv_ind = shap.Explainer(reg_predict, masker_ind)(X_sample) sv_part = shap.Explainer(reg_predict, masker_part)(X_sample) a, b = top_pair.index[0] print(f”nMean |φ| for top-correlated pair ({a}, {b}):”) print(f” Independent : {a}={np.abs(sv_ind[:,a].values).mean():.4f} {b}={np.abs(sv_ind[:,b].values).mean():.4f}”) print(f” Partition : {a}={np.abs(sv_part[:,a].values).mean():.4f} {b}={np.abs(sv_part[:,b].values).mean():.4f}”) print(“Partition redistributes credit across correlated features (on-manifold semantics).”) fig, axes = plt.subplots(1, 2, figsize=(13, 4)) plt.sca(axes[0]); shap.plots.bar(sv_ind, show=False); axes[0].set_title(“Independent masker”) plt.sca(axes[1]); shap.plots.bar(sv_part, show=False); axes[1].set_title(“Partition masker”) plt.tight_layout(); plt.show() We compare multiple SHAP explainers, including Tree, Exact, Permutation, and Kernel, on the same regression model and sample data. We measure each method by runtime, correlation with TreeExplainer, and maximum attribution difference to understand the trade-off between speed and approximation quality. We then study Independent and Partition maskers to see how correlated features receive different attribution credit under different masking assumptions. Copy CodeCopiedUse a different Browser print(“n” + “=”*72) print(“PART 3: Interaction decomposition”) print(“=”*72) inter = tree_expl.shap_interaction_values(X_te.iloc[:500]) inter_abs = np.abs(inter).mean(0) diag = np.diagonal(inter_abs).copy() off = inter_abs.copy(); np.fill_diagonal(off, 0) main_share = diag.sum() / (diag.sum() + off.sum()) print(f”Total attribution mass: {main_share*100:.1f}% main effects, ” f”{(1-main_share)*100:.1f}% interactions”) pairs = [(X.columns[i], X.columns[j], off[i, j]) for i in range(X.shape[1]) for j in range(i+1, X.shape[1])] pairs.sort(key=lambda t: -t[2]) print(“nTop 5 interaction pairs (mean |φ_ij|):”) for a, b, v in pairs[:5]: print(f” {a:10s} × {b:10s} → {v:.4f}”) fig, ax = plt.subplots(figsize=(7.5, 6)) im = ax.imshow(off, cmap=”viridis”) ax.set_xticks(range(X.shape[1])); ax.set_xticklabels(X.columns, rotation=45, ha=”right”) ax.set_yticks(range(X.shape[1])); ax.set_yticklabels(X.columns) plt.colorbar(im, label=”mean |φ_ij|”); plt.title(“Pairwise interaction strength”) plt.tight_layout(); plt.show() a, b, _ = pairs[0] i, j = X.columns.get_loc(a), X.columns.get_loc(b) xs = X_te.iloc[:500][a].values; cs = X_te.iloc[:500][b].values fig, axes = plt.subplots(1, 2, figsize=(13, 4), sharex=True) axes[0].scatter(xs, inter[:, i, i], c=cs, s=12, cmap=”coolwarm”) axes[0].set_title(f”Main effect of {a}”); axes[0].set_xlabel(a); axes[0].set_ylabel(“φ_{ii}”) sc = axes[1].scatter(xs, 2*inter[:, i, j], c=cs, s=12, cmap=”coolwarm”) axes[1].set_title(f”Interaction {a} × {b}”); axes[1].set_xlabel(a); axes[1].set_ylabel(“2·φ_{ij}”) plt.colorbar(sc, ax=axes[1], label=b); plt.tight_layout(); plt.show() print(“n” + “=”*72) print(“PART 4: Link functions — logit vs probability space”) print(“=”*72) cancer = load_breast_cancer() Xc = pd.DataFrame(cancer.data, columns=cancer.feature_names) yc = pd.Series(cancer.target) clf = xgb.XGBClassifier(n_estimators=300, max_depth=4, learning_rate=0.05, eval_metric=”logloss”, random_state=42).fit(Xc_tr, yc_tr) print(f”AUC = {roc_auc_score(yc_te, clf.predict_proba(Xc_te)[:,1]):.3f}”) expl_logit = shap.TreeExplainer(clf) sv_logit = expl_logit(Xc_te) expl_prob = shap.TreeExplainer(clf, Xc_tr.sample(100, random_state=42), model_output=”probability”) sv_prob = expl_prob(Xc_te) print(f”nSample 0 reconstruction (φ should sum to f – E[f]):”) print(f” log-odds : base + Σφ = {sv_logit.base_values[0] + sv_logit.values[0].sum():+.3f}”) print(f” prob : base + Σφ = {sv_prob.base_values[0] + sv_prob.values[0].sum():.3f} ” f”(model proba = {clf.predict_proba(Xc_te.iloc[[0]])[0,1]:.3f})”) fig, axes = plt.subplots(1, 2, figsize=(15, 5)) plt.sca(axes[0]); shap.plots.waterfall(sv_logit[0], max_display=8, show=False); axes[0].set_title(“Log-odds space”) plt.sca(axes[1]); shap.plots.waterfall(sv_prob[0], max_display=8, show=False); axes[1].set_title(“Probability space”) plt.tight_layout(); plt.show() We calculate SHAP interaction values to separate main feature effects from pairwise interaction effects in the housing model. We identify the strongest interaction pairs and visualize their attribution strength using heatmaps and scatter plots. We then move to a classification task and compare SHAP explanations in log-odds and probability spaces using a breast cancer classifier. Copy CodeCopiedUse a different Browser print(“n” + “=”*72) print(“PART 5: Owen values from a correlation-based feature hierarchy”) print(“=”*72) D = 1 – X_tr.corr().abs().values np.fill_diagonal(D, 0) condensed = D[np.triu_indices_from(D, k=1)] linkage = hierarchy.linkage(condensed, method=”average”) masker_owen = shap.maskers.Partition(X_tr, clustering=linkage, max_samples=100) sv_owen = shap.Explainer(reg_predict, masker_owen)(X_sample) fig, axes = plt.subplots(1, 2, figsize=(14, 4.5)) hierarchy.dendrogram(linkage, labels=X.columns.tolist(), ax=axes[0]) axes[0].set_title(“Feature hierarchy (1 − |ρ|)”) plt.sca(axes[1]); shap.plots.bar(sv_owen.abs.mean(0), show=False) axes[1].set_title(“Owen values (cluster-aware)”) plt.tight_layout(); plt.show() print(“n” + “=”*72) print(“PART 6: Cohort comparison with bootstrap CIs and hypothesis tests”) print(“=”*72) sv_all = tree_expl(X_te) q1, q3 = X_te[“MedInc”].quantile([0.25, 0.75]) low = (X_te[“MedInc”] <= q1).values high = (X_te[“MedInc”] >= q3).values def boot_ci(v, B=1000, seed=0):
