What is the difference between SHAP and LIME?

SHAP (SHapley Additive exPlanations) uses game theory to assign each feature a consistent contribution score across all predictions — global and local explanations are mathematically consistent. LIME (Local Interpretable Model-agnostic Explanations) fits a simple linear model around a single prediction — faster but less consistent across predictions. SHAP is the standard choice today.

Do I need to explain every model?

In regulated industries (credit, healthcare, hiring), explanation is often legally required. In other contexts, interpretability matters for debugging (why is the model wrong on these cases?), feature discovery (which inputs actually matter?), and trust-building. Tree-based models via SHAP are cheap to explain. Neural networks are more expensive but still feasible with SHAP.

Model Interpretability | SkillByExample

Real-World Scenario

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

Feature Importances (Tree Models)

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

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

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

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

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

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

SHAP Values

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

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

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

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

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

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

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

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

SHAP for Any Model (KernelExplainer)

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

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

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

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

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

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

Partial Dependence Plots (PDP)

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

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

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

# Compute partial dependence for the top feature
pd_result = partial_dependence(
    model, X_train,
    features=[0],          # feature index (income)
    kind="average",        # "average" = PDP, "individual" = ICE curves
    percentiles=(0.05, 0.95),
    grid_resolution=50,
)
values   = pd_result["grid_values"][0]
avg_pred = pd_result["average"][0]

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

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

LIME for Individual Predictions

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

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

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

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

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

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