Skills › Scientific Agent Skills › relsa-severity-assessment › kde_thresholds.py
kde_thresholds.py
388 lines · skills/relsa-severity-assessment/scripts/kde_thresholds.py at 5f2bc03
1#!/usr/bin/env python3
2"""Kernel-density thresholds on the RELSA scale (severity zones).
4Lutscher et al. (2026) define candidate severity thresholds by estimating the
5probability density of the RELSA scores observed in a model and taking the
6*minima* of that density — the sparsely populated valleys between clusters of
7scores. Two minima split the scale into three zones: normal, attention, danger.
9The published sepsis analysis (n = 7 mice, 239 scores) yields minima at
10RELSA = 0.337 and 0.643. This KDE uses the same Gaussian kernel, ``bw.nrd0``
11bandwidth and grid extent as R; direct SciPy evaluation differs from R's
12FFT approximation, so numerical parity is not promised.
14Thresholds are **model-specific and bandwidth-sensitive**, and they are not the
15severity categories of EU Directive 2010/63/EU. Report them as candidate zones
16for one model with the reference set and bandwidth stated, never as a
17regulatory grading.
18"""
20from __future__ import annotations
22import argparse
23import json
24import sys
25from dataclasses import dataclass, field
26from pathlib import Path
27from typing import Iterable, Sequence
29import numpy as np
30import pandas as pd
32DEFAULT_GRID = 512
33DEFAULT_CUT = 3.0
34IQR_TO_SIGMA = 1.349 # R's bw.nrd0 uses this literal, not 1.34898
37def bw_nrd0(values: np.ndarray) -> float:
38 """Silverman's rule of thumb exactly as R's ``bw.nrd0`` computes it.
40 ``0.9 * min(sd, IQR/1.349) * n^(-1/5)``, with R's fallback chain when the
41 spread estimate collapses to zero. scipy's own ``'silverman'`` factor is a
42 *different* formula, so passing it would shift every threshold.
43 """
44 values = np.asarray(values, dtype=float)
45 values = values[np.isfinite(values)]
46 if values.size < 2:
47 raise ValueError("need at least 2 finite values to estimate a bandwidth")
48 sd = float(np.std(values, ddof=1))
49 q75, q25 = np.percentile(values, [75, 25]) # linear interpolation == R type 7
50 lo = min(sd, float(q75 - q25) / IQR_TO_SIGMA)
51 if lo == 0:
52 lo = sd or abs(float(values[0])) or 1.0
53 return 0.9 * lo * values.size ** (-0.2)
56@dataclass
57class ThresholdResult:
58 """Density minima on the RELSA scale plus the zones they imply."""
60 thresholds: list[float]
61 modes: list[float]
62 bandwidth: float
63 n: int
64 grid: np.ndarray = field(repr=False)
65 density: np.ndarray = field(repr=False)
66 zone_counts: dict[str, int] = field(default_factory=dict)
68 def zone_edges(self) -> list[tuple[float, float]]:
69 edges = [0.0, *self.thresholds, float("inf")]
70 return list(zip(edges[:-1], edges[1:]))
72 def zone_names(self) -> list[str]:
73 n_zones = len(self.thresholds) + 1
74 if n_zones == 1:
75 return ["all"]
76 if n_zones == 2:
77 return ["normal", "danger"]
78 if n_zones == 3:
79 return ["normal", "attention", "danger"]
80 return [f"zone{i + 1}" for i in range(n_zones)]
82 def assign(self, values: Iterable[float]) -> list[str]:
83 """Label each RELSA score with the zone it falls in."""
84 names = self.zone_names()
85 out: list[str] = []
86 for value in values:
87 if not np.isfinite(value):
88 out.append("undefined")
89 continue
90 index = int(np.searchsorted(self.thresholds, value, side="right"))
91 out.append(names[index])
92 return out
94 def describe(self) -> str:
95 lines = [
96 f"KDE on {self.n} RELSA scores (bandwidth = {self.bandwidth:.4f})",
97 " candidate thresholds (density minima): "
98 + (", ".join(f"{t:.3f}" for t in self.thresholds) or "none found"),
99 " density modes: " + ", ".join(f"{m:.3f}" for m in self.modes),
100 ]
101 for name, (low, high) in zip(self.zone_names(), self.zone_edges()):
102 span = f"[{low:.3f}, {high:.3f})" if np.isfinite(high) else f">= {low:.3f}"
103 count = self.zone_counts.get(name)
104 tail = f" n={count} ({100 * count / self.n:.1f}%)" if count is not None else ""
105 lines.append(f" {name:<10}{span}{tail}")
106 return "\n".join(lines)
108 def as_dict(self) -> dict:
109 return {
110 "thresholds": [float(t) for t in self.thresholds],
111 "modes": [round(float(m), 4) for m in self.modes],
112 "bandwidth": round(float(self.bandwidth), 6),
113 "n": self.n,
114 "zones": {
115 name: {
116 "low": float(low),
117 "high": None if not np.isfinite(high) else float(high),
118 "n": self.zone_counts.get(name),
119 }
120 for name, (low, high) in zip(self.zone_names(), self.zone_edges())
121 },
122 }
125def density_curve(
126 values: Iterable[float],
127 bandwidth: float | None = None,
128 grid_size: int = DEFAULT_GRID,
129 cut: float = DEFAULT_CUT,
130) -> tuple[np.ndarray, np.ndarray, float]:
131 """Gaussian KDE on R's grid: ``[min - cut*bw, max + cut*bw]``, ``grid_size`` points."""
132 from scipy.stats import gaussian_kde
134 data = np.asarray(list(values), dtype=float)
135 if data.ndim != 1 or np.isinf(data).any():
136 raise ValueError("scores must be a vector of finite values or NaN")
137 data = data[np.isfinite(data)]
138 if data.size < 2:
139 raise ValueError("need at least 2 finite RELSA scores")
141 if np.any(data < 0):
142 raise ValueError("RELSA scores must be nonnegative")
143 bw = float(bandwidth) if bandwidth is not None else bw_nrd0(data)
144 if not np.isfinite(bw) or bw <= 0:
145 raise ValueError("bandwidth must be finite and positive")
146 if int(grid_size) != grid_size or grid_size < 3:
147 raise ValueError("grid_size must be an integer of at least three")
148 if not np.isfinite(cut) or cut <= 0:
149 raise ValueError("cut must be finite and positive")
150 # gaussian_kde scales its factor by the sample sd, so divide it back out.
151 sd = float(np.std(data, ddof=1))
152 grid = np.linspace(data.min() - cut * bw, data.max() + cut * bw, grid_size)
153 if sd == 0:
154 # A repeated constant has one Gaussian mode; SciPy covariance is singular.
155 density = np.exp(-0.5 * ((grid - data[0]) / bw) ** 2) / (bw * np.sqrt(2 * np.pi))
156 return grid, density, bw
157 kde = gaussian_kde(data, bw_method=bw / sd)
158 return grid, kde(grid), bw
161def _zone_counts(data: np.ndarray, thresholds: Sequence[float]) -> list[int]:
162 edges = [-np.inf, *thresholds, np.inf]
163 return [
164 int(np.sum((data >= low) & (data < high)))
165 for low, high in zip(edges[:-1], edges[1:])
166 ]
169def _prune_thin_zones(
170 data: np.ndarray,
171 minima: list[tuple[float, float]],
172 min_zone_fraction: float,
173) -> list[tuple[float, float]]:
174 """Drop thresholds that carve off a zone holding almost no observations.
176 A finite sample's density estimate wiggles in the tails, and a wiggle can
177 produce a minimum that separates one stray observation from the rest. That
178 is a property of the smoother, not a severity zone. Each surviving zone must
179 hold at least ``min_zone_fraction`` of the scores; when one does not, the
180 shallowest threshold bounding it is removed and the check repeats.
181 """
182 if min_zone_fraction <= 0 or not minima:
183 return minima
184 kept = sorted(minima)
185 floor = min_zone_fraction * data.size
186 while kept:
187 counts = _zone_counts(data, [m[0] for m in kept])
188 thinnest = int(np.argmin(counts))
189 if counts[thinnest] >= floor:
190 break
191 # Zone i is bounded by thresholds i-1 and i; drop the shallower one.
192 bounding = [j for j in (thinnest - 1, thinnest) if 0 <= j < len(kept)]
193 kept.pop(max(bounding, key=lambda j: kept[j][1]))
194 return kept
197def find_thresholds(
198 values: Iterable[float],
199 bandwidth: float | None = None,
200 grid_size: int = DEFAULT_GRID,
201 cut: float = DEFAULT_CUT,
202 n_thresholds: int | None = None,
203 within_data: bool = True,
204 min_zone_fraction: float = 0.02,
205) -> ThresholdResult:
206 """Locate density minima and turn them into severity zones.
208 ``n_thresholds`` keeps only the *deepest* k minima (the paper keeps two);
209 ``within_data`` discards minima outside the observed score range, which the
210 padded grid can otherwise produce; ``min_zone_fraction`` discards thresholds
211 that would isolate a near-empty zone (see ``_prune_thin_zones``).
213 When the density is unimodal there are no interior minima and the result
214 carries an empty threshold list — a real answer, meaning this cohort's
215 scores form one cluster and give no data-driven place to cut.
216 """
217 if n_thresholds is not None and (int(n_thresholds) != n_thresholds or n_thresholds < 0):
218 raise ValueError("n_thresholds must be a nonnegative integer")
219 if not np.isfinite(min_zone_fraction) or not 0 <= min_zone_fraction <= 1:
220 raise ValueError("min_zone_fraction must be between zero and one")
221 data = np.asarray(list(values), dtype=float)
222 if data.ndim != 1 or np.isinf(data).any():
223 raise ValueError("scores must be a vector of finite values or NaN")
224 data = data[np.isfinite(data)]
225 grid, dens, bw = density_curve(data, bandwidth, grid_size, cut)
227 lower, upper = float(data.min()), float(data.max())
228 minima: list[tuple[float, float]] = []
229 maxima: list[float] = []
230 for i in range(1, len(grid) - 1):
231 if dens[i] <= dens[i - 1] and dens[i] < dens[i + 1]:
232 if not within_data or lower <= grid[i] <= upper:
233 minima.append((float(grid[i]), float(dens[i])))
234 elif dens[i] >= dens[i - 1] and dens[i] > dens[i + 1]:
235 if not within_data or lower <= grid[i] <= upper:
236 maxima.append(float(grid[i]))
238 minima = _prune_thin_zones(data, minima, min_zone_fraction)
239 if n_thresholds is not None and len(minima) > n_thresholds:
240 minima = sorted(sorted(minima, key=lambda m: m[1])[:n_thresholds])
242 thresholds = [m[0] for m in sorted(minima)]
243 result = ThresholdResult(
244 thresholds=thresholds,
245 modes=maxima,
246 bandwidth=bw,
247 n=int(data.size),
248 grid=grid,
249 density=dens,
250 )
251 labels = result.assign(data)
252 result.zone_counts = {
253 name: int(sum(1 for label in labels if label == name))
254 for name in result.zone_names()
255 }
256 return result
259def plot_thresholds(
260 result: ThresholdResult,
261 values: Iterable[float],
262 path: str | Path,
263 title: str = "RELSA severity zones",
264) -> Path:
265 """Density curve with the minima marked, in the style of the paper's Figure 3."""
266 import matplotlib
268 matplotlib.use("Agg")
269 import matplotlib.pyplot as plt
271 data = np.asarray(list(values), dtype=float)
272 data = data[np.isfinite(data)]
274 fig, ax = plt.subplots(figsize=(5.5, 4.0), constrained_layout=True)
275 ax.plot(result.grid, result.density, color="#333333", lw=1.6)
276 ax.fill_between(result.grid, result.density, color="#333333", alpha=0.08)
277 colors = ["#4c9f70", "#e0a458", "#c1453b", "#7b3f9e"]
278 edges = [result.grid.min(), *result.thresholds, result.grid.max()]
279 for i, (low, high) in enumerate(zip(edges[:-1], edges[1:])):
280 band = (result.grid >= low) & (result.grid <= high)
281 ax.fill_between(
282 result.grid[band], result.density[band], color=colors[i % len(colors)], alpha=0.30
283 )
284 for threshold in result.thresholds:
285 ax.axvline(threshold, ls="--", lw=1.2, color="black")
286 ax.annotate(
287 f"{threshold:.3f}",
288 xy=(threshold, ax.get_ylim()[1]),
289 xytext=(3, -12),
290 textcoords="offset points",
291 fontsize=9,
292 )
293 ax.plot(data, np.full_like(data, -0.02 * result.density.max()), "|", color="#444444",
294 ms=6, alpha=0.6)
295 ax.set_xlabel("RELSA score")
296 ax.set_ylabel("density")
297 ax.set_title(f"{title} (n={result.n}, bw={result.bandwidth:.3f})", fontsize=10)
298 for spine in ("top", "right"):
299 ax.spines[spine].set_visible(False)
300 path = Path(path)
301 fig.savefig(path, dpi=200)
302 plt.close(fig)
303 return path
306def build_parser() -> argparse.ArgumentParser:
307 parser = argparse.ArgumentParser(
308 description="Find candidate RELSA severity thresholds by kernel density estimation.",
309 formatter_class=argparse.ArgumentDefaultsHelpFormatter,
310 )
311 parser.add_argument("scores", help="CSV of RELSA scores (output of relsa_score.py)")
312 parser.add_argument("--column", default="relsa", help="column holding the RELSA score")
313 parser.add_argument(
314 "--group", action="append", default=[], metavar="COL=VALUE",
315 help="restrict to a subset before estimating the density, repeatable",
316 )
317 parser.add_argument(
318 "--n-thresholds", type=int, default=None,
319 help="keep only the k deepest minima (the paper keeps 2)",
320 )
321 parser.add_argument("--bandwidth", type=float, default=None,
322 help="override the bw.nrd0 bandwidth")
323 parser.add_argument("--grid-size", type=int, default=DEFAULT_GRID)
324 parser.add_argument(
325 "--include-outside", action="store_true",
326 help="keep minima outside the observed score range",
327 )
328 parser.add_argument(
329 "--min-zone-fraction", type=float, default=0.02,
330 help="discard a threshold that would isolate a zone holding less than "
331 "this fraction of the scores (tail wiggle in the density estimate)",
332 )
333 parser.add_argument("--plot", help="write a density figure to this path")
334 parser.add_argument("--json", help="write the thresholds to this JSON file")
335 parser.add_argument(
336 "--label-out", help="write the input scores back out with a zone column"
337 )
338 return parser
341def main(argv: Sequence[str] | None = None) -> int:
342 args = build_parser().parse_args(argv)
343 frame = pd.read_csv(args.scores)
344 if args.column not in frame.columns:
345 raise SystemExit(f"column {args.column!r} not in {list(frame.columns)}")
347 for pair in args.group:
348 if "=" not in pair:
349 raise SystemExit(f"--group expects col=value, got {pair!r}")
350 col, value = pair.split("=", 1)
351 if col not in frame.columns:
352 raise SystemExit(f"--group column {col!r} not in data")
353 frame = frame[frame[col].astype(str) == value]
354 if frame.empty:
355 raise SystemExit("no rows left after --group filtering")
357 values = pd.to_numeric(frame[args.column], errors="coerce")
358 result = find_thresholds(
359 values,
360 bandwidth=args.bandwidth,
361 grid_size=args.grid_size,
362 n_thresholds=args.n_thresholds,
363 within_data=not args.include_outside,
364 min_zone_fraction=args.min_zone_fraction,
365 )
366 print(result.describe())
367 print(
368 "\nThese are candidate, model-specific zones on the RELSA scale. They are "
369 "not\nseverity categories under EU Directive 2010/63/EU and are not "
370 "comparable\nacross models or reference sets.",
371 file=sys.stderr,
372 )
374 if args.plot:
375 print(f"wrote {plot_thresholds(result, values, args.plot)}", file=sys.stderr)
376 if args.json:
377 Path(args.json).write_text(json.dumps(result.as_dict(), indent=2) + "\n")
378 print(f"wrote {args.json}", file=sys.stderr)
379 if args.label_out:
380 labelled = frame.copy()
381 labelled["zone"] = result.assign(values)
382 labelled.to_csv(args.label_out, index=False)
383 print(f"wrote {args.label_out}", file=sys.stderr)
384 return 0
387if __name__ == "__main__":
388 raise SystemExit(main())