relsa-severity-assessment skill

RELSA severity assessment and humane endpoint forecasting

Overview

Severity assessment in animal research is legally mandatory and scientifically load-bearing: it drives humane endpoint decisions, and poor welfare monitoring degrades reproducibility. The usual practice evaluates each readout in isolation — weight loss here, a clinical score there — which makes it hard to say how badly an individual animal is actually doing.

This skill implements two published procedures that address that:

per time point, expressed relative to a reference set of known burden. RELSA = 0 is baseline; RELSA = 1 means the animal has reached the reference set's maximum deviation.

• RELSA (Talbot et al., 2022) combines several outcome measures into one score per animal

trajectory and forecasts the next score with a 95% prediction interval, so animals heading for a humane endpoint can be identified before they get there. Kernel density estimation on the RELSA scale supplies candidate attention and danger zones for interpretation.

• foRcast (Lutscher et al., 2026) fits an ARIMA model to an individual animal's RELSA

The point is refinement: give at-risk animals attention earlier, and avoid euthanising animals that would have recovered. Both procedures are aids to severity assessment, not decision rules — see Boundaries.

When to use this skill

per-animal severity score

• Combining weight loss, temperature, clinical scoring, biomarkers, or telemetry into a single

the severity score at a coming time point

• Asking which animals in a cohort are at risk of reaching a humane endpoint, or predicting

relative scale

• Comparing severity between treatment groups, interventions, or animal models on a common

analysis, or an application under EU Directive 2010/63/EU

• Defining thresholds or zones on a severity scale from the data
• Writing the severity-assessment section of an animal welfare report, a 3Rs/refinement

For general forecasting of a time series that is not a severity score, use timesfm-forecasting or statsmodels. For study design and sample size, use experimental-design and statistical-power.

Installation

uv pip install "numpy>=1.26" "pandas>=2.0" "scipy>=1.11" "statsmodels>=0.14" matplotlib

relsascore.py and kdethresholds.py need only numpy/pandas/scipy; statsmodels is required for forecasting and matplotlib only for figures.

Data format

One row per animal per time point, in a CSV:

are optional labels used for grouping and for selecting the reference set.

• id and a time column (day, time, hour, …) are required; treatment and condition

convention codes the baseline time point as -1.

• Time may be days, hours, or minutes — just keep it monotonic per animal. The RELSA

first (the published models average heart rate, HRV, and temperature, and sum activity).

• One row per animal per time point. Average hourly telemetry to one value per interval

missing value treated as "no deviation" biases severity downward.

• Leave missing measurements empty. They are dropped from the score, never imputed — a

assets/example_cohort.csv is a small synthetic cohort (6 mice, 9 days, temperature, body weight, an 0–8 clinical score, and an IL-6-like biomarker) used by every command below, so each one is runnable as written.

The four decisions that determine the result

Make these explicitly and write them into the methods. Nothing else about the procedure matters as much.

1. Directionality — which variables rise under worsening? Falling is the default (body weight, activity, food intake, burrowing, wheel running). Variables that rise must be declared as --turned: clinical scores, inflammatory biomarkers, fever, tachycardia. Get this wrong and the variable contributes nothing at all, silently, because deviations in the "wrong" direction are floored at zero. Body temperature is model-dependent — it falls in sepsis and endotoxaemia, rises in fever models. Nothing in the data can settle this for you: in the published sepsis model activity legitimately swings further above baseline than below, so only a variable that never once moves the declared way is detectable, and build_reference() warns about exactly that case.

2. The reference set — relative to what? RELSA scores mean nothing without it. Use the group assumed to carry the greatest burden in your model (the published studies use the highest-dose or endpoint-reaching treatment group). Too mild a reference pushes every score above 1; too severe compresses everything toward 0. Save it with --save-reference and reuse it with --load-reference so later cohorts stay on the same scale.

3. Scores with a zero baseline. A clinical score of 0 in a healthy animal cannot be ratio-normalized — 0/0 is undefined. Use --score-scale score=8 to map the score's scale instead (healthy → 100%, worst possible → 200%), which also marks it as turned. This mapping is a modelling choice about how much one score point is worth relative to one percent of body weight; state it. The alternative is to keep the score out of RELSA and use it as an independent endpoint criterion.

