Python 64.6%
TeX 23.3%
HTML 11.6%
Makefile 0.5%
1"""2shap_analysis.py3================4SHAP-based interpretability analysis for hedonic pricing models.56Provides SHAP value computation, publication-quality summary / dependence /7interaction / waterfall plots, partial-dependence plots, accumulated local8effects (ALE), and hedonic premium interpretation. Gracefully degrades to9permutation importance and manual PDP when the ``shap`` library is not10installed.1112Author: Simon-Pierre Boucher13Date: 2026-05-2014"""1516from __future__ import annotations1718import logging19import os20import warnings21from pathlib import Path22from typing import Any, Dict, List, Optional, Tuple, Union2324import matplotlib25import matplotlib.pyplot as plt26import numpy as np27import pandas as pd28from sklearn.inspection import PartialDependenceDisplay, partial_dependence29from sklearn.metrics import r2_score3031logger = logging.getLogger(__name__)3233# Publication-quality defaults ------------------------------------------------34matplotlib.rcParams.update(35 {36 "figure.dpi": 150,37 "savefig.dpi": 300,38 "font.family": "serif",39 "font.size": 11,40 "axes.titlesize": 13,41 "axes.labelsize": 12,42 "xtick.labelsize": 10,43 "ytick.labelsize": 10,44 "legend.fontsize": 10,45 "figure.figsize": (10, 7),46 "savefig.bbox": "tight",47 "savefig.pad_inches": 0.1,48 }49)5051# Try importing shap ----------------------------------------------------------52_SHAP_AVAILABLE = False53try:54 import shap5556 _SHAP_AVAILABLE = True57except ImportError:58 logger.info(59 "shap library not installed. Falling back to permutation importance "60 "and manual PDP implementations."61 )626364# ---------------------------------------------------------------------------65# Utility helpers66# ---------------------------------------------------------------------------676869def _ensure_dir(path: str) -> None:70 """Create parent directories for *path* if they do not exist."""71 Path(path).parent.mkdir(parents=True, exist_ok=True)727374def _save_figure(fig: plt.Figure, save_path: str) -> None:75 """Save a matplotlib figure in both PNG and PDF."""76 _ensure_dir(save_path)77 base, _ = os.path.splitext(save_path)78 fig.savefig(f"{base}.png", format="png")79 fig.savefig(f"{base}.pdf", format="pdf")80 logger.info("Figure saved: %s.png / .pdf", base)81 plt.close(fig)828384def _to_dataframe(X: np.ndarray, feature_names: List[str]) -> pd.DataFrame:85 """Convert an ndarray to a DataFrame with proper column names."""86 return pd.DataFrame(X, columns=feature_names)878889# ---------------------------------------------------------------------------90# Permutation-importance fallback (when shap is unavailable)91# ---------------------------------------------------------------------------929394def _permutation_importance(95 model: Any,96 X: np.ndarray,97 y: Optional[np.ndarray] = None,98 n_repeats: int = 10,99 random_state: int = 42,100) -> np.ndarray:101 """Compute permutation importance as a (n_features,) array of mean102 absolute importance values.103104 When *y* is not provided, variance of prediction changes is used instead.105 """106 rng = np.random.RandomState(random_state)107 baseline = model.predict(X)108109 if y is not None:110 baseline_score = r2_score(y, baseline)111 else:112 baseline_score = None113114 importances = np.zeros(X.shape[1])115 for j in range(X.shape[1]):116 scores: List[float] = []117 for _ in range(n_repeats):118 X_perm = X.copy()119 X_perm[:, j] = rng.permutation(X_perm[:, j])120 preds_perm = model.predict(X_perm)121 if y is not None and baseline_score is not None:122 scores.append(baseline_score - r2_score(y, preds_perm))123 else:124 scores.append(float(np.mean(np.abs(baseline - preds_perm))))125 importances[j] = float(np.mean(scores))126 return importances127128129# ---------------------------------------------------------------------------130# SHAP value computation131# ---------------------------------------------------------------------------132133134def compute_shap_values(135 model: Any,136 X: np.ndarray,137 feature_names: List[str],138 max_samples: int = 500,139) -> Union["shap.Explanation", np.ndarray]:140 """Compute SHAP values for the given model and data.141142 Strategy:143 1. If ``shap`` is available, use ``TreeExplainer`` for tree-based models;144 otherwise fall back to ``KernelExplainer`` (on a subsample for speed).145 2. If ``shap`` is **not** installed, return a permutation-importance146 matrix shaped ``(n_samples, n_features)`` as a rough substitute.147148 Parameters149 ----------150 model : estimator151 Fitted sklearn-compatible model.152 X : np.ndarray153 Feature matrix.154 feature_names : list[str]155 Column names.156 max_samples : int157 Maximum number of background samples for ``KernelExplainer``.158159 Returns160 -------161 shap.Explanation | np.ndarray162 SHAP values object or permutation-based importance matrix.163 """164 if _SHAP_AVAILABLE:165 model_type = type(model).__name__166 tree_types = {167 "RandomForestRegressor",168 "GradientBoostingRegressor",169 "XGBRegressor",170 "LGBMRegressor",171 "CatBoostRegressor",172 }173174 X_df = _to_dataframe(X, feature_names)175176 if model_type in tree_types:177 try:178 explainer = shap.TreeExplainer(model)179 sv = explainer(X_df)180 logger.info(181 "SHAP TreeExplainer computed for %s (%d samples).",182 model_type,183 X.shape[0],184 )185 return sv186 except Exception:187 logger.warning(188 "TreeExplainer failed for %s — falling back to KernelExplainer.",189 model_type,190 )191192 # KernelExplainer fallback.193 n_bg = min(max_samples, X.shape[0])194 bg = shap.sample(X_df, n_bg)195 explainer = shap.KernelExplainer(model.predict, bg)196 sv = explainer(X_df)197 logger.info(198 "SHAP KernelExplainer computed for %s (%d samples, %d bg).",199 model_type,200 X.shape[0],201 n_bg,202 )203 return sv204 else:205 # Permutation-importance fallback: tile per-feature importance across206 # samples to produce an (n_samples, n_features) matrix.207 logger.warning(208 "shap not installed — computing permutation importance as fallback."209 )210 imp = _permutation_importance(model, X)211 # Normalise to approximate per-sample attribution.212 imp_matrix = np.tile(imp, (X.shape[0], 1))213 return imp_matrix214215216# ---------------------------------------------------------------------------217# SHAP Plots218# ---------------------------------------------------------------------------219220221def _extract_values_array(222 shap_values: Union[Any, np.ndarray],223) -> np.ndarray:224 """Get a plain (n_samples, n_features) ndarray from shap Explanation or raw array."""225 if hasattr(shap_values, "values"):226 return np.asarray(shap_values.values)227 return np.asarray(shap_values)228229230def plot_shap_summary(231 shap_values: Union[Any, np.ndarray],232 X: np.ndarray,233 feature_names: List[str],234 save_path: str,235 max_display: int = 20,236) -> None:237 """SHAP beeswarm summary plot (or bar chart fallback).238239 Parameters240 ----------241 shap_values : shap.Explanation | np.ndarray242 SHAP values.243 X : np.ndarray244 Feature matrix (for colouring dots).245 feature_names : list[str]246 Column names.247 save_path : str248 Output file path (without extension; both .png and .pdf are saved).249 max_display : int250 Maximum number of features to show.251 """252 X_df = _to_dataframe(X, feature_names)253254 if _SHAP_AVAILABLE and hasattr(shap_values, "values"):255 fig, ax = plt.subplots(figsize=(10, max(7, max_display * 0.35)))256 shap.summary_plot(257 shap_values,258 X_df,259 max_display=max_display,260 show=False,261 )262 fig = plt.gcf()263 fig.suptitle("SHAP Feature Importance", y=1.02, fontweight="bold")264 _save_figure(fig, save_path)265 else:266 # Fallback: horizontal bar chart of mean |values|.267 vals = _extract_values_array(shap_values)268 mean_abs = np.mean(np.abs(vals), axis=0)269 order = np.argsort(mean_abs)[::-1][:max_display]270271 fig, ax = plt.subplots(figsize=(10, max(6, max_display * 0.35)))272 ax.barh(273 range(len(order)),274 mean_abs[order][::-1],275 color="#1f77b4",276 edgecolor="white",277 )278 ax.set_yticks(range(len(order)))279 ax.set_yticklabels([feature_names[i] for i in order][::-1])280 ax.set_xlabel("Mean |SHAP value| (or permutation importance)")281 ax.set_title("Feature Importance", fontweight="bold")282 fig.tight_layout()283 _save_figure(fig, save_path)284285286def plot_shap_dependence(287 shap_values: Union[Any, np.ndarray],288 X: np.ndarray,289 feature: str,290 feature_names: List[str],291 save_path: str,292 interaction_feature: Optional[str] = None,293) -> None:294 """SHAP dependence plot for a single feature.295296 Parameters297 ----------298 shap_values : shap.Explanation | np.ndarray299 SHAP values.300 X : np.ndarray301 Feature matrix.302 feature : str303 Target feature name.304 feature_names : list[str]305 All column names.306 save_path : str307 Output file path.308 interaction_feature : str | None309 Optional second feature for colouring.310 """311 if feature not in feature_names:312 raise ValueError(f"Feature '{feature}' not found in feature_names.")313314 j = feature_names.index(feature)315 vals = _extract_values_array(shap_values)316317 fig, ax = plt.subplots(figsize=(8, 6))318319 if interaction_feature and interaction_feature in feature_names:320 k = feature_names.index(interaction_feature)321 scatter = ax.scatter(322 X[:, j],323 vals[:, j],324 c=X[:, k],325 cmap="coolwarm",326 alpha=0.6,327 s=10,328 edgecolors="none",329 )330 cbar = fig.colorbar(scatter, ax=ax)331 cbar.set_label(interaction_feature)332 else:333 ax.scatter(X[:, j], vals[:, j], alpha=0.4, s=10, color="#1f77b4", edgecolors="none")334335 ax.set_xlabel(feature)336 ax.set_ylabel(f"SHAP value for {feature}")337 ax.set_title(f"SHAP Dependence: {feature}", fontweight="bold")338 ax.axhline(0, color="grey", linewidth=0.8, linestyle="--")339 fig.tight_layout()340 _save_figure(fig, save_path)341342343def plot_shap_interaction(344 shap_values: Union[Any, np.ndarray],345 X: np.ndarray,346 feature1: str,347 feature2: str,348 save_path: str,349 feature_names: Optional[List[str]] = None,350) -> None:351 """Interaction plot between two features coloured by SHAP contribution.352353 Parameters354 ----------355 shap_values : shap.Explanation | np.ndarray356 SHAP values.357 X : np.ndarray358 Feature matrix.359 feature1, feature2 : str360 Feature names.361 save_path : str362 Output file path.363 feature_names : list[str] | None364 Column names.365 """366 if feature_names is None:367 raise ValueError("feature_names must be provided.")368369 j1 = feature_names.index(feature1)370 j2 = feature_names.index(feature2)371 vals = _extract_values_array(shap_values)372373 fig, ax = plt.subplots(figsize=(8, 6))374 scatter = ax.scatter(375 X[:, j1],376 X[:, j2],377 c=vals[:, j1] + vals[:, j2],378 cmap="RdBu_r",379 alpha=0.6,380 s=12,381 edgecolors="none",382 )383 cbar = fig.colorbar(scatter, ax=ax)384 cbar.set_label("Combined SHAP value")385 ax.set_xlabel(feature1)386 ax.set_ylabel(feature2)387 ax.set_title(388 f"SHAP Interaction: {feature1} x {feature2}",389 fontweight="bold",390 )391 fig.tight_layout()392 _save_figure(fig, save_path)393394395# ---------------------------------------------------------------------------396# SHAP Importance Table397# ---------------------------------------------------------------------------398399400def get_shap_importance(401 shap_values: Union[Any, np.ndarray],402 feature_names: List[str],403) -> pd.DataFrame:404 """Compute mean |SHAP| per feature and return a ranked DataFrame.405406 Parameters407 ----------408 shap_values : shap.Explanation | np.ndarray409 SHAP values.410 feature_names : list[str]411 Column names.412413 Returns414 -------415 pd.DataFrame416 Columns: ``Feature``, ``mean_abs_shap``, sorted descending.417 """418 vals = _extract_values_array(shap_values)419 mean_abs = np.mean(np.abs(vals), axis=0)420 df = (421 pd.DataFrame({"Feature": feature_names, "mean_abs_shap": mean_abs})422 .sort_values("mean_abs_shap", ascending=False)423 .reset_index(drop=True)424 )425 return df426427428# ---------------------------------------------------------------------------429# Partial Dependence Plots (sklearn-based)430# ---------------------------------------------------------------------------431432433def plot_partial_dependence(434 model: Any,435 X: np.ndarray,436 features: List[str],437 feature_names: List[str],438 save_path: str,439 grid_resolution: int = 50,440) -> None:441 """Partial dependence plots for selected features using sklearn.442443 Parameters444 ----------445 model : estimator446 Fitted model.447 X : np.ndarray448 Feature matrix.449 features : list[str]450 Feature names to plot (max 9 per figure).451 feature_names : list[str]452 All column names.453 save_path : str454 Output file path.455 grid_resolution : int456 Number of grid points.457 """458 indices = [feature_names.index(f) for f in features if f in feature_names]459 if not indices:460 logger.warning("No valid features found for partial dependence plot.")461 return462463 X_df = _to_dataframe(X, feature_names)464465 n_feats = len(indices)466 ncols = min(3, n_feats)467 nrows = int(np.ceil(n_feats / ncols))468469 fig, axes = plt.subplots(470 nrows,471 ncols,472 figsize=(5 * ncols, 4 * nrows),473 constrained_layout=True,474 )475 if n_feats == 1:476 axes = np.array([axes])477 axes = np.atleast_2d(axes)478479 for ax_idx, feat_idx in enumerate(indices):480 row, col = divmod(ax_idx, ncols)481 ax = axes[row, col]482483 try:484 pd_result = partial_dependence(485 model,486 X_df,487 features=[feat_idx],488 grid_resolution=grid_resolution,489 kind="average",490 )491 grid_vals = pd_result["grid_values"][0]492 avg_pred = pd_result["average"][0]493494 ax.plot(grid_vals, avg_pred, color="#1f77b4", linewidth=2)495 ax.set_xlabel(feature_names[feat_idx])496 ax.set_ylabel("Partial dependence")497 ax.set_title(feature_names[feat_idx], fontweight="bold")498 ax.grid(True, alpha=0.3)499 except Exception as exc:500 logger.warning("PDP failed for feature %s: %s", feature_names[feat_idx], exc)501 ax.set_visible(False)502503 # Hide unused axes.504 for ax_idx in range(n_feats, nrows * ncols):505 row, col = divmod(ax_idx, ncols)506 axes[row, col].set_visible(False)507508 fig.suptitle("Partial Dependence Plots", fontweight="bold", fontsize=14)509 _save_figure(fig, save_path)510511512# ---------------------------------------------------------------------------513# Accumulated Local Effects (ALE) — manual implementation514# ---------------------------------------------------------------------------515516517def plot_ale(518 model: Any,519 X: np.ndarray,520 feature: str,521 feature_names: List[str],522 save_path: str,523 n_bins: int = 40,524) -> None:525 """Accumulated Local Effects plot for a single feature.526527 ALE is computed manually following Apley & Zhu (2020): the feature range528 is binned, and the average finite-difference in predictions within each529 bin is accumulated.530531 Parameters532 ----------533 model : estimator534 Fitted model.535 X : np.ndarray536 Feature matrix.537 feature : str538 Target feature name.539 feature_names : list[str]540 All column names.541 save_path : str542 Output file path.543 n_bins : int544 Number of quantile bins.545 """546 j = feature_names.index(feature)547 x_feat = X[:, j].copy()548549 # Quantile bin edges.550 percentiles = np.linspace(0, 100, n_bins + 1)551 edges = np.unique(np.percentile(x_feat, percentiles))552 if len(edges) < 3:553 logger.warning(554 "Feature '%s' has too few unique values for ALE — skipping.",555 feature,556 )557 return558559 n_actual_bins = len(edges) - 1560 ale_values = np.zeros(n_actual_bins)561 bin_counts = np.zeros(n_actual_bins)562563 for b in range(n_actual_bins):564 lower, upper = edges[b], edges[b + 1]565 if b < n_actual_bins - 1:566 mask = (x_feat >= lower) & (x_feat < upper)567 else:568 mask = (x_feat >= lower) & (x_feat <= upper)569570 if mask.sum() == 0:571 continue572573 X_low = X[mask].copy()574 X_high = X[mask].copy()575 X_low[:, j] = lower576 X_high[:, j] = upper577578 pred_diff = model.predict(X_high) - model.predict(X_low)579 ale_values[b] = float(np.mean(pred_diff))580 bin_counts[b] = mask.sum()581582 # Accumulate.583 ale_cumulative = np.cumsum(ale_values)584 # Centre so that the weighted mean is zero.585 total = bin_counts.sum()586 if total > 0:587 ale_cumulative -= np.sum(ale_cumulative * bin_counts) / total588589 bin_centres = (edges[:-1] + edges[1:]) / 2.0590591 fig, ax = plt.subplots(figsize=(8, 5))592 ax.plot(bin_centres, ale_cumulative, color="#d62728", linewidth=2)593 ax.fill_between(bin_centres, ale_cumulative, alpha=0.15, color="#d62728")594 ax.axhline(0, color="grey", linewidth=0.8, linestyle="--")595 ax.set_xlabel(feature)596 ax.set_ylabel("ALE")597 ax.set_title(f"Accumulated Local Effects: {feature}", fontweight="bold")598 ax.grid(True, alpha=0.3)599 fig.tight_layout()600 _save_figure(fig, save_path)601602603# ---------------------------------------------------------------------------604# Hedonic Premium Interpretation605# ---------------------------------------------------------------------------606607608def interpret_hedonic_premiums(609 shap_values: Union[Any, np.ndarray],610 X: np.ndarray,611 feature_names: List[str],612) -> pd.DataFrame:613 """Convert SHAP values into interpretable hedonic premiums.614615 For each feature the function reports:616617 * **mean_shap** -- average SHAP contribution (signed; positive = premium,618 negative = discount).619 * **mean_abs_shap** -- average absolute contribution (importance).620 * **std_shap** -- variability of the SHAP contribution across listings.621 * **median_shap** -- median contribution.622 * **premium_pct** -- percentage of listings for which the feature adds623 a positive premium.624 * **feature_mean** -- mean raw feature value (for context).625 * **feature_std** -- standard deviation of the raw feature.626627 Parameters628 ----------629 shap_values : shap.Explanation | np.ndarray630 SHAP values.631 X : np.ndarray632 Feature matrix.633 feature_names : list[str]634 Column names.635636 Returns637 -------638 pd.DataFrame639 One row per feature, sorted by ``mean_abs_shap`` descending.640 """641 vals = _extract_values_array(shap_values)642 records: List[Dict[str, Any]] = []643644 for j, fname in enumerate(feature_names):645 sv = vals[:, j]646 records.append(647 {648 "Feature": fname,649 "mean_shap": float(np.mean(sv)),650 "mean_abs_shap": float(np.mean(np.abs(sv))),651 "std_shap": float(np.std(sv)),652 "median_shap": float(np.median(sv)),653 "premium_pct": float(np.mean(sv > 0) * 100),654 "feature_mean": float(np.mean(X[:, j])),655 "feature_std": float(np.std(X[:, j])),656 }657 )658659 df = (660 pd.DataFrame(records)661 .sort_values("mean_abs_shap", ascending=False)662 .reset_index(drop=True)663 )664 logger.info("Hedonic premiums computed for %d features.", len(feature_names))665 return df666667668# ---------------------------------------------------------------------------669# Single-Listing Explanation670# ---------------------------------------------------------------------------671672673def explain_single_listing(674 model: Any,675 shap_values: Union[Any, np.ndarray],676 X: np.ndarray,677 idx: int,678 feature_names: List[str],679 save_path: Optional[str] = None,680 top_n: int = 15,681) -> pd.DataFrame:682 """Explain the predicted price of a single listing via a waterfall chart.683684 Parameters685 ----------686 model : estimator687 Fitted model (used only to obtain the base prediction).688 shap_values : shap.Explanation | np.ndarray689 SHAP values.690 X : np.ndarray691 Feature matrix.692 idx : int693 Row index of the listing to explain.694 feature_names : list[str]695 Column names.696 save_path : str | None697 If provided, save the waterfall plot.698 top_n : int699 Number of top features to display.700701 Returns702 -------703 pd.DataFrame704 Per-feature SHAP contributions for the listing, sorted by absolute705 magnitude.706 """707 vals = _extract_values_array(shap_values)708 sv_row = vals[idx]709710 explanation_df = (711 pd.DataFrame(712 {713 "Feature": feature_names,714 "Feature_Value": X[idx],715 "SHAP_Value": sv_row,716 }717 )718 .assign(abs_shap=lambda d: np.abs(d["SHAP_Value"]))719 .sort_values("abs_shap", ascending=False)720 .drop(columns="abs_shap")721 .reset_index(drop=True)722 )723724 predicted = float(model.predict(X[idx : idx + 1])[0])725 base_value: float = predicted - float(np.sum(sv_row))726 logger.info(727 "Listing %d: predicted=%.2f base=%.2f sum(SHAP)=%.2f",728 idx,729 predicted,730 base_value,731 np.sum(sv_row),732 )733734 # Waterfall / bar chart.735 if _SHAP_AVAILABLE and hasattr(shap_values, "values"):736 try:737 fig = plt.figure(figsize=(10, 7))738 single_exp = shap_values[idx]739 shap.waterfall_plot(single_exp, max_display=top_n, show=False)740 fig = plt.gcf()741 fig.suptitle(742 f"Listing {idx} Explanation (predicted=${predicted:,.0f})",743 y=1.02,744 fontweight="bold",745 )746 if save_path:747 _save_figure(fig, save_path)748 else:749 plt.close(fig)750 return explanation_df751 except Exception:752 logger.warning("SHAP waterfall_plot failed — using manual bar chart.")753754 # Manual waterfall-style horizontal bar chart.755 top = explanation_df.head(top_n)756 fig, ax = plt.subplots(figsize=(10, max(6, top_n * 0.4)))757 colours = ["#d62728" if v < 0 else "#2ca02c" for v in top["SHAP_Value"]]758 y_pos = range(len(top) - 1, -1, -1)759 ax.barh(760 y_pos,761 top["SHAP_Value"].values,762 color=colours,763 edgecolor="white",764 height=0.7,765 )766 labels = [767 f"{row.Feature} = {row.Feature_Value:.2g}"768 for _, row in top.iterrows()769 ]770 ax.set_yticks(list(y_pos))771 ax.set_yticklabels(labels)772 ax.axvline(0, color="grey", linewidth=0.8)773 ax.set_xlabel("SHAP value (impact on prediction)")774 ax.set_title(775 f"Listing {idx} Explanation | Base: ${base_value:,.0f} | Predicted: ${predicted:,.0f}",776 fontweight="bold",777 )778 fig.tight_layout()779780 if save_path:781 _save_figure(fig, save_path)782 else:783 plt.close(fig)784785 return explanation_df786