Skills › Scientific Agent Skills › relsa-severity-assessment › _common.py
_common.py
333 lines · skills/relsa-severity-assessment/scripts/_common.py at 5f2bc03
1"""Shared I/O, validation, and variable-preparation helpers for the RELSA workflow.
3The RELSA long format is one row per animal per time point:
5 id treatment condition time var_1 ... var_n
7``id`` and ``time`` are required; ``treatment`` and ``condition`` are optional
8grouping labels carried through untouched. Time may be days, hours, or minutes
9as long as it increases monotonically per animal. The RELSA convention codes the
10baseline time point as -1, but any value works when it is named explicitly.
11"""
13from __future__ import annotations
15import warnings
16import csv
17from dataclasses import dataclass
18from pathlib import Path
19from typing import Iterable, Sequence
21import numpy as np
22import pandas as pd
24ID_COL = "id"
25TIME_COL = "time"
26META_COLS = ("treatment", "condition")
28# Time-column aliases accepted on input; all are renamed to ``time`` internally.
29TIME_ALIASES = ("time", "day", "days", "hour", "hours", "timepoint", "t")
32class RelsaDataError(ValueError):
33 """Raised when input data cannot support a RELSA calculation."""
36# --------------------------------------------------------------------------- #
37# loading
38# --------------------------------------------------------------------------- #
39def read_relsa_table(
40 path: str | Path,
41 sep: str | None = None,
42 id_col: str = ID_COL,
43 time_col: str | None = None,
44) -> pd.DataFrame:
45 """Read a RELSA-format table from CSV/TSV and return it with canonical names.
47 ``sep=None`` chooses tab for .txt/.tsv/.tab, comma otherwise. An unnamed index column (as written
48 by R's ``write.table``, and present in the published RELSA raw data) is
49 dropped. The time column is detected from ``TIME_ALIASES`` unless named.
50 """
51 path = Path(path)
52 if sep is None:
53 sep = "\t" if path.suffix.lower() in {".txt", ".tsv", ".tab"} else ","
54 with path.open(newline="") as handle:
55 header = next(csv.reader(handle, delimiter=sep))
56 if len(header) != len(set(header)):
57 raise RelsaDataError("duplicate column names in input")
58 frame = pd.read_csv(path, sep=sep, dtype={id_col: "string"},
59 keep_default_na=False, na_values=[""])
61 unnamed = [c for c in frame.columns if str(c).startswith("Unnamed:")]
62 frame = frame.drop(columns=unnamed)
64 return canonicalize(frame, id_col=id_col, time_col=time_col)
67def canonicalize(
68 frame: pd.DataFrame,
69 id_col: str = ID_COL,
70 time_col: str | None = None,
71) -> pd.DataFrame:
72 """Rename the id/time columns to ``id``/``time`` and sort by animal and time."""
73 frame = frame.copy()
74 if not frame.columns.is_unique:
75 raise RelsaDataError("duplicate column names")
77 if id_col != ID_COL:
78 if id_col not in frame.columns:
79 raise RelsaDataError(f"id column {id_col!r} not in {list(frame.columns)}")
80 if ID_COL in frame.columns:
81 raise RelsaDataError("renaming the id column would overwrite 'id'")
82 frame = frame.rename(columns={id_col: ID_COL})
83 if ID_COL not in frame.columns:
84 raise RelsaDataError(f"no {ID_COL!r} column in {list(frame.columns)}")
86 if time_col is None:
87 found = [c for c in frame.columns if str(c).lower() in TIME_ALIASES]
88 if not found:
89 raise RelsaDataError(
90 "no time column found; expected one of "
91 f"{TIME_ALIASES} or an explicit time_col"
92 )
93 if len(found) != 1:
94 raise RelsaDataError("multiple time columns; specify time_col explicitly")
95 time_col = found[0]
96 elif time_col not in frame.columns:
97 raise RelsaDataError(f"time column {time_col!r} not in {list(frame.columns)}")
98 if time_col != TIME_COL:
99 if TIME_COL in frame.columns:
100 raise RelsaDataError("renaming the time column would overwrite 'time'")
101 frame = frame.rename(columns={time_col: TIME_COL})
103 frame[TIME_COL] = pd.to_numeric(frame[TIME_COL], errors="coerce")
104 if not np.isfinite(frame[TIME_COL]).all():
105 raise RelsaDataError("time column must contain finite numeric values")
106 if frame[ID_COL].isna().any() or frame[ID_COL].astype(str).str.strip().eq("").any():
107 raise RelsaDataError("animal IDs must be nonempty")
109 return frame.sort_values([ID_COL, TIME_COL], kind="stable").reset_index(drop=True)
112def variable_columns(frame: pd.DataFrame) -> list[str]:
113 """Every numeric measurement column: not id, time, or a metadata label."""
114 reserved = {ID_COL, TIME_COL, *META_COLS}
115 return [c for c in frame.columns if c not in reserved]
118# --------------------------------------------------------------------------- #
119# validation
120# --------------------------------------------------------------------------- #
121def validate(frame: pd.DataFrame, variables: Sequence[str]) -> None:
122 """Check the invariants RELSA depends on, raising or warning as appropriate.
124 Duplicate (id, time) rows are fatal: normalization and forecasting both
125 assume one measurement per animal per time point. Everything else is a
126 warning, because RELSA is explicitly designed to tolerate missing data.
127 """
128 if not frame.columns.is_unique or not frame.index.is_unique:
129 raise RelsaDataError("column names and row index must be unique")
130 if not variables or len(set(variables)) != len(variables):
131 raise RelsaDataError("variables must be nonempty and unique")
132 if ID_COL not in frame or TIME_COL not in frame:
133 raise RelsaDataError("id and time columns are required")
134 if frame[ID_COL].isna().any() or frame[ID_COL].astype(str).str.strip().eq("").any():
135 raise RelsaDataError("animal IDs must be nonempty")
136 if not np.isfinite(pd.to_numeric(frame[TIME_COL], errors="coerce")).all():
137 raise RelsaDataError("time column must contain finite numeric values")
138 missing = [v for v in variables if v not in frame.columns]
139 if missing:
140 raise RelsaDataError(f"variables not in data: {missing}")
142 dup = frame.duplicated([ID_COL, TIME_COL], keep=False)
143 if dup.any():
144 offenders = (
145 frame.loc[dup, [ID_COL, TIME_COL]].astype(str).agg(" @ ".join, axis=1).unique()
146 )
147 raise RelsaDataError(
148 "multiple rows share an (id, time) pair, so normalization would be "
149 f"ambiguous: {list(offenders)[:5]}. Aggregate per time point first "
150 "(the published models average hourly telemetry to one value per day)."
151 )
153 for var in variables:
154 col = pd.to_numeric(frame[var], errors="coerce")
155 if (frame[var].notna() & col.isna()).any() or np.isinf(col).any():
156 raise RelsaDataError(f"variable {var!r} has nonnumeric or infinite measurements")
157 if col.notna().sum() == 0:
158 warnings.warn(f"variable {var!r} is entirely missing", stacklevel=2)
160 per_animal = frame.groupby(ID_COL)[TIME_COL].size()
161 if per_animal.nunique() > 1:
162 warnings.warn(
163 "animals have different numbers of time points; RELSA handles this, "
164 "but check that gaps are genuinely missing data and not misaligned time axes",
165 stacklevel=2,
166 )
169# --------------------------------------------------------------------------- #
170# variable preparation
171# --------------------------------------------------------------------------- #
172def score_to_percent(
173 values: Iterable[float],
174 max_score: float,
175 baseline_score: float = 0.0,
176) -> np.ndarray:
177 """Map an ordinal severity score onto the RELSA percent scale.
179 RELSA normalizes by dividing each measurement by its own baseline, which is
180 undefined for a clinical score whose healthy baseline is 0. This maps the
181 score's *scale* instead of its ratio: ``baseline_score`` becomes 100 and
182 ``max_score`` becomes 200, so the variable behaves like any other "turned"
183 parameter (rises above 100 as the animal worsens) and one full score point
184 is worth ``100 / |max_score - baseline_score|`` percent.
186 ``max_score`` may lie *below* ``baseline_score`` for scales where a lower
187 number is worse — a nesting score where a well-built nest scores 5 and no
188 nest scores 0 is ``score_to_percent(nesting, max_score=0,
189 baseline_score=5)``. Either way the mapped variable is "turned".
191 Pass the resulting column straight to the RELSA calculation, list it in
192 ``turned``, and leave it out of ``normalize`` — it is already normalized.
193 """
194 values = np.asarray(list(values), dtype=float)
195 span = float(max_score) - float(baseline_score)
196 if not np.isfinite([max_score, baseline_score]).all() or span == 0:
197 raise RelsaDataError(
198 f"max_score ({max_score}) must differ from baseline_score ({baseline_score})"
199 )
200 finite = values[np.isfinite(values)]
201 lo, hi = sorted((float(max_score), float(baseline_score)))
202 if np.isinf(values).any() or np.any((finite < lo) | (finite > hi)):
203 raise RelsaDataError("ordinal score lies outside its declared scale")
204 return 100.0 + 100.0 * (values - float(baseline_score)) / span
207def percent_of_baseline(
208 frame: pd.DataFrame,
209 variables: Sequence[str],
210 baseline_time: float | Sequence[float] | None = None,
211) -> pd.DataFrame:
212 """Express each variable as a percentage of that animal's own baseline.
214 ``baseline_time`` selects the baseline: a single time value, several time
215 values (averaged, i.e. a baseline window), or ``None`` for each animal's
216 first time point. Animals whose baseline is missing or zero yield all-NaN
217 for that variable, with a warning — a zero baseline makes the ratio
218 undefined, which is what ``score_to_percent`` exists to avoid.
219 """
220 validate(frame, variables)
221 out = frame.copy()
222 if baseline_time is None:
223 window: list[float] | None = None
224 elif np.isscalar(baseline_time):
225 window = [float(baseline_time)] # type: ignore[arg-type]
226 else:
227 window = [float(t) for t in baseline_time] # type: ignore[union-attr]
228 if window is not None and (not window or not np.isfinite(window).all()):
229 raise RelsaDataError("baseline window must contain finite times")
231 for var in variables:
232 out[var] = pd.to_numeric(out[var], errors="raise").astype(float)
234 problems: list[str] = []
235 for animal, block in out.groupby(ID_COL, sort=False):
236 if window is None:
237 rows = block.index[block[TIME_COL] == block[TIME_COL].min()]
238 else:
239 rows = block.index[block[TIME_COL].isin(window)]
240 if len(rows) == 0:
241 problems.append(f"{animal} (no baseline time point)")
242 out.loc[block.index, list(variables)] = np.nan
243 continue
244 for var in variables:
245 window_values = out.loc[rows, var].to_numpy(dtype=float)
246 finite = window_values[np.isfinite(window_values)]
247 base = float(finite.mean()) if finite.size else np.nan
248 if not np.isfinite(base) or base <= 0:
249 problems.append(f"{animal}/{var} (baseline {base})")
250 out.loc[block.index, var] = np.nan
251 continue
252 out.loc[block.index, var] = out.loc[block.index, var] / base * 100.0
254 if problems:
255 warnings.warn(
256 "baseline missing or nonpositive, variable set to NaN for: "
257 + ", ".join(problems[:8])
258 + ("..." if len(problems) > 8 else "")
259 + ". For scores whose healthy baseline is 0, use score_to_percent().",
260 stacklevel=2,
261 )
262 return out
265def parse_list(value: str | None) -> list[str]:
266 """Split a comma-separated CLI option into a clean list."""
267 if not value:
268 return []
269 return [item.strip() for item in value.split(",") if item.strip()]
272# --------------------------------------------------------------------------- #
273# metrics
274# --------------------------------------------------------------------------- #
275@dataclass(frozen=True)
276class ForecastMetrics:
277 """The three error metrics reported in Lutscher et al. (2026), Table 1."""
279 n: int
280 rmse: float
281 picp: float
282 mpiw: float
283 n_interval: int = 0
285 def as_dict(self) -> dict[str, float]:
286 return {"n": self.n, "n_interval": self.n_interval,
287 "rmse": self.rmse, "picp": self.picp, "mpiw": self.mpiw}
290def forecast_metrics(
291 actual: Iterable[float],
292 predicted: Iterable[float],
293 lower: Iterable[float] | None = None,
294 upper: Iterable[float] | None = None,
295) -> ForecastMetrics:
296 """RMSE, prediction-interval coverage probability, and mean interval width.
298 RMSE and PICP answer different questions and are reported together on
299 purpose: a model can widen its intervals until PICP hits 100% without
300 predicting anything, which is why MPIW (the mean width, in RELSA units)
301 has to be read alongside the coverage.
302 """
303 a = np.asarray(list(actual), dtype=float)
304 p = np.asarray(list(predicted), dtype=float)
305 if a.ndim != 1 or a.shape != p.shape:
306 raise ValueError(f"actual {a.shape} and predicted {p.shape} differ in length")
308 ok = np.isfinite(a) & np.isfinite(p)
309 rmse = float(np.sqrt(np.mean((a[ok] - p[ok]) ** 2))) if ok.any() else float("nan")
311 picp = float("nan")
312 mpiw = float("nan")
313 n_interval = 0
314 if (lower is None) != (upper is None):
315 raise ValueError("supply both lower and upper intervals")
316 if lower is not None and upper is not None:
317 lo = np.asarray(list(lower), dtype=float)
318 hi = np.asarray(list(upper), dtype=float)
319 if lo.shape != a.shape or hi.shape != a.shape:
320 raise ValueError("interval arrays must match actual and predicted")
321 if np.any(lo > hi):
322 raise ValueError("lower interval exceeds upper interval")
323 band = ok & np.isfinite(lo) & np.isfinite(hi)
324 n_interval = int(band.sum())
325 if band.any():
326 mpiw = float(np.mean(hi[band] - lo[band]))
327 cov = band & np.isfinite(a)
328 if cov.any():
329 inside = (a[cov] >= lo[cov]) & (a[cov] <= hi[cov])
330 picp = float(100.0 * np.mean(inside))
332 return ForecastMetrics(n=int(ok.sum()), rmse=rmse, picp=picp, mpiw=mpiw,
333 n_interval=n_interval)