4. Which variables are measured throughout. Because the score averages over whichever variables are available, a variable that appears or disappears mid-trajectory moves the score by itself. In the published sepsis data, adding body weight — recorded only on the day of euthanasia — drops that animal's endpoint score from 0.93 to 0.83 for no biological reason. relsa_scores() warns when composition changes; score the variables present throughout.

Workflow

Step 1 — compute RELSA scores

python scripts/relsa_score.py assets/example_cohort.csv \
    --variables weight,temp,score,il6 \
    --normalize weight,temp,il6 \
    --turned il6 \
    --score-scale score=8 \
    --baseline-time -1 \
    --reference-group condition=endpoint \
    --save-reference reference.json \
    --out relsa_scores.csv

The reference model is echoed so the scale is auditable:

reference model: assets/example_cohort.csv [condition=endpoint]
  animals=2  rows=18  baseline_time=-1.0
  variable      turned   max reached   max delta
  weight            no         82.40       17.60
  temp              no         92.79        7.21
  score            yes        187.50       87.50
  il6              yes        797.72      697.72

relsa_scores.csv holds each variable's weight alongside the score, which is what makes a score explainable — here M01 deteriorating to its endpoint, M03 peaking on day 3 and recovering:

id  time  weight  temp  score  il6  n_vars  relsa
M01     1    0.46  0.49   0.57 0.52       4   0.51
M01     3    0.84  0.76   1.00 0.89       4   0.88
M01     5    1.00  1.00   1.00 1.00       4   1.00
M03     3    0.56  0.44   0.57 0.54       4   0.53
M03     5    0.35  0.26   0.43 0.32       4   0.35
M03     7    0.12  0.06   0.14 0.11       4   0.11

A weight of 1.00 means that variable hit the reference maximum; n_vars is how many variables entered the score at that time point.

Same thing from Python, when you need the objects:

import sys; sys.path.insert(0, "scripts")
from _common import read_relsa_table, score_to_percent
from relsa_score import prepare, build_reference, relsa_scores

frame = read_relsa_table("assets/example_cohort.csv")
frame["score"] = score_to_percent(frame["score"], max_score=8)   # 0-8 clinical score
VARS, TURNED = ["weight", "temp", "score", "il6"], ["score", "il6"]

prepared  = prepare(frame, normalize=["weight", "temp", "il6"], baseline_time=-1)
reference = build_reference(prepared[prepared.condition == "endpoint"],
                           variables=VARS, turned=TURNED, baseline_time=-1,
                           label="endpoint-reaching animals")
scores    = relsa_scores(prepared, reference)

Step 2 — forecast the endpoint

Train on everything up to the time point before the endpoint, predict the score at the endpoint, and score the prediction:

python scripts/forecast_relsa.py relsa_scores.csv \
    --animals M01,M02 --endpoints M01=5 --endpoints M02=6 \
    --group-col condition --plot-dir figs --endpoint-line 1.0

id  time  predicted    lower    upper        model  actual
M01   5.0   0.932585 0.670443 1.194728 ARIMA(1,1,0)    1.00
M02   6.0   0.955696 0.748309 1.163084 ARIMA(1,1,0)    0.94

   group             id        model  n   rmse  picp  mpiw
endpoint            M01 ARIMA(1,1,0)  1 0.0674 100.0 0.524
endpoint            M02 ARIMA(1,1,0)  1 0.0157 100.0 0.415
endpoint -- endpoint --               2 0.0489 100.0 0.470
                OVERALL               2 0.0489 100.0 0.470

Report all three metrics together. RMSE is point accuracy, PICP the percentage of actual values inside the interval, and MPIW the mean interval width in RELSA units — a model can reach PICP = 100% by making the interval so wide it says nothing, which is exactly what the paper's pancreatic cancer row (PICP 100%, MPIW 7.35, i.e. 735% of the RELSA range) shows.

For live monitoring, forecast one step ahead at every time point instead:

python scripts/forecast_relsa.py relsa_scores.csv --mode rolling --animals M03

Two things to know before trusting a forecast:

day is far too sparse for ARIMA. It buys usable model selection and narrower intervals at the cost of honest uncertainty. Set --interpolate-step 0 when measurement frequency allows.

• Interpolation is on by default (--interpolate-step 0.1), because one measurement per

collapse in the last hours before an endpoint will not be forecast from a smooth prior trajectory — the paper's own failure case. Act on the upper bound of the interval, and never let a low forecast override an animal that looks unwell.