shap_analysis.py

"""
SHAP analysis + publication-quality plots.
Requires: data/nba_props_dataset.csv + data/models/gb_{stat}.pkl
Outputs:  data/plots/shap_*.png, data/plots/accuracy_comparison.png
"""
import os, json, pickle, logging, warnings
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")   # headless
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import shap

warnings.filterwarnings("ignore")
logging.basicConfig(level=logging.INFO, format="%(asctime)s %(message)s")
logger = logging.getLogger(__name__)

DATA_DIR  = os.path.join(os.path.dirname(__file__), "data")
MODEL_DIR = os.path.join(DATA_DIR, "models")
PLOT_DIR  = os.path.join(DATA_DIR, "plots")
os.makedirs(PLOT_DIR, exist_ok=True)

FEATURE_LABELS = {
    "l1":           "L1 avg",
    "l3":           "L3 avg",
    "l5":           "L5 avg",
    "l10":          "L10 avg (line proxy)",
    "l20":          "L20 avg",
    "season_avg":   "Season avg",
    "momentum_l3":  "L3 momentum",
    "momentum_l5":  "L5 momentum",
    "acceleration": "Acceleration",
    "cv_5":         "CV (5-game)",
    "cv_10":        "CV (10-game)",
    "min_l5":       "Avg minutes (L5)",
    "min_trend":    "Minutes trend",
    "hr_10":        "Hit rate (L10)",
    "pts_l5":       "Points L5",
    "reb_l5":       "Rebounds L5",
    "ast_l5":       "Assists L5",
    "is_playoff":   "Playoff game",
    "game_num":     "Career game #",
}

STATS = ["pts", "ast", "reb", "threes", "blk", "stl"]
STAT_LABELS = {
    "pts": "Points", "ast": "Assists", "reb": "Rebounds",
    "threes": "3-Pointers Made", "blk": "Blocks", "stl": "Steals"
}


def load_model(stat):
    path = os.path.join(MODEL_DIR, f"gb_{stat}.pkl")
    if not os.path.exists(path):
        return None
    with open(path, "rb") as f:
        return pickle.load(f)


def shap_summary_per_stat(stat, df):
    pkg = load_model(stat)
    if pkg is None:
        logger.warning(f"No model found for {stat}, skipping")
        return

    model        = pkg["model"]
    feature_cols = pkg["feature_cols"]
    df_stat      = df[df["stat"] == stat].copy()

    if len(df_stat) < 100:
        return

    X = df_stat[feature_cols].fillna(0).values
    # Use a sample of 2000 for SHAP (speed)
    idx = np.random.choice(len(X), min(2000, len(X)), replace=False)
    X_sample = X[idx]

    logger.info(f"Computing SHAP for {stat} ({len(X_sample)} samples)...")
    explainer   = shap.TreeExplainer(model)
    shap_values = explainer.shap_values(X_sample)
    # For binary classifiers, shap_values may be list[2]; take class 1
    if isinstance(shap_values, list):
        sv = shap_values[1]
    else:
        sv = shap_values

    # Plot
    fig, ax = plt.subplots(figsize=(9, 6))
    feature_names = [FEATURE_LABELS.get(c, c) for c in feature_cols]

    mean_abs = np.abs(sv).mean(axis=0)
    order    = np.argsort(mean_abs)[::-1][:15]

    colors = plt.cm.RdYlGn(np.linspace(0.15, 0.85, len(order)))
    ax.barh(
        [feature_names[i] for i in order[::-1]],
        mean_abs[order[::-1]],
        color=colors,
        edgecolor="white", linewidth=0.5,
    )
    ax.set_xlabel("Mean |SHAP value|", fontsize=11)
    ax.set_title(f"Feature Importance — {STAT_LABELS.get(stat, stat)}\n"
                 f"(GradientBoosting, SHAP, n={len(X_sample):,})", fontsize=12)
    ax.spines[["top", "right"]].set_visible(False)
    plt.tight_layout()
    out = os.path.join(PLOT_DIR, f"shap_{stat}.png")
    fig.savefig(out, dpi=150, bbox_inches="tight")
    plt.close(fig)
    logger.info(f"  Saved → {out}")


