Skills › Scientific Agent Skills › clinical-decision-support › model_biomarker_evaluation.py
model_biomarker_evaluation.py
362 lines · skills/clinical-decision-support/scripts/model_biomarker_evaluation.py at 5f2bc03
1#!/usr/bin/env python3
2"""Create a bounded report from aggregate model/biomarker validation counts."""
4from __future__ import annotations
6import argparse
7import math
8import sys
9from typing import Any
11from _common import (
12 InputError,
13 IssueLog,
14 finite_number,
15 load_json_object,
16 nonnegative_int,
17 print_report,
18 require_list,
19 require_nonempty_text,
20 source_ids,
21 write_json,
22)
24MAX_GROUPS = 12
25MAX_BINS = 30
26MAX_GROUP_N = 10_000_000
29def _wilson(successes: int, total: int) -> dict[str, Any]:
30 if total <= 0:
31 return {"estimate": None, "ci95": [None, None]}
32 z = 1.959963984540054
33 proportion = successes / total
34 denominator = 1.0 + (z * z / total)
35 center = (proportion + z * z / (2.0 * total)) / denominator
36 margin = (
37 z
38 * math.sqrt(
39 proportion * (1.0 - proportion) / total
40 + z * z / (4.0 * total * total)
41 )
42 / denominator
43 )
44 return {
45 "estimate": round(proportion, 6),
46 "ci95": [round(max(0.0, center - margin), 6), round(min(1.0, center + margin), 6)],
47 }
50def _small_nonzero(values: list[int], minimum: int) -> bool:
51 return any(0 < value < minimum for value in values)
54def _metric(successes: int, total: int) -> dict[str, Any]:
55 return _wilson(successes, total)
58def _calibration_summary(
59 bins: list[Any], group_name: str, minimum: int, expected_events: int,
60 expected_n: int,
61) -> tuple[dict[str, Any], bool]:
62 total = 0
63 events = 0
64 weighted_predicted = 0.0
65 weighted_absolute_gap = 0.0
66 sensitive = False
67 for index, entry in enumerate(bins):
68 field = f"groups[{group_name}].calibration_bins[{index}]"
69 if not isinstance(entry, dict):
70 raise InputError(f"{field} must be an object")
71 count = nonnegative_int(entry.get("n"), f"{field}.n")
72 observed = nonnegative_int(
73 entry.get("observed_events"), f"{field}.observed_events"
74 )
75 predicted = finite_number(
76 entry.get("mean_predicted_probability"),
77 f"{field}.mean_predicted_probability",
78 )
79 if count == 0 or observed > count:
80 raise InputError(f"{field} has inconsistent counts")
81 if not 0.0 <= predicted <= 1.0:
82 raise InputError(f"{field}.mean_predicted_probability must be in [0, 1]")
83 observed_rate = observed / count
84 total += count
85 events += observed
86 weighted_predicted += count * predicted
87 weighted_absolute_gap += count * abs(predicted - observed_rate)
88 sensitive = sensitive or _small_nonzero(
89 [count, observed, count - observed], minimum
90 )
91 if events != expected_events:
92 raise InputError(
93 f"Calibration observed events for {group_name} do not match confusion counts"
94 )
95 if total == 0:
96 raise InputError(f"Calibration bins for {group_name} are empty")
97 if total != expected_n:
98 raise InputError(f"Calibration bins for {group_name} must sum to group n")
99 if sensitive:
100 return (
101 {
102 "suppressed": True,
103 "reason": "One or more calibration cells are below the disclosure threshold",
104 },
105 True,
106 )
107 observed_rate = events / total
108 mean_predicted = weighted_predicted / total
109 return (
110 {
111 "suppressed": False,
112 "n": total,
113 "observed_rate": round(observed_rate, 6),
114 "mean_predicted_probability": round(mean_predicted, 6),
115 "calibration_in_the_large_gap": round(mean_predicted - observed_rate, 6),
116 "weighted_absolute_calibration_gap": round(
117 weighted_absolute_gap / total, 6
118 ),
119 "limitations": (
120 "Aggregate-bin calibration is approximate and cannot estimate a "
121 "calibration slope or replace individual-level validation."
122 ),
123 },
124 False,
125 )
128def evaluate(
129 document: dict[str, Any], minimum: int
130) -> tuple[IssueLog, dict[str, Any]]:
131 log = IssueLog()
132 report: dict[str, Any] = {
133 "report_type": "aggregate_model_biomarker_evaluation",
134 "individual_output_generated": False,
135 "clinical_classification_generated": False,
136 "recommendation_generated": False,
137 "groups": [],
138 "subgroup_performance_differences": {},
139 }
140 metric_values: dict[str, list[tuple[str, float]]] = {}
142 try:
143 require_nonempty_text(document.get("schema_version"), "schema_version")
144 known_sources = source_ids(document)
145 metadata = document.get("metadata")
146 if not isinstance(metadata, dict):
147 raise InputError("metadata must be an object")
148 for field in (
149 "evaluation_id",
150 "title",
151 "evaluation_target",
152 "target_version",
153 "purpose",
154 "population",
155 "setting",
156 "outcome",
157 "outcome_horizon",
158 "threshold",
159 "threshold_source_id",
160 "validation_dataset",
161 "data_cut_date",
162 ):
163 require_nonempty_text(metadata.get(field), f"metadata.{field}")
164 if not isinstance(metadata.get("data_level"), str) or metadata.get("data_level") not in {"aggregate", "synthetic"}:
165 log.errors.append("metadata.data_level must be aggregate or synthetic")
166 if metadata.get("person_level_output") is not False:
167 log.errors.append("metadata.person_level_output must be false")
168 if metadata.get("threshold_pre_specified") is not True:
169 log.errors.append("metadata.threshold_pre_specified must be true")
170 if metadata.get("threshold_source_id") not in known_sources:
171 log.errors.append("metadata.threshold_source_id must reference sources")
172 if metadata.get("external_validation") is not True:
173 log.warnings.append("Independent external validation is not documented")
174 if metadata.get("human_review_required") is not True:
175 log.errors.append("metadata.human_review_required must be true")
177 review = document.get("human_review")
178 if not isinstance(review, dict):
179 raise InputError("human_review must be an object")
180 if review.get("required") is not True:
181 log.errors.append("human_review.required must be true")
182 roles = require_list(review.get("roles"), "human_review.roles", maximum=20)
183 if not roles:
184 log.errors.append("At least one human-review role is required")
185 for index, role in enumerate(roles):
186 require_nonempty_text(role, f"human_review.roles[{index}]")
187 if review.get("completed") is not True:
188 log.warnings.append("Human review is not recorded as complete")
190 governance = document.get("governance")
191 if not isinstance(governance, dict):
192 raise InputError("governance must be an object")
193 for field in (
194 "owner",
195 "version_and_change_control",
196 "monitoring",
197 "auditability",
198 "human_factors",
199 "rollback_and_retirement",
200 ):
201 require_nonempty_text(governance.get(field), f"governance.{field}")
203 if not log.ok:
204 raise InputError("Evaluation stopped before processing aggregate counts")
206 groups = require_list(document.get("groups"), "groups", maximum=MAX_GROUPS)
207 if not groups:
208 raise InputError("At least one aggregate group is required")
209 names: set[str] = set()
210 for index, group in enumerate(groups):
211 field = f"groups[{index}]"
212 if not isinstance(group, dict):
213 raise InputError(f"{field} must be an object")
214 name = require_nonempty_text(group.get("name"), f"{field}.name", max_length=100)
215 if name in names:
216 raise InputError(f"Duplicate group name: {name}")
217 names.add(name)
218 n = nonnegative_int(group.get("n"), f"{field}.n")
219 if n == 0 or n > MAX_GROUP_N:
220 raise InputError(f"{field}.n must be between 1 and {MAX_GROUP_N}")
221 confusion = group.get("confusion")
222 if not isinstance(confusion, dict):
223 raise InputError(f"{field}.confusion must be an object")
224 tp = nonnegative_int(confusion.get("tp"), f"{field}.confusion.tp")
225 fp = nonnegative_int(confusion.get("fp"), f"{field}.confusion.fp")
226 tn = nonnegative_int(confusion.get("tn"), f"{field}.confusion.tn")
227 fn = nonnegative_int(confusion.get("fn"), f"{field}.confusion.fn")
228 cells = [tp, fp, tn, fn]
229 if sum(cells) != n:
230 raise InputError(f"{field}.confusion counts must sum to n")
232 group_result: dict[str, Any] = {"name": name}
233 if n < minimum or _small_nonzero(cells, minimum):
234 group_result.update(
235 {
236 "suppressed": True,
237 "reason": (
238 "Group or contributing confusion cell is below the "
239 "disclosure threshold"
240 ),
241 }
242 )
243 else:
244 metrics = {
245 "sensitivity": _metric(tp, tp + fn),
246 "specificity": _metric(tn, tn + fp),
247 "positive_predictive_value": _metric(tp, tp + fp),
248 "negative_predictive_value": _metric(tn, tn + fn),
249 "accuracy": _metric(tp + tn, n),
250 "prevalence": _metric(tp + fn, n),
251 }
252 sensitivity = metrics["sensitivity"]["estimate"]
253 specificity = metrics["specificity"]["estimate"]
254 metrics["balanced_accuracy"] = {
255 "estimate": (
256 round((sensitivity + specificity) / 2.0, 6)
257 if sensitivity is not None and specificity is not None
258 else None
259 ),
260 "ci95": None,
261 }
262 if sensitivity is None or specificity is None:
263 log.warnings.append(
264 f"Group {name} lacks one observed outcome class; balanced "
265 "accuracy and the corresponding class metric are undefined"
266 )
267 group_result.update(
268 {"suppressed": False, "n": n, "metrics": metrics}
269 )
270 for metric_name, metric in metrics.items():
271 estimate = metric.get("estimate")
272 if estimate is not None:
273 metric_values.setdefault(metric_name, []).append(
274 (name, float(estimate))
275 )
277 bins = require_list(
278 group.get("calibration_bins"),
279 f"{field}.calibration_bins",
280 maximum=MAX_BINS,
281 )
282 calibration, _ = _calibration_summary(
283 bins, name, minimum, expected_events=tp + fn, expected_n=n
284 )
285 group_result["calibration"] = calibration
286 report["groups"].append(group_result)
288 for metric_name, values in metric_values.items():
289 if len(values) < 2:
290 continue
291 estimates = [value for _, value in values]
292 low_name, low = min(values, key=lambda item: item[1])
293 high_name, high = max(values, key=lambda item: item[1])
294 report["subgroup_performance_differences"][metric_name] = {
295 "maximum_absolute_difference": round(high - low, 6),
296 "lowest_group": low_name,
297 "highest_group": high_name,
298 "interpretation": (
299 "Descriptive difference only; it is not a fairness judgment "
300 "and may reflect case mix, measurement, sampling, or model behavior."
301 ),
302 }
303 if len(groups) < 2:
304 log.warnings.append("Only one group supplied; subgroup comparisons are unavailable")
305 if not report["subgroup_performance_differences"]:
306 log.warnings.append(
307 "Subgroup differences were not estimable after disclosure suppression"
308 )
309 except InputError as exc:
310 log.errors.append(str(exc))
312 if not log.ok:
313 report["groups"] = []
314 report["subgroup_performance_differences"] = {}
316 report["disclosure_threshold"] = minimum
317 report["limitations"] = [
318 "Input is aggregate and cannot support person-level output.",
319 "Metrics do not establish clinical validity, clinical utility, safety, or fairness.",
320 "Wilson intervals do not account for clustering, repeated observations, censoring, or verification bias.",
321 "Calibration bins are lossy summaries and cannot estimate calibration slope.",
322 ]
323 if log.ok:
324 log.info.append("Aggregate evaluation completed without person-level output")
325 return log, report
328def main() -> int:
329 parser = argparse.ArgumentParser(
330 description=(
331 "Report bounded descriptive performance from local aggregate counts. "
332 "Never emits a person-level class or clinical recommendation."
333 )
334 )
335 parser.add_argument("input", help="Local aggregate JSON")
336 parser.add_argument("-o", "--output", help="Optional local JSON report")
337 parser.add_argument(
338 "--min-cell-size",
339 type=int,
340 default=11,
341 help="Suppress contributing nonzero cells below this value (default: 11)",
342 )
343 args = parser.parse_args()
344 if not 2 <= args.min_cell_size <= 1_000:
345 parser.error("--min-cell-size must be between 2 and 1000")
347 try:
348 document = load_json_object(args.input)
349 log, details = evaluate(document, args.min_cell_size)
350 result = log.as_dict()
351 result["evaluation"] = details
352 if args.output:
353 write_json(args.output, result)
354 print_report(result)
355 except InputError as exc:
356 print_report(IssueLog(errors=[str(exc)]).as_dict())
357 return 2
358 return 0 if result["status"] == "pass" else 1
361if __name__ == "__main__":
362 sys.exit(main())