name: relsa-severity-assessment
description: Multivariate severity assessment and humane endpoint prediction for laboratory animal studies using the RELSA (RELative Severity Assessment) score and ARIMA-based foRcast forecasting. Use when combining welfare readouts — body weight or weight loss, body temperature, clinical or nesting scores, biomarkers, activity, heart rate, burrowing, wheel running — into one severity score per animal per day, when asking which animals are at risk of reaching a humane endpoint or when one will be reached, when defining attention/danger zones or thresholds on a severity scale by kernel density estimation, or when reporting severity for a 3Rs, refinement, animal-welfare, or EU Directive 2010/63/EU severity-assessment context. Covers directionality ("turned" variables), baseline normalization, reference sets, RELSA weights, ARIMA prediction intervals, and RMSE/PICP/MPIW evaluation.
license: MIT
allowed-tools: Read Write Edit Bash
compatibility: Requires Python >=3.10 with numpy, pandas, and scipy; statsmodels >=0.14 for forecasting and matplotlib for figures. Tested with numpy 2.5, pandas 3.0, scipy 1.18, statsmodels 0.14.6. No network access needed.
metadata:
    version: "1.0"
    skill-author: K-Dense Inc.

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:

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

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.

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

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.

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

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

per-animal severity score

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

the severity score at a coming time point

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

relative scale

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

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

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

relsa_score.py and kde_thresholds.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:

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

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

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

convention codes the baseline time point as -1.

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

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

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

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

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:

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

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.

• ARIMA cannot predict a cliff. It assumes stationarity and linearity, so an abrupt

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.

Step 3 — put the score in context with severity zones

python scripts/kde_thresholds.py relsa_scores.csv \
    --group treatment=treated --n-thresholds 2 --plot zones.png --json zones.json

KDE on 33 RELSA scores  (bandwidth = 0.1502)
  candidate thresholds (density minima): 0.703
  density modes: 0.264, 0.866
  normal    [0.000, 0.703)  n=25 (75.8%)
  danger    >= 0.703  n=8 (24.2%)

Thresholds are the *minima* of the score density — the sparse valleys between clusters of

scores. Include endpoint animals, survivors, and shams: the zones are meant to separate

those states, so all of them must be represented.

Check the bandwidth before believing a threshold. On the published sepsis data this

implementation finds minima at 0.355 and 0.655 (published: 0.337 and 0.643) — but a 10%

larger bandwidth removes both minima entirely. Run the sweep in

references/thresholds-and-zones.md and report the sweep, not a bare pair of numbers. An

empty threshold list is a legitimate answer: the scores form one cluster and there is no

data-driven place to cut.

Boundaries: state these when you report

• RELSA is an aid to severity assessment, not a decisive parameter. An animal with a low

RELSA score that shows other signs of distress must still be handled accordingly. Neither

procedure is a validated predictor of death.

• KDE zones are not regulatory severity gradings. EU Directive 2010/63/EU's categories

(non-recovery, mild, moderate, severe) are assigned prospectively by a different process.

The paper is explicit that its thresholds "should not be confused with regulatory severity

gradings" and are not directly translatable to them.

• Scores are not comparable across reference sets or models. RELSA is relative by

construction, and clinical scoring is not harmonized between laboratories. Always report

the reference set with the score.

• The published evidence is a proof of concept: 13 animals across seven models, five of

those rows resting on one or two animals. The overall RMSE of 0.069 and PICP of 96% come

from 13 endpoint predictions.

• An underestimated score is the dangerous error, because it discourages attention and can

delay a euthanasia decision, whereas an overestimate merely prompts extra care.

Reporting checklist

A severity analysis is reproducible only if all of this is stated:

1. Outcome measures, their units, and their directionality (which were turned, and why).
2. The baseline time point or window, and which variables were normalized.
3. Any score mapping applied to ordinal variables, with its scale.
4. The reference set: which animals, which group, how many, and why they are assumed to

carry the greatest burden.