def accuracy_comparison_plot():
    results_path = os.path.join(DATA_DIR, "model_results.csv")
    if not os.path.exists(results_path):
        logger.warning("No model_results.csv found — run ml_pipeline.py first")
        return

    df = pd.read_csv(results_path)

    # Use last test season only for clean comparison
    last_split = df[df["split"].str.endswith(str(max([int(s.split("test")[1]) for s in df["split"].unique()])))].copy()

    models    = ["LogReg", "SVM", "RandomForest", "GradBoost", "Ensemble"]
    stats_ord = ["pts", "ast", "reb", "threes", "blk", "stl"]

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Panel 1: accuracy by model (bar chart)
    ax1  = axes[0]
    summary = last_split.groupby("model")["accuracy"].mean().reindex(models).dropna()
    colors  = ["#3498db", "#e74c3c", "#2ecc71", "#f39c12", "#9b59b6"]
    bars    = ax1.bar(summary.index, summary.values * 100, color=colors[:len(summary)],
                      edgecolor="white", linewidth=0.8)
    ax1.axhline(52.4, color="gray",  linestyle="--", linewidth=1, label="Break-even (52.4%)")
    ax1.axhline(56.0, color="orange", linestyle="--", linewidth=1, label="Rule-based baseline (56%)")
    ax1.set_ylabel("Accuracy (%)", fontsize=11)
    ax1.set_title("Model Accuracy Comparison\n(Walk-forward, final test season)", fontsize=11)
    ax1.set_ylim(48, 72)
    ax1.legend(fontsize=9)
    ax1.spines[["top", "right"]].set_visible(False)
    for bar, val in zip(bars, summary.values):
        ax1.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.3,
                 f"{val*100:.1f}%", ha="center", va="bottom", fontsize=9)

    # Panel 2: GradBoost accuracy by stat
    ax2  = axes[1]
    gb   = last_split[last_split["model"] == "GradBoost"].set_index("stat")
    stat_accs = [gb.loc[s, "accuracy"] * 100 if s in gb.index else None for s in stats_ord]
    stat_names = [STAT_LABELS.get(s, s) for s in stats_ord]
    colors2 = ["#2ecc71" if (a and a >= 60) else "#e67e22" if (a and a >= 56) else "#e74c3c"
               for a in stat_accs]
    valid   = [(n, a, c) for n, a, c in zip(stat_names, stat_accs, colors2) if a is not None]
    if valid:
        ns, accs, cols = zip(*valid)
        bars2 = ax2.bar(ns, accs, color=cols, edgecolor="white", linewidth=0.8)
        ax2.axhline(52.4, color="gray", linestyle="--", linewidth=1)
        ax2.set_ylabel("Accuracy (%)", fontsize=11)
        ax2.set_title("GradientBoosting Accuracy by Stat\n(Walk-forward, final test season)", fontsize=11)
        ax2.set_ylim(48, 72)
        ax2.spines[["top", "right"]].set_visible(False)
        ax2.tick_params(axis="x", rotation=15)
        for bar, val in zip(bars2, accs):
            ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.3,
                     f"{val:.1f}%", ha="center", va="bottom", fontsize=9)

    plt.tight_layout()
    out = os.path.join(PLOT_DIR, "accuracy_comparison.png")
    fig.savefig(out, dpi=150, bbox_inches="tight")
    plt.close(fig)
    logger.info(f"Accuracy comparison plot saved → {out}")


def feature_correlation_heatmap(stat="pts"):
    dataset_path = os.path.join(DATA_DIR, "nba_props_dataset.csv")
    if not os.path.exists(dataset_path):
        return

    df     = pd.read_csv(dataset_path)
    df_s   = df[df["stat"] == stat]
    corr   = df_s[["l1","l3","l5","l10","season_avg","momentum_l3","cv_10","hr_10","label"]].corr()

    fig, ax = plt.subplots(figsize=(8, 7))
    im = ax.imshow(corr.values, cmap="RdBu_r", vmin=-1, vmax=1)
    ax.set_xticks(range(len(corr.columns)))
    ax.set_yticks(range(len(corr.columns)))
    labs = [FEATURE_LABELS.get(c, c) for c in corr.columns]
    ax.set_xticklabels(labs, rotation=35, ha="right", fontsize=9)
    ax.set_yticklabels(labs, fontsize=9)
    for i in range(len(corr)):
        for j in range(len(corr)):
            ax.text(j, i, f"{corr.values[i,j]:.2f}", ha="center", va="center", fontsize=7)
    plt.colorbar(im, ax=ax, fraction=0.04)
    ax.set_title(f"Feature Correlation — {STAT_LABELS.get(stat, stat)}", fontsize=12)
    plt.tight_layout()
    out = os.path.join(PLOT_DIR, f"correlation_{stat}.png")
    fig.savefig(out, dpi=150, bbox_inches="tight")
    plt.close(fig)
    logger.info(f"Correlation heatmap saved → {out}")


def run_all():
    dataset_path = os.path.join(DATA_DIR, "nba_props_dataset.csv")
    if not os.path.exists(dataset_path):
        logger.error("Dataset not found. Run dataset_builder.py first.")
        return

    df = pd.read_csv(dataset_path)
    logger.info(f"Loaded {len(df):,} rows for SHAP analysis")

    for stat in STATS:
        shap_summary_per_stat(stat, df)

    accuracy_comparison_plot()
    feature_correlation_heatmap("pts")
    feature_correlation_heatmap("ast")

    logger.info(f"\nAll plots saved to {PLOT_DIR}/")


if __name__ == "__main__":
    run_all()