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

import numpy as np
from sklearn.inspection import permutation_importance, PartialDependenceDisplay
import lime.lime_tabular
import shap
 
# Model-agnostic methods work with ANY model
# All they need is a predict() or predict_proba() function
 
class AnyModel:
    """This could be ANY model - the methods don't care"""
    def __init__(self, model):
        self.model = model
    
    def predict(self, X):
        return self.model.predict(X)
    
    def predict_proba(self, X):
        return self.model.predict_proba(X)
 
# Wrap your model (could be sklearn, tensorflow, pytorch, custom, etc.)
wrapped_model = AnyModel(your_model)
 
# ============================================
# Method 1: Permutation Importance (Global, Model-Agnostic)
# ============================================
# Only needs: model.predict() or model.score()
print("📊 Permutation Importance")
result = permutation_importance(
    wrapped_model.model,  # Any model with predict
    X_test, 
    y_test,
    n_repeats=30,
    scoring='accuracy',
    random_state=42
)
 
for name, importance in sorted(
    zip(feature_names, result.importances_mean), 
    key=lambda x: -x[1]
)[:5]:
    print(f"  {name}: {importance:.4f}")
 
# ============================================
# Method 2: LIME (Local, Model-Agnostic)
# ============================================
# Only needs: model.predict_proba()
print("
📊 LIME Local Explanation")
lime_explainer = lime.lime_tabular.LimeTabularExplainer(
    X_train,
    feature_names=feature_names,
    class_names=['Class 0', 'Class 1'],
    mode='classification'
)
 
# Explain any single instance
instance = X_test[0]
explanation = lime_explainer.explain_instance(
    instance,
    wrapped_model.predict_proba,  # Just needs predict_proba
    num_features=5
)
 
print(f"Instance prediction: {wrapped_model.predict([instance])[0]}")
for feature, weight in explanation.as_list():
    print(f"  {feature}: {weight:+.4f}")
 
# ============================================
# Method 3: KernelSHAP (Local+Global, Model-Agnostic)
# ============================================
# Only needs: model.predict() or model.predict_proba()
print("
📊 KernelSHAP Explanation")
kernel_explainer = shap.KernelExplainer(
    wrapped_model.predict_proba,  # Just needs callable
    shap.sample(X_train, 100)     # Background samples
)
 
# Note: This is SLOW compared to model-specific SHAP variants
shap_values = kernel_explainer.shap_values(X_test[:10], nsamples=500)
 
print("Top features by mean |SHAP|:")
importance = np.abs(shap_values[1]).mean(axis=0)
for name, imp in sorted(zip(feature_names, importance), key=lambda x: -x[1])[:5]:
    print(f"  {name}: {imp:.4f}")
 
# ============================================
# Key Point: Same code works for ANY model
# ============================================
# To use with a different model, just change the wrapped_model
# Everything else stays exactly the same
 
# Example models this code works with unchanged:
# - RandomForest, GradientBoosting, XGBoost, LightGBM, CatBoost
# - LogisticRegression, SVM, KNN
# - Neural Networks (sklearn, keras, pytorch)
# - Custom models, API-based models, ensemble models
# - Pre-trained models you can't access internals of

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

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
import statsmodels.api as sm
 
# ============================================
# Method 1: Raw Coefficients
# ============================================
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
 
log_reg = LogisticRegression(random_state=42)
log_reg.fit(X_train_scaled, y_train)
 
print("📊 Logistic Regression Coefficients (Standardized Features)")
print("="*60)
print(f"Intercept: {log_reg.intercept_[0]:.4f}")
print("
Feature coefficients (log-odds scale):")
for name, coef in sorted(
    zip(feature_names, log_reg.coef_[0]), 
    key=lambda x: -abs(x[1])
):
    odds_ratio = np.exp(coef)
    effect = "increases" if coef > 0 else "decreases"
    print(f"  {name}: β={coef:+.4f}, OR={odds_ratio:.3f}")
    print(f"    → 1σ increase in {name} {effect} odds by {abs(odds_ratio-1)*100:.1f}%")
 
# ============================================
# Method 2: Statistical Inference with Statsmodels
# ============================================
print("
📊 Full Statistical Analysis")
print("="*60)
 
# Add constant for intercept
X_train_sm = sm.add_constant(X_train_scaled)
logit_model = sm.Logit(y_train, X_train_sm)
result = logit_model.fit(disp=0)
 
print(result.summary2().tables[1].to_string())
 
# Key outputs:
# - Coef: coefficient value
# - Std.Err: standard error (uncertainty)
# - z: z-statistic (coefficient / std.err)
# - P>|z|: p-value for significance test
# - [0.025, 0.975]: 95% confidence interval
 
# ============================================
# Method 3: Feature Importance Decomposition
# ============================================
print("
📊 Prediction Decomposition for Single Instance")
print("="*60)
 
instance_idx = 0
instance_scaled = X_test_scaled[instance_idx]
prediction = log_reg.predict_proba([instance_scaled])[0, 1]
 
contributions = log_reg.coef_[0] * instance_scaled
log_odds = log_reg.intercept_[0] + contributions.sum()
 
print(f"Instance {instance_idx}:")
print(f"  P(y=1) = {prediction:.4f}")
print(f"  Log-odds = {log_odds:.4f}")
print(f"  Intercept contribution: {log_reg.intercept_[0]:.4f}")
print("
  Feature contributions to log-odds:")
for name, value, contrib in sorted(
    zip(feature_names, instance_scaled, contributions),
    key=lambda x: -abs(x[2])
)[:5]:
    direction = "↑" if contrib > 0 else "↓"
    print(f"    {name} (z={value:.2f}): {contrib:+.4f} {direction}")
 
# Verify: contributions sum to prediction
predicted_log_odds = log_reg.intercept_[0] + contributions.sum()
predicted_prob = 1 / (1 + np.exp(-predicted_log_odds))
print(f"
  Verification: predicted P = {predicted_prob:.4f} (matches)")
 
# ============================================
# Advantages of Linear Interpretability
# ============================================
# 1. EXACT feature contributions (no approximation)
# 2. Statistical inference (p-values, confidence intervals)
# 3. Global interpretation = local interpretation (linear = additive)
# 4. Comparison across features (standardized coefficients)
# 5. Established theory, regulatory acceptance

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

import shap
import numpy as np
from sklearn.tree import DecisionTreeClassifier, export_text
from sklearn.ensemble import GradientBoostingClassifier
import xgboost as xgb
 
# ============================================
# Single Tree: Direct Interpretation
# ============================================
print("📊 Single Decision Tree Interpretation")
print("="*60)
 
tree = DecisionTreeClassifier(max_depth=4, random_state=42)
tree.fit(X_train, y_train)
 
# Print tree rules
print("Decision Rules:")
print(export_text(tree, feature_names=list(feature_names), max_depth=3))
 
# Path for specific prediction
instance_idx = 0
instance = X_test[instance_idx:instance_idx+1]
prediction = tree.predict(instance)[0]
proba = tree.predict_proba(instance)[0]
 
# Get decision path
path = tree.decision_path(instance)
node_indicators = path.toarray()[0]
active_nodes = np.where(node_indicators)[0]
 
print(f"
Prediction path for instance {instance_idx}:")
print(f"  Final prediction: {prediction} (P={proba[prediction]:.3f})")
print("  Decision path:")
for node_id in active_nodes:
    if tree.tree_.children_left[node_id] != -1:  # Not a leaf
        feat = feature_names[tree.tree_.feature[node_id]]
        threshold = tree.tree_.threshold[node_id]
        value = instance[0, tree.tree_.feature[node_id]]
        direction = "≤" if value <= threshold else ">"
        print(f"    Node {node_id}: {feat} = {value:.2f} {direction} {threshold:.2f}")
    else:
        print(f"    Leaf {node_id}: predict class {np.argmax(tree.tree_.value[node_id])}")
 
# ============================================
# Tree Ensemble: TreeSHAP
# ============================================
print("
📊 TreeSHAP for Gradient Boosting")
print("="*60)
 
# Train ensemble
gb_model = GradientBoostingClassifier(n_estimators=100, random_state=42)
gb_model.fit(X_train, y_train)
 
# TreeSHAP - MUCH faster than KernelSHAP
import time
explainer = shap.TreeExplainer(gb_model)
 
start = time.time()
shap_values = explainer.shap_values(X_test[:100])
tree_time = time.time() - start
print(f"TreeSHAP time for 100 instances: {tree_time:.2f}s")
 
# Compare with KernelSHAP (would be ~100x slower)
print(f"Estimated KernelSHAP time: ~{tree_time * 100:.0f}s")
 
# Global importance from TreeSHAP
print("
Global Feature Importance (mean |SHAP|):")
importance = np.abs(shap_values[1]).mean(axis=0)
for name, imp in sorted(zip(feature_names, importance), key=lambda x: -x[1])[:5]:
    print(f"  {name}: {imp:.4f}")
 
# ============================================
# Interaction Effects
# ============================================
print("
📊 SHAP Interaction Analysis")
print("="*60)
 
# Compute interaction values (more expensive)
interaction_values = explainer.shap_interaction_values(X_test[:50])
 
# Sum absolute interaction values
n_features = len(feature_names)
interaction_matrix = np.abs(interaction_values[1]).mean(axis=0)
 
print("Top feature interactions:")
pairs = []
for i in range(n_features):
    for j in range(i+1, n_features):
        pairs.append((feature_names[i], feature_names[j], interaction_matrix[i, j]))
 
for f1, f2, strength in sorted(pairs, key=lambda x: -x[2])[:5]:
    print(f"  {f1} × {f2}: {strength:.4f}")
    
# Interpretation: high values mean these features interact
# Their joint effect differs from sum of individual effects

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

import torch
import torch.nn as nn
import numpy as np
 
# ============================================
# Setup: Simple neural network
# ============================================
class SimpleNN(nn.Module):
    def __init__(self, input_dim, hidden_dim, output_dim):
        super().__init__()
        self.fc1 = nn.Linear(input_dim, hidden_dim)
        self.fc2 = nn.Linear(hidden_dim, hidden_dim)
        self.fc3 = nn.Linear(hidden_dim, output_dim)
        self.relu = nn.ReLU()
    
    def forward(self, x):
        x = self.relu(self.fc1(x))
        x = self.relu(self.fc2(x))
        return self.fc3(x)
 
model = SimpleNN(input_dim=20, hidden_dim=64, output_dim=2)
model.load_state_dict(torch.load('model.pt'))
model.eval()
 
# ============================================
# Method 1: Simple Gradient Saliency
# ============================================
def simple_gradient_saliency(model, input_tensor, target_class):
    """Compute gradient of output w.r.t. input"""
    input_tensor = input_tensor.clone().requires_grad_(True)
    output = model(input_tensor)
    
    # Compute gradient of target class output w.r.t. input
    model.zero_grad()
    output[0, target_class].backward()
    
    # Gradient magnitude as saliency
    saliency = input_tensor.grad.abs()
    return saliency.detach().numpy()
 
instance = torch.tensor(X_test[0:1], dtype=torch.float32)
target = model(instance).argmax().item()
saliency = simple_gradient_saliency(model, instance, target)
 
print("📊 Gradient Saliency")
for name, sal in sorted(zip(feature_names, saliency[0]), key=lambda x: -x[1])[:5]:
    print(f"  {name}: {sal:.4f}")
 
# ============================================
# Method 2: Integrated Gradients
# ============================================
def integrated_gradients(model, input_tensor, baseline, target_class, steps=100):
    """Compute integrated gradients from baseline to input"""
    # Generate interpolated inputs
    scaled_inputs = [baseline + (float(i) / steps) * (input_tensor - baseline) 
                     for i in range(steps + 1)]
    scaled_inputs = torch.cat(scaled_inputs, dim=0)
    scaled_inputs.requires_grad_(True)
    
    # Forward pass
    outputs = model(scaled_inputs)
    
    # Compute gradients
    model.zero_grad()
    outputs[:, target_class].sum().backward()
    gradients = scaled_inputs.grad
    
    # Integrate (average gradients × input difference)
    avg_gradients = gradients.mean(dim=0, keepdim=True)
    integrated_grads = (input_tensor - baseline) * avg_gradients
    
    return integrated_grads.detach().numpy()
 
baseline = torch.zeros_like(instance)  # Common choice: zero baseline
ig = integrated_gradients(model, instance, baseline, target)
 
print("
📊 Integrated Gradients")
for name, val in sorted(zip(feature_names, ig[0]), key=lambda x: -abs(x[1]))[:5]:
    direction = "↑" if val > 0 else "↓"
    print(f"  {name}: {val:+.4f} {direction}")
 
# ============================================
# Method 3: Using Captum Library
# ============================================
from captum.attr import IntegratedGradients, DeepLift, LayerGradCam
 
# Integrated Gradients with Captum
ig_captum = IntegratedGradients(model)
attributions = ig_captum.attribute(instance, baseline, target=target)
 
print("
📊 Integrated Gradients (Captum)")
for name, attr in sorted(zip(feature_names, attributions[0].numpy()), 
                         key=lambda x: -abs(x[1]))[:5]:
    print(f"  {name}: {attr:+.4f}")
 
# DeepLift
dl = DeepLift(model)
dl_attr = dl.attribute(instance, baseline, target=target)
 
print("
📊 DeepLift")
for name, attr in sorted(zip(feature_names, dl_attr[0].numpy()), 
                         key=lambda x: -abs(x[1]))[:5]:
    print(f"  {name}: {attr:+.4f}")
 
# ============================================
# CNN: GradCAM for image classification
# ============================================
# For CNNs, GradCAM provides spatial localization
# from captum.attr import LayerGradCam
# 
# grad_cam = LayerGradCam(cnn_model, cnn_model.layer4)  # Target conv layer
# attributions = grad_cam.attribute(image_tensor, target=predicted_class)
# # Result: heatmap showing which image regions matter for this prediction

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

import shap
import numpy as np
from functools import lru_cache
import joblib
 
class ProductionExplainer:
    """
    Production-ready explainer with caching, batching, and efficiency optimizations.
    """
    
    def __init__(self, model, X_background, feature_names):
        self.model = model
        self.feature_names = feature_names
        
        # Pre-compute SHAP explainer (expensive)
        print("Initializing TreeExplainer (one-time cost)...")
        self.explainer = shap.TreeExplainer(model)
        
        # Pre-compute expected value
        self.expected_value = self.explainer.expected_value
        
        # Cache for frequently requested explanations
        self.cache = {}
    
    def _instance_hash(self, instance):
        """Create hashable key for instance caching"""
        return tuple(instance.flatten().tolist())
    
    def explain_single(self, instance, use_cache=True):
        """Explain single instance with optional caching"""
        key = self._instance_hash(instance)
        
        if use_cache and key in self.cache:
            return self.cache[key]
        
        shap_values = self.explainer.shap_values(instance.reshape(1, -1))
        
        # Handle binary classification
        if isinstance(shap_values, list):
            shap_values = shap_values[1]  # Positive class
        
        result = {
            'shap_values': shap_values[0],
            'base_value': float(self.expected_value[1] if isinstance(self.expected_value, list) 
                               else self.expected_value),
            'features': {
                name: {'value': float(instance[i]), 'shap': float(shap_values[0, i])}
                for i, name in enumerate(self.feature_names)
            },
            'top_features': self._get_top_features(shap_values[0], instance)
        }
        
        if use_cache:
            self.cache[key] = result
        
        return result
    
    def explain_batch(self, instances, progress=False):
        """Efficient batch explanation"""
        shap_values = self.explainer.shap_values(instances)
        
        if isinstance(shap_values, list):
            shap_values = shap_values[1]
        
        results = []
        for i in range(len(instances)):
            results.append({
                'shap_values': shap_values[i],
                'top_features': self._get_top_features(shap_values[i], instances[i])
            })
        
        return results
    
    def _get_top_features(self, shap_values, instance, n=5):
        """Get top n features by SHAP magnitude"""
        sorted_idx = np.argsort(-np.abs(shap_values))[:n]
        return [
            {
                'feature': self.feature_names[i],
                'value': float(instance[i]),
                'shap': float(shap_values[i]),
                'direction': 'positive' if shap_values[i] > 0 else 'negative'
            }
            for i in sorted_idx
        ]
    
    def get_global_importance(self, X_sample):
        """Compute global feature importance from sample"""
        shap_values = self.explainer.shap_values(X_sample)
        
        if isinstance(shap_values, list):
            shap_values = shap_values[1]
        
        importance = np.abs(shap_values).mean(axis=0)
        
        return {
            name: float(imp)
            for name, imp in sorted(
                zip(self.feature_names, importance),
                key=lambda x: -x[1]
            )
        }
    
    def save(self, path):
        """Serialize explainer for deployment"""
        joblib.dump({
            'model': self.model,
            'feature_names': self.feature_names,
            'expected_value': self.expected_value
        }, path)
 
# Usage
explainer = ProductionExplainer(model, X_train[:100], feature_names)
 
# Single prediction with explanation
instance = X_test[0]
prediction = model.predict([instance])[0]
explanation = explainer.explain_single(instance)
 
print(f"Prediction: {prediction}")
print("Top contributing features:")
for feat in explanation['top_features']:
    print(f"  {feat['feature']}: {feat['shap']:+.4f} ({feat['direction']})")