1. Humane endpoint criteria actually applied in the study, separately from the RELSA score.
2. For forecasts: interpolation step, the selected ARIMA order per animal, and RMSE, PICP,

*and* MPIW.

1. For thresholds: the bandwidth, the number of scores, and a bandwidth sensitivity sweep.
2. Software versions, and the statement that thresholds are model-specific and not regulatory

gradings.

Common pitfalls

1. Wrong directionality — a rising variable not listed in --turned contributes exactly

zero, silently, and no warning is possible unless it never once falls. Check the reference

model table yourself: max reached should be below 100 for a falling variable and above 100

for a turned one, and max delta should be a plausible size for that measure.

1. Normalizing a percentage twice — bwc [%] and mapped scores are already on the percent

scale; passing them to --normalize flattens them.

1. A zero baseline — a clinical score of 0 makes the ratio undefined; the variable becomes

all-NaN with a warning. Use --score-scale.

1. A reference set that does not express the burden — a variable that never deviates in it

raises an error rather than dividing by zero, and one that barely deviates inflates every

score.

1. Changing variable composition along a trajectory — see decision 4 above.
2. Reading MPIW as a good thing — a wide interval raises PICP while destroying the

forecast's usefulness.

1. Reporting a KDE threshold without its bandwidth — thresholds can vanish under a 10%

bandwidth change.

1. Treating the forecast as permission to wait — the model cannot see abrupt

deterioration, and the humane endpoint criteria of the protocol always take precedence.

1. Comparing RELSA scores between models — only valid within one reference frame.

Resources

Scripts

• scripts/relsa_score.py — the RELSA procedure: prepare(), build_reference(),

relsa_scores(), relsa_weights(), and a ReferenceModel that serialises to JSON.

Reproduces the R package's published worked example to two decimals.

• scripts/forecast_relsa.py — the foRcast tool: auto_arima() (Hyndman–Khandakar stepwise

AICc selection), forecast_animal(), predict_endpoint(), rolling_forecast(),

forecast_indirect(), summarize(), and Figure-1-style plots.

• scripts/kde_thresholds.py — severity zones: bw_nrd0() (R's bandwidth), density_curve(),

find_thresholds(), zone assignment, and Figure-3-style density plots.

• scripts/_common.py — RELSA-format I/O, validation, score_to_percent(),

percent_of_baseline(), and forecast_metrics() (RMSE/PICP/MPIW).

References

• references/relsa-method.md — the four steps in full, the score/zero-baseline problem, the

variable-composition trap, parity notes against the R package, and the outcome measures and

endpoint criteria of all seven published models.

• references/forecasting.md — ARIMA selection, why interpolation is a distortion, direct vs

indirect prediction, the metrics, the published Table 1, and what this port reproduces.

• references/thresholds-and-zones.md — KDE method, published thresholds, the bandwidth

sensitivity sweep, the regulatory boundary, and alternatives when KDE gives nothing.

Assets

• assets/example_cohort.csv — synthetic 6-mouse cohort with temperature, body weight, a

clinical score, and a biomarker; illustrative only, not real data.

Related skills

• experimental-design, statistical-power — designing the study and sizing the groups.
• statsmodels, timesfm-forecasting — general time-series modelling.
• statistical-analysis, scientific-visualization — group comparisons and figures.

Key references

• Talbot, S. R. et al. (2022). RELSA — a multidimensional procedure for the comparative

assessment of well-being and the quantitative determination of severity in experimental

procedures. *Front. Vet. Sci.* 9:937711. R package: <https://github.com/mytalbot/RELSA>

• Lutscher, S. et al. (2026). Refining humane endpoint detection by time-series forecasting

and threshold definition using a multivariate severity score. *Front. Physiol.* 17:1869563.

• Hyndman, R. J. & Khandakar, Y. (2008). Automatic time series forecasting: the forecast

package for R. *J. Stat. Softw.* 27, 1–22.

• EU Commission (2010). Directive 2010/63/EU on the protection of animals used for scientific

purposes.