Model Interpretability and Explainability Questions

import numpy as np
import pandas as pd

def _is_categorical(series):
    return pd.api.types.is_categorical_dtype(series) or series.dtype == object

def partial_dependence_plot(model, X, features, grid_resolution=50, kind='average', batch_size=10000):
    """
    Compute PDP values for 1 or 2 features.
    - model: scikit-learn-like estimator (predict or predict_proba)
    - X: pandas DataFrame or numpy array (if numpy, features should be indices)
    - features: single feature or tuple/list of 2 features (column names or indices)
    - grid_resolution: number of points for continuous features
    - kind: 'average' (averaged prediction) - placeholder to match signature
    Returns: (grid, pdp_values)
      - grid: array of grid values (1D for single feature, tuple(mesh_x, mesh_y) for two)
      - pdp_values: averaged prediction per grid point (1D array for single feature, 2D array for two)
    """
    # Normalize X to DataFrame for convenience
    is_numpy = isinstance(X, np.ndarray)
    if is_numpy:
        X = pd.DataFrame(X)
    else:
        X = X.copy()
    # Normalize features input
    if isinstance(features, (list, tuple)) and len(features) == 2:
        feat = list(features)
        two = True
    else:
        feat = [features] if not isinstance(features, (list, tuple)) else features
        two = False if len(feat) == 1 else True

    # Build grid values per feature
    grids = []
    is_cat = []
    for f in feat:
        col = X[f] if f in X.columns else X.iloc[:, f]
        if _is_categorical(col):
            vals = np.asarray(pd.Series(col).astype('category').cat.categories)
            is_cat.append(True)
        else:
            vals = np.linspace(col.min(), col.max(), grid_resolution)
            is_cat.append(False)
        grids.append(vals)

    # Create list of grid points
    if not two:
        grid_vals = grids[0]
        G = len(grid_vals)
        # We'll create G modified datasets each with N rows, stack them -> (G*N, n_features)
        N = len(X)
        modified_blocks = []
        for v in grid_vals:
            Xmod = X.copy()
            Xmod[feat[0]] = v
            modified_blocks.append(Xmod)
        stacked = pd.concat(modified_blocks, ignore_index=True)
        # Predict in batches to avoid memory/time blowup
        preds = []
        predict_func = None
        if hasattr(model, "predict_proba"):
            def predict_func(arr):
                proba = model.predict_proba(arr)
                # choose positive class if binary, else average across classes weighted? we'll pick class 1 if exists
                if proba.shape[1] >= 2:
                    return proba[:, 1]
                return proba.mean(axis=1)
        else:
            predict_func = lambda arr: model.predict(arr)
        # Batch prediction
        total = len(stacked)
        for start in range(0, total, batch_size):
            end = start + batch_size
            preds.append(predict_func(stacked.iloc[start:end]))
        preds = np.concatenate([np.asarray(p) for p in preds])
        preds = preds.reshape(G, N)  # each row: grid point predictions across N rows
        pdp_vals = preds.mean(axis=1)
        return grid_vals, pdp_vals

    else:
        gx, gy = grids[0], grids[1]
        Gx, Gy = len(gx), len(gy)
        N = len(X)
        # prepare modified stacked datasets for all (gx,gy) pairs
        modified_blocks = []
        for vx in gx:
            for vy in gy:
                Xmod = X.copy()
                Xmod[feat[0]] = vx
                Xmod[feat[1]] = vy
                modified_blocks.append(Xmod)
        stacked = pd.concat(modified_blocks, ignore_index=True)
        # predict in batches
        if hasattr(model, "predict_proba"):
            def predict_func(arr):
                proba = model.predict_proba(arr)
                if proba.shape[1] >= 2:
                    return proba[:, 1]
                return proba.mean(axis=1)
        else:
            predict_func = lambda arr: model.predict(arr)
        preds = []
        total = len(stacked)
        for start in range(0, total, batch_size):
            end = start + batch_size
            preds.append(predict_func(stacked.iloc[start:end]))
        preds = np.concatenate([np.asarray(p) for p in preds])
        preds = preds.reshape(Gx * Gy, N)
        pdp_vals = preds.mean(axis=1).reshape(Gx, Gy)
        # return meshgrid to facilitate plotting
        mesh_x, mesh_y = np.meshgrid(gx, gy, indexing='xy')
        return (mesh_x, mesh_y), pdp_vals.T  # transpose so dims align with meshgrid order

import shap
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Preconditions: model is a trained tree-based classifier (RandomForestClassifier or XGBClassifier)
# X_test is a pandas.DataFrame (same columns/order as training set preprocessing)
# If you used a preprocessing pipeline, apply it to raw test data first and keep feature names.

explainer = shap.TreeExplainer(model)            # fast, model-aware explainer for trees
shap_values = explainer.shap_values(X_test)     # note: shape depends on model & shap version

# Normalize returned shape: for binary classifiers SKLearn often returns list [neg, pos]
if isinstance(shap_values, list):
    # shap_values[class_index] has shape (n_samples, n_features)
    shap_vals_class1 = shap_values[1]           # class 1 (positive)
else:
    # XGBoost sklearn API sometimes returns array (n_samples, n_features) already for prob
    shap_vals_class1 = shap_values

# 1) Summary plot (global feature importance for class 1)
shap.summary_plot(shap_vals_class1, X_test, show=False)
plt.title("SHAP summary (class 1)")
plt.tight_layout()
plt.show()

# 2) Extract top 5 features contributing to class 1 predictions
mean_abs_shap = np.abs(shap_vals_class1).mean(axis=0)
feat_importance = pd.Series(mean_abs_shap, index=X_test.columns).sort_values(ascending=False)
top5 = feat_importance.head(5)
print("Top 5 features by mean |SHAP| for class 1:\n", top5)