Skills › Scientific Agent Skills › relsa-severity-assessment › relsa_score.py
relsa_score.py
534 lines · skills/relsa-severity-assessment/scripts/relsa_score.py at 5f2bc03
1#!/usr/bin/env python3
2"""RELSA (RELative Severity Assessment) score — an independent Python implementation.
4Implements the four-step procedure of Talbot et al. (2022), Front. Vet. Sci.
59:937711, as coded in the R package ``mytalbot/RELSA``:
71. directionality — variables that *rise* under duress are "turned";
82. normalization — each variable divided by that animal's own baseline (100%);
93. reference set — the cohort of assumed greatest burden fixes the scale;
104. weights and score — per-variable deviation relative to the reference
11 maximum, combined by a root-mean-square:
13 delta_i(t) = 100 - x_i(t) (turned: x_i(t) - 100), floored at 0
14 RW_i(t) = delta_i(t) / |100 - max_i,ref|
15 RELSA(t) = sqrt( mean_i RW_i(t)^2 ) over non-missing i
17RELSA = 0 means no measured deviation in the declared worsening directions.
18RELSA = 1 is a unit RMS of per-variable reference-scaled deviations, not a
19universal endpoint or necessarily the reference cohort's observed maximum. Missing variables are
20dropped from the mean rather than imputed, so a score is defined whenever at
21least one variable was measured.
23Numerically reproduces the R package's published worked example (surgery
24dataset, animal Ca_001) to the two decimals the package prints; see
25``tests/relsa-severity-assessment/``.
26"""
28from __future__ import annotations
30import argparse
31import json
32import sys
33import warnings
34from dataclasses import dataclass, field
35from pathlib import Path
36from typing import Sequence
38import numpy as np
39import pandas as pd
41sys.path.insert(0, str(Path(__file__).resolve().parent))
43from _common import ( # noqa: E402
44 ID_COL,
45 TIME_COL,
46 RelsaDataError,
47 parse_list,
48 percent_of_baseline,
49 read_relsa_table,
50 score_to_percent,
51 validate,
52 variable_columns,
53)
55BASELINE_PCT = 100.0
58# --------------------------------------------------------------------------- #
59# reference model
60# --------------------------------------------------------------------------- #
61@dataclass
62class ReferenceModel:
63 """The severity scale RELSA scores are relative to.
65 ``maxsev`` is the most extreme normalized value each variable reached
66 anywhere in the reference cohort (minimum, or maximum for turned
67 variables); ``maxdelta`` is its distance from baseline, the denominator of
68 every weight. A RELSA score is only interpretable together with the
69 reference set that produced it, so record ``label`` in any report.
70 """
72 variables: tuple[str, ...]
73 turned: tuple[str, ...]
74 maxsev: dict[str, float]
75 maxdelta: dict[str, float]
76 n_animals: int
77 n_rows: int
78 label: str = ""
79 baseline_time: float | list[float] | None = None
80 notes: dict[str, str] = field(default_factory=dict)
81 preprocessing: dict | None = None
83 def __post_init__(self) -> None:
84 if not self.variables or len(set(self.variables)) != len(self.variables):
85 raise RelsaDataError("reference variables must be nonempty and unique")
86 if set(self.turned) - set(self.variables):
87 raise RelsaDataError("turned variables must belong to the reference")
88 if set(self.maxsev) != set(self.variables) or set(self.maxdelta) != set(self.variables):
89 raise RelsaDataError("reference extrema must cover exactly its variables")
90 for var in self.variables:
91 delta = self.maxdelta[var]
92 if not np.isfinite([self.maxsev[var], delta]).all() or delta <= 0:
93 raise RelsaDataError(f"invalid reference denominator for {var!r}")
94 if not np.isclose(delta, abs(BASELINE_PCT - self.maxsev[var])):
95 raise RelsaDataError(f"inconsistent reference denominator for {var!r}")
97 def to_json(self, path: str | Path) -> None:
98 payload = {
99 "variables": list(self.variables),
100 "turned": list(self.turned),
101 "maxsev": self.maxsev,
102 "maxdelta": self.maxdelta,
103 "n_animals": self.n_animals,
104 "n_rows": self.n_rows,
105 "label": self.label,
106 "baseline_time": self.baseline_time,
107 "notes": self.notes,
108 "preprocessing": self.preprocessing,
109 }
110 Path(path).write_text(json.dumps(payload, indent=2, sort_keys=True) + "\n")
112 @classmethod
113 def from_json(cls, path: str | Path) -> "ReferenceModel":
114 payload = json.loads(Path(path).read_text())
115 return cls(
116 variables=tuple(payload["variables"]),
117 turned=tuple(payload["turned"]),
118 maxsev={k: float(v) for k, v in payload["maxsev"].items()},
119 maxdelta={k: float(v) for k, v in payload["maxdelta"].items()},
120 n_animals=int(payload["n_animals"]),
121 n_rows=int(payload["n_rows"]),
122 label=payload.get("label", ""),
123 baseline_time=payload.get("baseline_time"),
124 notes=payload.get("notes", {}),
125 preprocessing=payload.get("preprocessing"),
126 )
128 def describe(self) -> str:
129 rows = [
130 f"reference model: {self.label or '(unlabelled)'}",
131 f" animals={self.n_animals} rows={self.n_rows} "
132 f"baseline_time={self.baseline_time}",
133 f" {'variable':<12}{'turned':>8}{'max reached':>14}{'max delta':>12}",
134 ]
135 for var in self.variables:
136 rows.append(
137 f" {var:<12}{'yes' if var in self.turned else 'no':>8}"
138 f"{self.maxsev[var]:>14.2f}{self.maxdelta[var]:>12.2f}"
139 )
140 return "\n".join(rows)
143def build_reference(
144 frame: pd.DataFrame,
145 variables: Sequence[str],
146 turned: Sequence[str] = (),
147 baseline_time: float | Sequence[float] | None = None,
148 label: str = "",
149) -> ReferenceModel:
150 """Derive the reference scale from an already-normalized cohort.
152 ``frame`` must be on the percent scale (baseline = 100) — run
153 ``prepare`` or ``_common.percent_of_baseline`` first. ``turned`` names the
154 variables whose *increase* signals worsening.
155 """
156 variables = list(variables)
157 if set(turned) - set(variables):
158 raise RelsaDataError("turned variables must belong to variables")
159 turned = list(turned)
160 validate(frame, variables)
162 maxsev: dict[str, float] = {}
163 maxdelta: dict[str, float] = {}
164 degenerate: list[str] = []
165 for var in variables:
166 col = pd.to_numeric(frame[var], errors="coerce").to_numpy(dtype=float)
167 if not np.isfinite(col).any():
168 raise RelsaDataError(
169 f"variable {var!r} has no finite values in the reference set"
170 )
171 extreme = float(np.nanmax(col)) if var in turned else float(np.nanmin(col))
172 maxsev[var] = extreme
173 maxdelta[var] = abs(BASELINE_PCT - extreme)
174 if maxdelta[var] == 0:
175 degenerate.append(var)
177 if degenerate:
178 raise RelsaDataError(
179 f"variables {degenerate} never deviate from baseline in the reference "
180 "set, so their weight would divide by zero. Drop them, or choose a "
181 "reference cohort that actually expresses the burden (the published "
182 "models use the treatment group of assumed greatest severity)."
183 )
185 wrong_way = [
186 v for v in variables if (maxsev[v] > BASELINE_PCT) != (v in turned)
187 ]
188 if wrong_way:
189 warnings.warn(
190 f"variables {wrong_way} deviate in the opposite direction to their "
191 "declared directionality — check the `turned` list. A variable that "
192 "only ever rises must be turned; one that only falls must not be.",
193 stacklevel=2,
194 )
196 return ReferenceModel(
197 variables=tuple(variables),
198 turned=tuple(turned),
199 maxsev=maxsev,
200 maxdelta=maxdelta,
201 n_animals=int(frame[ID_COL].nunique()),
202 n_rows=int(len(frame)),
203 label=label,
204 baseline_time=(
205 None
206 if baseline_time is None
207 else (float(baseline_time) if np.isscalar(baseline_time) else [float(t) for t in baseline_time]) # type: ignore[arg-type]
208 ),
209 )
212# --------------------------------------------------------------------------- #
213# score
214# --------------------------------------------------------------------------- #
215def relsa_weights(
216 frame: pd.DataFrame,
217 reference: ReferenceModel,
218 drop: Sequence[str] = (),
219 round_digits: int | None = 2,
220) -> pd.DataFrame:
221 """Per-variable RELSA weights (RW) for every animal and time point.
223 Rounding to two decimals is the R package's behaviour and is applied to the
224 deltas and to the weights, before the root-mean-square. Pass
225 ``round_digits=None`` for full precision (scores then differ from R in the
226 third decimal).
227 """
228 if set(drop) - set(reference.variables):
229 raise RelsaDataError("drop contains unknown reference variables")
230 used = [v for v in reference.variables if v not in set(drop)]
231 if not used:
232 raise RelsaDataError("every reference variable was dropped")
233 missing = [v for v in used if v not in frame.columns]
234 if missing:
235 raise RelsaDataError(
236 f"variables {missing} are in the reference model but not in the data. "
237 "The test set must carry the same variable names as the reference."
238 )
239 validate(frame, used)
241 out = frame[[ID_COL, TIME_COL]].copy()
242 for var in used:
243 values = pd.to_numeric(frame[var], errors="coerce").to_numpy(dtype=float)
244 delta = (
245 values - BASELINE_PCT if var in reference.turned else BASELINE_PCT - values
246 )
247 delta = np.where(np.isfinite(delta), np.clip(delta, 0.0, None), np.nan)
248 if round_digits is not None:
249 delta = np.round(delta, round_digits)
250 weight = delta / reference.maxdelta[var]
251 if round_digits is not None:
252 weight = np.round(weight, round_digits)
253 out[var] = weight
254 return out
257def _warn_on_composition_change(scored: pd.DataFrame, used: Sequence[str]) -> None:
258 """Warn when the set of available variables changes along an animal's trajectory.
260 RELSA averages over whichever variables were measured, so a variable that
261 appears or disappears mid-trajectory moves the score even when the animal's
262 state is unchanged. In the published sepsis data, body weight is recorded
263 only on the day of euthanasia; including it drops that animal's endpoint
264 score from 0.93 to 0.83 purely by changing the mean's composition. Score the
265 variables that are measured throughout, and keep the rest as separate
266 endpoint criteria.
267 """
268 offenders: list[str] = []
269 for animal, block in scored.groupby(ID_COL, sort=False):
270 present = block[list(used)].notna()
271 signatures = {tuple(row) for row in present.itertuples(index=False)}
272 signatures.discard(tuple([False] * len(used))) # fully missing rows are fine
273 if len(signatures) > 1:
274 varying = [
275 var for var in used if present[var].nunique(dropna=False) > 1
276 ]
277 offenders.append(f"{animal}: {', '.join(varying)}")
278 if offenders:
279 warnings.warn(
280 "the set of measured variables changes across time points, so RELSA "
281 "scores along these trajectories are not strictly comparable — a "
282 "variable appearing or disappearing shifts the score by itself: "
283 + "; ".join(offenders[:5])
284 + ("..." if len(offenders) > 5 else ""),
285 stacklevel=3,
286 )
289def relsa_scores(
290 frame: pd.DataFrame,
291 reference: ReferenceModel,
292 drop: Sequence[str] = (),
293 round_digits: int | None = 2,
294 keep_meta: bool = True,
295) -> pd.DataFrame:
296 """RELSA score per animal per time point, with the weights that produced it.
298 Returns ``id``, ``time``, any metadata columns, one column per variable
299 holding its weight, ``n_vars`` (variables available at that time point), and
300 ``relsa``. Rows where every variable is missing get ``relsa = NaN``.
301 """
302 weights = relsa_weights(frame, reference, drop=drop, round_digits=round_digits)
303 used = [c for c in weights.columns if c not in (ID_COL, TIME_COL)]
305 matrix = weights[used].to_numpy(dtype=float)
306 available = np.isfinite(matrix).sum(axis=1)
307 with np.errstate(invalid="ignore", divide="ignore"):
308 mean_square = np.nansum(matrix**2, axis=1) / available
309 score = np.sqrt(mean_square)
310 score = np.where(available == 0, np.nan, score)
311 if round_digits is not None:
312 score = np.round(score, round_digits)
314 out = weights.copy()
315 out["n_vars"] = available
316 out["relsa"] = score
318 _warn_on_composition_change(out, used)
320 if keep_meta:
321 meta = [c for c in ("treatment", "condition") if c in frame.columns]
322 if meta:
323 labels = frame[[ID_COL, TIME_COL, *meta]]
324 out = out.merge(labels, on=[ID_COL, TIME_COL], how="left")
325 order = [ID_COL, TIME_COL, *meta, *used, "n_vars", "relsa"]
326 out = out[order]
327 return out
330def prepare(
331 frame: pd.DataFrame,
332 normalize: Sequence[str] = (),
333 baseline_time: float | Sequence[float] | None = None,
334) -> pd.DataFrame:
335 """Normalize the named variables to each animal's own baseline (= 100%).
337 Variables *not* named are left untouched — that is the RELSA convention and
338 it matters: body-weight change (bwc) and scores mapped with
339 ``score_to_percent`` already sit on the percent scale, and normalizing baseline-100 values again is redundant and can change the
340 anchor if a different baseline window is selected.
341 """
342 if not normalize:
343 return frame.copy()
344 return percent_of_baseline(frame, list(normalize), baseline_time=baseline_time)
347# --------------------------------------------------------------------------- #
348# CLI
349# --------------------------------------------------------------------------- #
350def _filter_group(frame: pd.DataFrame, pairs: Sequence[str]) -> pd.DataFrame:
351 for pair in pairs:
352 if "=" not in pair:
353 raise SystemExit(f"--reference-group expects col=value, got {pair!r}")
354 col, value = pair.split("=", 1)
355 if col not in frame.columns:
356 raise SystemExit(f"--reference-group column {col!r} not in data")
357 frame = frame[frame[col].astype(str) == value]
358 if frame.empty:
359 raise SystemExit(f"--reference-group {pair!r} selected no rows")
360 return frame
363def _parse_score_scale(spec: str, frame: pd.DataFrame) -> tuple[str, tuple[float, float]]:
364 """Parse ``COL=MAX`` or ``COL=MAX:BASELINE`` from --score-scale."""
365 if "=" not in spec:
366 raise SystemExit(f"--score-scale expects COL=MAX[:BASELINE], got {spec!r}")
367 column, bounds = spec.split("=", 1)
368 if column not in frame.columns:
369 raise SystemExit(f"--score-scale column {column!r} not in data")
370 parts = bounds.split(":")
371 if len(parts) > 2:
372 raise SystemExit("--score-scale accepts exactly MAX or MAX:BASELINE")
373 try:
374 max_score = float(parts[0])
375 baseline = float(parts[1]) if len(parts) > 1 else 0.0
376 except ValueError:
377 raise SystemExit(f"--score-scale bounds must be numeric, got {bounds!r}") from None
378 return column, (max_score, baseline)
381def build_parser() -> argparse.ArgumentParser:
382 parser = argparse.ArgumentParser(
383 description="Compute RELSA severity scores for a RELSA-format table.",
384 formatter_class=argparse.ArgumentDefaultsHelpFormatter,
385 )
386 parser.add_argument("data", help="CSV/TSV in RELSA long format (id, time, variables)")
387 parser.add_argument(
388 "--variables",
389 help="comma-separated variables to score (default: every measurement column)",
390 )
391 parser.add_argument(
392 "--turned",
393 help="comma-separated variables that RISE under worsening severity "
394 "(clinical scores, biomarkers, fever, tachycardia)",
395 )
396 parser.add_argument(
397 "--normalize",
398 help="comma-separated variables to divide by their own baseline; omit any "
399 "variable already on a percent scale (bwc, mapped scores)",
400 )
401 parser.add_argument(
402 "--score-scale",
403 action="append",
404 default=[],
405 metavar="COL=MAX[:BASELINE]",
406 help="map an ordinal severity score onto the percent scale, where MAX is "
407 "the worst possible score and BASELINE the healthy one (default 0). "
408 "Required for scores whose baseline is 0, which cannot be ratio-normalized. "
409 "Implies the variable is turned and not normalized. Repeatable.",
410 )
411 parser.add_argument(
412 "--baseline-time",
413 help="time value(s) defining the baseline, comma-separated to average a "
414 "window (default: each animal's first time point)",
415 )
416 parser.add_argument(
417 "--reference-group",
418 action="append",
419 default=[],
420 metavar="COL=VALUE",
421 help="restrict the reference set, repeatable (e.g. treatment=CLP)",
422 )
423 parser.add_argument(
424 "--reference-data",
425 help="separate table for the reference set (default: the same file)",
426 )
427 parser.add_argument(
428 "--load-reference", help="reuse a reference model saved by --save-reference"
429 )
430 parser.add_argument("--save-reference", help="write the reference model to JSON")
431 parser.add_argument("--drop", help="comma-separated variables to exclude from scoring")
432 parser.add_argument(
433 "--full-precision",
434 action="store_true",
435 help="skip the 2-decimal rounding the R package applies",
436 )
437 parser.add_argument("--out", help="write scores to this CSV (default: stdout)")
438 parser.add_argument("--id-col", default=ID_COL)
439 parser.add_argument("--time-col", help="name of the time column, if not auto-detected")
440 return parser
443def main(argv: Sequence[str] | None = None) -> int:
444 args = build_parser().parse_args(argv)
445 digits = None if args.full_precision else 2
447 reference = ReferenceModel.from_json(args.load_reference) if args.load_reference else None
448 if reference is not None and (args.reference_data or args.reference_group):
449 raise SystemExit("--load-reference cannot be combined with reference selection")
450 saved = reference.preprocessing if reference is not None else None
451 if saved:
452 supplied = {
453 "normalize": parse_list(args.normalize) if args.normalize is not None else None,
454 "score_scale": args.score_scale or None,
455 "baseline_time": args.baseline_time,
456 }
457 for key, value in supplied.items():
458 if value is not None and value != saved[key]:
459 raise SystemExit(f"{key} conflicts with the saved preprocessing contract")
460 if args.full_precision and saved["round_digits"] is not None:
461 raise SystemExit("--full-precision conflicts with the saved reference")
462 normalize = saved["normalize"]
463 score_specs = saved["score_scale"]
464 baseline_spec = saved["baseline_time"]
465 digits = saved["round_digits"]
466 else:
467 normalize = parse_list(args.normalize)
468 score_specs = args.score_scale
469 baseline_spec = args.baseline_time
470 if reference is not None:
471 warnings.warn("legacy reference has no preprocessing contract; explicitly repeat "
472 "the original normalization and score mappings", stacklevel=2)
473 frame = read_relsa_table(args.data, id_col=args.id_col, time_col=args.time_col)
474 variables = parse_list(args.variables) or (list(reference.variables) if reference else variable_columns(frame))
475 if reference and tuple(variables) != reference.variables:
476 raise SystemExit("--variables conflicts with the loaded reference")
477 turned = parse_list(args.turned) if args.turned is not None else (list(reference.turned) if reference else [])
478 baseline_time: float | list[float] | None = None
479 if baseline_spec:
480 times = [float(t) for t in parse_list(baseline_spec)]
481 baseline_time = times[0] if len(times) == 1 else times
483 def transform(raw: pd.DataFrame) -> pd.DataFrame:
484 result = raw.copy()
485 for spec in score_specs:
486 column, bounds = _parse_score_scale(spec, result)
487 result[column] = score_to_percent(result[column], *bounds)
488 return prepare(result, normalize=normalize, baseline_time=baseline_time)
490 for spec in score_specs:
491 column, _ = _parse_score_scale(spec, frame)
492 if column not in turned:
493 turned.append(column)
494 if column in normalize:
495 normalize = [v for v in normalize if v != column]
496 if column not in variables:
497 variables.append(column)
498 if reference and set(turned) != set(reference.turned):
499 raise SystemExit("--turned conflicts with the loaded reference")
500 if set(normalize) - set(variables):
501 raise SystemExit("--normalize contains unscored variables")
502 validate(frame, variables)
503 prepared = transform(frame)
504 if reference is None:
505 ref_frame = transform(read_relsa_table(args.reference_data, id_col=args.id_col,
506 time_col=args.time_col)) if args.reference_data else prepared
507 ref_frame = _filter_group(ref_frame, args.reference_group)
508 label = args.reference_data or args.data
509 if args.reference_group:
510 label = f"{label} [{', '.join(args.reference_group)}]"
511 reference = build_reference(ref_frame, variables=variables, turned=turned,
512 baseline_time=baseline_time, label=label)
513 reference.preprocessing = {
514 "normalize": normalize, "score_scale": score_specs,
515 "baseline_time": baseline_spec, "round_digits": digits,
516 }
517 if args.save_reference:
518 reference.to_json(args.save_reference)
520 print(reference.describe(), file=sys.stderr)
522 scores = relsa_scores(
523 prepared, reference, drop=parse_list(args.drop), round_digits=digits
524 )
525 if args.out:
526 scores.to_csv(args.out, index=False)
527 print(f"wrote {len(scores)} rows to {args.out}", file=sys.stderr)
528 else:
529 scores.to_csv(sys.stdout, index=False)
530 return 0
533if __name__ == "__main__":
534 raise SystemExit(main())