Skills › Scientific Agent Skills › shap › tabular_report.py
tabular_report.py
348 lines · skills/shap/scripts/tabular_report.py at 5f2bc03
1#!/usr/bin/env python3
2"""Generate a validated SHAP 0.52 tabular report from a built-in dataset.
4The script is intentionally self-contained: it downloads no data and loads no
5serialized model. Use it as a modern template for Explanation slicing,
6probability-space TreeExplainer output, additivity validation, and figure
7export.
8"""
10from __future__ import annotations
12import argparse
13import json
14import math
15import platform
16from pathlib import Path
17from typing import Any
19import matplotlib
21matplotlib.use("Agg")
23import matplotlib.pyplot as plt
24import numpy as np
25import pandas as pd
26import shap
27import sklearn
28from sklearn.datasets import load_breast_cancer
29from sklearn.ensemble import RandomForestClassifier
30from sklearn.model_selection import train_test_split
33def positive_int(value: str) -> int:
34 """Parse a strictly positive integer."""
35 parsed = int(value)
36 if parsed <= 0:
37 raise argparse.ArgumentTypeError("value must be greater than zero")
38 return parsed
41def nonnegative_int(value: str) -> int:
42 """Parse a non-negative integer."""
43 parsed = int(value)
44 if parsed < 0:
45 raise argparse.ArgumentTypeError("value must be zero or greater")
46 return parsed
49def positive_float(value: str) -> float:
50 """Parse a finite, strictly positive tolerance."""
51 parsed = float(value)
52 if not math.isfinite(parsed) or parsed <= 0:
53 raise argparse.ArgumentTypeError("value must be finite and greater than zero")
54 return parsed
57def parse_args() -> argparse.Namespace:
58 """Build and parse command-line arguments."""
59 parser = argparse.ArgumentParser(
60 description=(
61 "Train a deterministic random forest on scikit-learn's breast "
62 "cancer dataset and export a validated probability-space SHAP report."
63 )
64 )
65 parser.add_argument(
66 "--output-dir",
67 type=Path,
68 default=Path("shap-report"),
69 help="Directory for CSV, JSON, and PNG outputs (default: shap-report).",
70 )
71 parser.add_argument(
72 "--class-index",
73 type=nonnegative_int,
74 default=1,
75 help="Class probability: 0=malignant, 1=benign (default: 1).",
76 )
77 parser.add_argument(
78 "--background-size",
79 type=positive_int,
80 default=100,
81 help="Maximum number of training rows in the background (default: 100).",
82 )
83 parser.add_argument(
84 "--explain-size",
85 type=positive_int,
86 default=100,
87 help="Maximum number of held-out rows to explain (default: 100).",
88 )
89 parser.add_argument(
90 "--max-display",
91 type=positive_int,
92 default=15,
93 help="Maximum features displayed in plots (default: 15).",
94 )
95 parser.add_argument(
96 "--seed",
97 type=int,
98 default=7,
99 help="Random seed for splitting, fitting, and background sampling.",
100 )
101 parser.add_argument(
102 "--atol",
103 type=positive_float,
104 default=1e-6,
105 help="Absolute tolerance for additive reconstruction (default: 1e-6).",
106 )
107 parser.add_argument(
108 "--rtol",
109 type=positive_float,
110 default=1e-5,
111 help="Relative tolerance for additive reconstruction (default: 1e-5).",
112 )
113 return parser.parse_args()
116def select_output(
117 explanation: shap.Explanation,
118 class_index: int,
119) -> shap.Explanation:
120 """Select one output from a tabular Explanation."""
121 if class_index < 0:
122 raise ValueError("class_index must be zero or greater")
123 values = np.asarray(explanation.values)
125 if values.ndim == 2:
126 if class_index != 0:
127 raise ValueError(
128 "The explainer produced one output; use --class-index 0 or "
129 "configure a model output with multiple classes."
130 )
131 return explanation
133 if values.ndim != 3:
134 raise ValueError(
135 "Expected tabular values shaped (samples, features) or "
136 f"(samples, features, outputs), received {values.shape}."
137 )
139 output_count = values.shape[-1]
140 if class_index >= output_count:
141 raise ValueError(
142 f"--class-index {class_index} is outside {output_count} outputs."
143 )
145 return explanation[..., class_index]
148def save_axis(axis: Any, path: Path, title: str | None = None) -> None:
149 """Save and close the figure associated with a SHAP plot axis."""
150 figure = axis.figure if axis is not None else plt.gcf()
151 if title is not None:
152 figure.suptitle(title)
153 figure.tight_layout()
154 figure.savefig(path, dpi=200, bbox_inches="tight")
155 plt.close(figure)
158def output_name(explanation: shap.Explanation, class_index: int) -> str:
159 """Return a stable display name for the selected output."""
160 name = explanation.output_names
161 if name is None:
162 return f"class_{class_index}"
163 if isinstance(name, str):
164 return name
165 if np.ndim(name) == 0:
166 return str(name)
167 return f"class_{class_index}"
170def main() -> int:
171 """Run the complete example and return a process exit code."""
172 args = parse_args()
173 args.output_dir.mkdir(parents=True, exist_ok=True)
175 dataset = load_breast_cancer(as_frame=True)
176 # Match the float32 input representation used by sklearn tree prediction.
177 X, y = dataset.data.astype(np.float32), dataset.target
178 X_train, X_test, y_train, y_test = train_test_split(
179 X,
180 y,
181 test_size=0.25,
182 stratify=y,
183 random_state=args.seed,
184 )
186 model = RandomForestClassifier(
187 n_estimators=200,
188 min_samples_leaf=3,
189 random_state=args.seed,
190 n_jobs=-1,
191 )
192 model.fit(X_train, y_train)
194 if args.class_index >= len(model.classes_):
195 raise ValueError(
196 f"--class-index {args.class_index} is outside model classes "
197 f"{model.classes_.tolist()}."
198 )
200 background_size = min(args.background_size, len(X_train))
201 background = shap.sample(
202 X_train,
203 background_size,
204 random_state=args.seed,
205 )
206 X_explain = X_test.iloc[: min(args.explain_size, len(X_test))].copy()
208 explainer = shap.TreeExplainer(
209 model,
210 # Passing a bare frame silently resamples backgrounds larger than 100.
211 data=shap.maskers.Independent(background, max_samples=len(background)),
212 feature_perturbation="interventional",
213 model_output="probability",
214 )
215 all_outputs = explainer(X_explain)
216 all_outputs.output_names = [
217 str(dataset.target_names[int(label)]) for label in model.classes_
218 ]
219 explanation = select_output(all_outputs, args.class_index)
221 values = np.asarray(explanation.values, dtype=float)
222 base_values = np.asarray(explanation.base_values, dtype=float)
223 predictions = model.predict_proba(X_explain)[:, args.class_index]
224 reconstructed = base_values + values.sum(axis=1)
225 errors = reconstructed - predictions
226 max_abs_error = float(np.max(np.abs(errors)))
228 if not all(np.isfinite(array).all() for array in (values, base_values, predictions)):
229 raise ValueError("Non-finite model outputs or attributions cannot be reported")
231 np.testing.assert_allclose(
232 reconstructed,
233 predictions,
234 rtol=args.rtol,
235 atol=args.atol,
236 )
238 feature_names = (
239 list(explanation.feature_names)
240 if explanation.feature_names is not None
241 else X_explain.columns.tolist()
242 )
243 importance = pd.DataFrame(
244 {
245 "feature": feature_names,
246 "mean_abs_shap": np.abs(values).mean(axis=0),
247 "mean_signed_shap": values.mean(axis=0),
248 }
249 ).sort_values("mean_abs_shap", ascending=False, ignore_index=True)
250 importance.to_csv(args.output_dir / "feature_importance.csv", index=False)
252 local_order = np.argsort(np.abs(values[0]))[::-1]
253 local = pd.DataFrame(
254 {
255 "feature": np.asarray(feature_names)[local_order],
256 "feature_value": X_explain.iloc[0].to_numpy()[local_order],
257 "shap_value": values[0, local_order],
258 }
259 )
260 local.to_csv(args.output_dir / "first_row_contributions.csv", index=False)
262 predictions_frame = pd.DataFrame(
263 {
264 "row_index": X_explain.index.astype(str),
265 "prediction": predictions,
266 "base_value": base_values,
267 "reconstructed_output": reconstructed,
268 "reconstruction_error": errors,
269 }
270 )
271 predictions_frame.to_csv(
272 args.output_dir / "prediction_reconstruction.csv",
273 index=False,
274 )
276 plot_title = f"SHAP for {dataset.target_names[int(model.classes_[args.class_index])]} probability"
277 axis = shap.plots.bar(
278 explanation,
279 max_display=args.max_display,
280 show=False,
281 )
282 save_axis(axis, args.output_dir / "bar.png", title=plot_title)
284 axis = shap.plots.beeswarm(
285 explanation,
286 max_display=args.max_display,
287 show=False,
288 )
289 save_axis(axis, args.output_dir / "beeswarm.png", title=plot_title)
291 axis = shap.plots.waterfall(
292 explanation[0],
293 max_display=args.max_display,
294 show=False,
295 )
296 save_axis(axis, args.output_dir / "waterfall-first-row.png", title=plot_title)
298 top_feature = str(importance.loc[0, "feature"])
299 axis = shap.plots.scatter(
300 explanation[:, top_feature],
301 color=explanation,
302 alpha=0.6,
303 show=False,
304 )
305 save_axis(axis, args.output_dir / "scatter-top-feature.png", title=plot_title)
307 metadata = {
308 "python_version": platform.python_version(),
309 "shap_version": shap.__version__,
310 "numpy_version": np.__version__,
311 "scikit_learn_version": sklearn.__version__,
312 "model": type(model).__name__,
313 "model_score": float(model.score(X_test, y_test)),
314 "seed": args.seed,
315 "background_rows": len(explainer.data),
316 "background_requested_rows": args.background_size,
317 "background_mean_prediction": float(
318 model.predict_proba(background)[:, args.class_index].mean()
319 ),
320 "explained_rows": len(X_explain),
321 "feature_count": X_explain.shape[1],
322 "model_input_dtype": "float32",
323 "all_output_shape": list(np.shape(all_outputs.values)),
324 "selected_output_shape": list(values.shape),
325 "selected_class_index": args.class_index,
326 "selected_class_label": str(model.classes_[args.class_index]),
327 "selected_class_meaning": str(dataset.target_names[int(model.classes_[args.class_index])]),
328 "selected_output_name": output_name(explanation, args.class_index),
329 "output_units": "class probability",
330 "feature_perturbation": "interventional",
331 "max_abs_additivity_error": max_abs_error,
332 "additivity_atol": args.atol,
333 "additivity_rtol": args.rtol,
334 "top_feature": top_feature,
335 "interpretation": "Predictive attribution demo; no causal or clinical validity claim.",
336 }
337 (args.output_dir / "metadata.json").write_text(
338 json.dumps(metadata, indent=2) + "\n",
339 encoding="utf-8",
340 )
342 print(json.dumps(metadata, indent=2))
343 print(f"Report written to {args.output_dir.resolve()}")
344 return 0
347if __name__ == "__main__":
348 raise SystemExit(main())