Local vs Global Interpretability

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105

import shap
import lime
import lime.lime_tabular
import numpy as np
from sklearn.ensemble import GradientBoostingClassifier
 
# Train model
model = GradientBoostingClassifier(n_estimators=100, random_state=42)
model.fit(X_train, y_train)
 
# Select instance to explain
instance_idx = 42
instance = X_test[instance_idx:instance_idx+1]
prediction = model.predict(instance)[0]
probability = model.predict_proba(instance)[0, 1]
 
print(f"Instance {instance_idx}: Prediction = {prediction}, P(positive) = {probability:.3f}")
print("="*60)
 
# Method 1: SHAP Local Explanation
explainer_shap = shap.TreeExplainer(model)
shap_values = explainer_shap.shap_values(instance)
 
print("
📊 SHAP Local Explanation:")
print(f"Base value (expected): {explainer_shap.expected_value[1]:.4f}")
print(f"Sum of SHAP values: {shap_values[1][0].sum():.4f}")
print(f"Prediction log-odds: {explainer_shap.expected_value[1] + shap_values[1][0].sum():.4f}")
print("
Top 5 feature contributions:")
contributions = list(zip(feature_names, instance[0], shap_values[1][0]))
for name, val, shap_val in sorted(contributions, key=lambda x: -abs(x[2]))[:5]:
    direction = "↑" if shap_val > 0 else "↓"
    print(f"  {name}={val:.2f}: {shap_val:+.4f} {direction}")
 
# Method 2: LIME Local Explanation
explainer_lime = lime.lime_tabular.LimeTabularExplainer(
    X_train, 
    feature_names=feature_names,
    class_names=['Negative', 'Positive'],
    mode='classification'
)
 
lime_exp = explainer_lime.explain_instance(
    instance[0], 
    model.predict_proba,
    num_features=5
)
 
print("
📊 LIME Local Explanation:")
print(f"Local model intercept: {lime_exp.intercept[1]:.4f}")
print(f"Local prediction: {lime_exp.local_pred[0]:.4f}")
print("
Top 5 local linear coefficients:")
for feature_rule, weight in lime_exp.as_list()[:5]:
    direction = "↑" if weight > 0 else "↓"
    print(f"  {feature_rule}: {weight:+.4f} {direction}")
 
# Method 3: Counterfactual Explanation (using DiCE)
try:
    import dice_ml
    
    # Setup DiCE
    data = dice_ml.Data(
        dataframe=train_df,
        continuous_features=continuous_features,
        outcome_name='target'
    )
    dice_model = dice_ml.Model(model=model, backend='sklearn')
    exp = dice_ml.Dice(data, dice_model, method='random')
    
    # Generate counterfactuals
    query_instance = pd.DataFrame([instance[0]], columns=feature_names)
    cf = exp.generate_counterfactuals(
        query_instance,
        total_CFs=3,
        desired_class='opposite'
    )
    
    print("
📊 Counterfactual Explanations:")
    print("To flip the prediction, consider these changes:")
    cf.visualize_as_dataframe()
except ImportError:
    print("
(DiCE not installed - skipping counterfactual example)")
 
# Comparison: Do methods agree?
print("
" + "="*60)
print("📊 Method Agreement Analysis")
print("="*60)
 
# Get top 3 features from each method
shap_top3 = set([name for name, _, _ in sorted(contributions, key=lambda x: -abs(x[2]))[:3]])
lime_top3 = set([feat.split()[0] for feat, _ in lime_exp.as_list()[:3]])
 
overlap = shap_top3.intersection(lime_top3)
print(f"SHAP top 3: {shap_top3}")
print(f"LIME top 3: {lime_top3}")
print(f"Agreement: {len(overlap)}/3 features overlap")
 
# Note: Some disagreement is expected due to different methodologies
# High disagreement may indicate instability or complex local behavior

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

import shap
import numpy as np
import matplotlib.pyplot as plt
from sklearn.inspection import permutation_importance, PartialDependenceDisplay
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.tree import DecisionTreeClassifier
 
# Train complex model
model = GradientBoostingClassifier(n_estimators=100, random_state=42)
model.fit(X_train, y_train)
 
# ============================================
# Method 1: Permutation Importance
# ============================================
print("📊 Permutation Feature Importance")
print("="*50)
 
perm_imp = permutation_importance(
    model, X_test, y_test, 
    n_repeats=30, 
    random_state=42
)
 
for name, imp, std in sorted(
    zip(feature_names, perm_imp.importances_mean, perm_imp.importances_std),
    key=lambda x: -x[1]
):
    print(f"  {name}: {imp:.4f} ± {std:.4f}")
 
# ============================================
# Method 2: SHAP Global Importance
# ============================================
print("
📊 SHAP Global Feature Importance")
print("="*50)
 
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)
 
# Mean absolute SHAP value per feature
shap_importance = np.abs(shap_values[1]).mean(axis=0)
for name, imp in sorted(zip(feature_names, shap_importance), key=lambda x: -x[1]):
    print(f"  {name}: {imp:.4f}")
 
# SHAP summary plot (beeswarm) - shows distribution of effects
shap.summary_plot(shap_values[1], X_test, feature_names=feature_names, show=False)
plt.title("SHAP Summary Plot")
plt.tight_layout()
plt.savefig("shap_summary.png", dpi=150)
plt.close()
 
# ============================================
# Method 3: Partial Dependence Plots
# ============================================
print("
📊 Generating Partial Dependence Plots...")
 
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
top_features = [name for name, _ in sorted(
    zip(feature_names, shap_importance), key=lambda x: -x[1]
)[:6]]
 
for ax, feature in zip(axes.flatten(), top_features):
    feature_idx = list(feature_names).index(feature)
    PartialDependenceDisplay.from_estimator(
        model, X_test, [feature_idx],
        feature_names=feature_names,
        ax=ax
    )
    ax.set_title(f"PDP: {feature}")
 
plt.tight_layout()
plt.savefig("pdp_plots.png", dpi=150)
plt.close()
print("  Saved to pdp_plots.png")
 
# ============================================
# Method 4: Global Surrogate Model
# ============================================
print("
📊 Global Surrogate Model Analysis")
print("="*50)
 
# Use black-box predictions as target
bb_predictions = model.predict(X_train)
 
# Train interpretable surrogate
surrogate = DecisionTreeClassifier(max_depth=5, random_state=42)
surrogate.fit(X_train, bb_predictions)
 
# Evaluate fidelity
train_fidelity = (surrogate.predict(X_train) == bb_predictions).mean()
test_predictions = model.predict(X_test)
test_fidelity = (surrogate.predict(X_test) == test_predictions).mean()
 
print(f"  Surrogate fidelity (train): {train_fidelity:.3f}")
print(f"  Surrogate fidelity (test): {test_fidelity:.3f}")
 
if test_fidelity > 0.85:
    print("  ✓ High fidelity - surrogate explanations are reasonably faithful")
    # Export tree rules
    from sklearn.tree import export_text
    print("
  Decision tree rules (simplified global explanation):")
    print(export_text(surrogate, feature_names=list(feature_names), max_depth=3))
else:
    print("  ⚠ Low fidelity - interpret surrogate with caution")
 
# ============================================
# Comparing Global Methods
# ============================================
print("
📊 Cross-Method Agreement")
print("="*50)
 
# Get top 5 from each method
perm_top5 = set([name for name, _, _ in sorted(
    zip(feature_names, perm_imp.importances_mean, perm_imp.importances_std),
    key=lambda x: -x[1]
)[:5]])
 
shap_top5 = set([name for name, _ in sorted(
    zip(feature_names, shap_importance), key=lambda x: -x[1]
)[:5]])
 
print(f"  Permutation top 5: {perm_top5}")
print(f"  SHAP top 5: {shap_top5}")
print(f"  Agreement: {len(perm_top5.intersection(shap_top5))}/5 features")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85

import shap
import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
 
# Compute SHAP values for all test instances
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)[1]  # For positive class
 
# ============================================
# Aggregation 1: Mean Absolute SHAP (Global Importance)
# ============================================
print("📊 Global Feature Importance (Mean |SHAP|)")
global_importance = np.abs(shap_values).mean(axis=0)
for name, imp in sorted(zip(feature_names, global_importance), key=lambda x: -x[1]):
    print(f"  {name}: {imp:.4f}")
 
# ============================================
# Aggregation 2: SHAP Distribution Analysis
# ============================================
print("
📊 SHAP Value Distribution Analysis")
print("(Heterogeneity indicator: std/mean)")
 
for i, name in enumerate(feature_names):
    mean_abs = np.abs(shap_values[:, i]).mean()
    std = np.abs(shap_values[:, i]).std()
    heterogeneity = std / mean_abs if mean_abs > 0 else 0
    print(f"  {name}: mean={mean_abs:.4f}, std={std:.4f}, CV={heterogeneity:.2f}")
 
# High CV = heterogeneous effect (feature matters more for some instances)
# Low CV = homogeneous effect (feature matters similarly across instances)
 
# ============================================
# Aggregation 3: Subgroup Analysis
# ============================================
print("
📊 Subgroup SHAP Analysis (by demographic)")
 
# Assume we have a demographic variable
# This reveals if model treats groups differently
for group in ['A', 'B']:
    mask = X_test_df['demographic'] == group
    group_shap = shap_values[mask]
    group_importance = np.abs(group_shap).mean(axis=0)
    
    print(f"
  Group {group} (n={mask.sum()}):")
    for name, imp in sorted(zip(feature_names, group_importance), key=lambda x: -x[1])[:3]:
        print(f"    {name}: {imp:.4f}")
 
# If top features differ between groups, model uses different decision logic
 
# ============================================
# Aggregation 4: Explanation Clustering
# ============================================
print("
📊 Explanation Clustering (finding decision modes)")
 
# Cluster instances by their SHAP profiles
n_clusters = 3
kmeans = KMeans(n_clusters=n_clusters, random_state=42)
cluster_labels = kmeans.fit_predict(shap_values)
 
for c in range(n_clusters):
    mask = cluster_labels == c
    cluster_shap = shap_values[mask]
    mean_shap = cluster_shap.mean(axis=0)
    
    print(f"
  Cluster {c} (n={mask.sum()}, {100*mask.mean():.1f}% of data):")
    # Top features for this cluster
    top_idx = np.argsort(-np.abs(mean_shap))[:3]
    for idx in top_idx:
        direction = "↑" if mean_shap[idx] > 0 else "↓"
        print(f"    {feature_names[idx]}: {mean_shap[idx]:+.4f} {direction}")
    
    # Characterize cluster
    cluster_X = X_test[mask]
    print(f"    Avg prediction: {model.predict_proba(cluster_X)[:, 1].mean():.3f}")
 
# Each cluster represents a distinct "explanation mode"
# Cluster 0 might be "approved due to high income"
# Cluster 1 might be "denied due to poor credit history"
# etc.