foRcast: ARIMA forecasting of RELSA trajectories

scripts/forecast_relsa.py ports the foRcast tool of Lutscher et al. (2026), Front. Physiol. 17:1869563 — an ARIMA model fitted per animal to its own RELSA trajectory, forecasting the score at the next time point (or at the humane endpoint) with a 95% prediction interval.

The purpose is triage, not automation: identify the individuals at risk of reaching a humane endpoint so handling personnel give them attention, while avoiding euthanising animals that would have recovered. It is a proof of concept on 13 animals across seven models, not a validated clinical tool.

Why ARIMA

ARIMA(p, d, q) combines an autoregressive part (p lags of the series), differencing (d, to remove trend and reach stationarity), and a moving-average part (q lags of the forecast errors). It needs nothing but the animal's own history, which suits single-animal severity assessment where each individual is its own control.

Model selection follows Hyndman & Khandakar (2008), i.e. forecast::auto.arima:

1. Choose d by successive KPSS tests (null = stationary; difference while it is rejected).
2. Fit four seed models — (2,d,2), (0,d,0), (1,d,0), (0,d,1) — with and without a constant/drift term.
3. Hill-climb from the best of those over neighbouring (p, q) and the drift term until AICc stops improving.

auto_arima(..., stepwise=False) searches the full p × q grid instead. Both are bounded by max_p, max_q, max_d; the paper notes that the globally best model could lie outside that range, which is a limitation of the approach rather than of one implementation.

Interpolation: the necessary distortion

Animal experiments typically produce one measurement per animal per day. ARIMA is conventionally said to want ~50 observations (Box et al., 2016), a number recently challenged (Hassouna & Al-Sahili, 2020) but still far above what a 7-day study yields. The paper's workaround is to interpolate linearly between observed values at 0.1-day increments and fit the model to that denser series, and it is explicit that this is an alteration of the method, not a free improvement:

• It raises autocorrelation and partial autocorrelation, which is what lets automatic order selection work at all on such short series.
• It narrows the prediction interval, improving coverage (PICP) at the cost of honestly representing uncertainty. The paper identifies interpolation as necessary "to minimize errors while maximizing prediction interval coverage with narrower boundaries".
• It adds no information. Interpolated points are a smoothness assumption, and a trajectory that actually moved non-linearly between measurements is misrepresented.

interpolate_step=None / --interpolate-step 0 fits the observed series directly. Prefer it whenever measurement frequency allows — with automated home-cage or telemetry monitoring the interpolation step becomes unnecessary, which is the paper's own outlook.

Forecast directly, not variable-by-variable

Two routes to a predicted RELSA score:

• Direct — forecast the RELSA series itself. forecast_animal(), predict_endpoint().
• Indirect — forecast each outcome measure, then compute RELSA from the forecasts. forecast_indirect().

The paper compared them in the sepsis model and direct won clearly: median deviation from the actual score −0.002 (direct) versus −0.240 (indirect), a large effect (d = 1.42, 95% CI [1.03, 1.81]). The reason is error propagation — each variable's forecast error accumulates through the score, whereas the direct forecast carries only its own error.

Use direct. forecast_indirect() exists to reproduce the comparison and to inspect which variable is driving a forecast.

Metrics

Reported together, because each hides a failure the others catch (_common.forecast_metrics):

MPIW is read against the RELSA scale, which normally spans about 0–1: the paper's overall MPIW of 1.69 means the average interval covered 169% of the RELSA range, and the pancreatic cancer model's 7.35 means 735% — a technically perfect PICP with almost no information in it. Always report MPIW next to PICP.

Published performance (Table 1)

Predicting the RELSA score at the (pre-)humane endpoint from all measurements up to the time point immediately before it:

Five of the seven rows rest on one or two animals. The overall PICP of 96% comes from 13 endpoint predictions.

What this port reproduces

Using the public sepsis data (tm_sepsis.txt, 7 mice) with the paper's four telemetry variables, no turned variables, and the CLP animals as reference set:

• Mouse ID_801 (the paper's Figure 1A): predicted RELSA 0.94 at the endpoint hour against an actual 0.93, RMSE 0.010, actual value inside the 95% interval. The published sepsis row is RMSE 0.009 over two animals.
• PICP 100% for both endpoint animals, matching the published row.
• MPIW 0.42–0.46 against a published 0.30 — this port's intervals are wider. The exact width depends on the interpolation step, the fitted variance, and the state-space implementation (statsmodels SARIMAX versus R's arima), so treat MPIW comparisons across implementations as approximate.

The paper's exact reference set and baseline window per model are in its Supplementary Table S2, which is not bundled here; small differences in those choices shift every score slightly.

Limits that matter more than the metrics

• ARIMA cannot predict a cliff. The model assumes stationarity and linearity. An abrupt collapse in the last hours before an endpoint is not forecastable from a smooth prior trajectory — this is the paper's own failure case (Figure 1C, the DSS blood-sampling mouse whose pre-endpoint score rose sharply and fell outside the 95% bounds). For sudden change, the paper points to Bayesian online changepoint detection (Adams & MacKay, 2007) or Markov switching models (Hamilton, 2020) as alternatives.
• An underestimated score is the dangerous error. An overestimate merely prompts extra attention; an underestimate discourages personnel from giving an animal the attention it needs and can delay a euthanasia decision. Asymmetric consequences deserve asymmetric handling: act on the upper bound of the interval.
• RELSA is a severity-assessment aid, not a decision rule. An animal with a low RELSA score that shows other signs of distress must still be handled accordingly. The paper is explicit that RELSA is "intended as an aid to severity assessment rather than a decisive parameter", and the RELSA package's own documentation states it is not a predictor of death.
• Two prior measurements are not enough. The paper's largest direct-prediction errors (Δ = 0.76 and 0.74) came from forecasts made at the earliest possible time point with only two prior observations. forecast_animal() records a warning below four observed points.
• Parameter volatility hurts. Activity forecast worst of the sepsis variables, being both intrinsically volatile and measured at low frequency. Including a noisy variable in the multivariate RELSA score mitigates its noise — one argument for the composite over single-parameter forecasting.

Key references

• Hyndman, R. J. & Khandakar, Y. (2008). Automatic time series forecasting: the forecast package for R. J. Stat. Softw. 27, 1–22.
• Hyndman, R. J. & Athanasopoulos, G. (2021). Forecasting: Principles and Practice, 3rd ed.
• Khosravi, A. et al. (2011). Comprehensive review of neural network-based prediction intervals. IEEE Trans. Neural Netw. 22, 1341. (PICP/MPIW)
• Pang, J. et al. (2018). Optimize the coverage probability of prediction interval for anomaly detection of sensor-based monitoring series. Sensors 18, 967.
• Petrică, A. et al. (2016). Limitation of ARIMA models in financial and monetary economics. Theor. Appl. Econ. 23, 19–42.