relsa-severity-assessment

RELSA 严重程度评估与人道终点预测

概述

在动物研究中，严重程度评估不仅是法律要求，也是科学研究的重要支撑： 它会影响人道终点的决策，而福利监测不佳则会降低研究的可重复性。 通常的做法是单独评估每项读数——这里看体重下降，那里看临床评分—— 因此很难判断某只动物的实际状况有多严重。

此技能实现了两种已发表的方法来解决这一问题：

• RELSA（Talbot 等，2022）将多项结果指标合并为每只动物在每个时间点的一个评分， 并以已知负担的参考集为基准进行表示。RELSA = 0 表示基线；RELSA = 1 表示动物已达到 参考集中的最大偏离程度。
• foRcast（Lutscher 等，2026）为单只动物的 RELSA 轨迹拟合 ARIMA 模型，并预测下一个评分及其 95% 预测区间， 从而能在动物到达人道终点前识别出可能即将到达终点的动物。通过在 RELSA 尺度上进行核密度估计， 可得出用于解读的候选关注区间和危险区间。

重点在于优化：更早关注有风险的动物，并避免对原本能够恢复的动物实施安乐死。这两种方法都是严重程度评估的辅助工具， 而非决策规则——请参阅适用边界。

何时使用此技能

• 将体重下降、体温、临床评分、生物标志物或遥测数据合并为单只动物的严重程度评分
• 判断队列中的哪些动物有可能到达人道终点，或预测未来某个时间点的严重程度评分
• 在共同的相对尺度上比较不同治疗组、干预措施或动物模型的严重程度
• 根据数据定义严重程度尺度上的阈值或区间
• 撰写动物福利报告、3Rs/优化分析或依据 EU Directive 2010/63/EU 提交的申请中的严重程度评估部分

对于不属于严重程度评分的时间序列的一般预测，请使用 timesfm-forecasting 或 statsmodels。对于研究设计和样本量，请使用 experimental-design 和 statistical-power。

安装

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

relsa_score.py 和 kde_thresholds.py 仅需 numpy/pandas/scipy；预测需要 statsmodels，绘图则仅需 matplotlib。

数据格式

CSV 文件中每只动物每个时间点占一行：

id	treatment	condition	day	temp	weight	score	il6
M01	treated	endpoint	-1	37.15	25.17	0	35.1
M01	treated	endpoint	0	37.26	25.25	0	39.5
M01	treated	endpoint	1	35.83	23.12	4	162.0

• 必须包含 id 和时间列（day、time、hour 等）；treatment 和 condition 是可选标签，可用于分组和选择参考集。
• 时间可以用天、小时或分钟表示——只需确保每只动物的时间都单调递增。RELSA 的惯例是将基线时间点编码为 -1。
• **每只动物每个时间点占一行。**请先将每小时遥测数据平均为每个区间的一个值（已发表的模型会对心率、HRV 和体温取平均，并对活动量求和）。
• 缺失测量值请留空。它们会从评分中剔除，而不会进行插补——将缺失值视为“无偏离”会导致严重程度被低估。

assets/example_cohort.csv 是一个小型合成队列（6 只小鼠，9 天，包含体温、体重、0–8 临床评分和类似 IL-6 的生物标志物）， 供以下每条命令使用，因此所有命令均可直接运行。

决定结果的四个要点

明确作出这些决定，并将其写入方法中。程序的其他方面都没有这些重要。

1. 方向性 — 哪些变量会在状况恶化时升高？ 默认是下降（体重、活动量、食物摄入量、掘洞行为、跑轮活动）。必须声明会升高的变量为 --turned：临床评分、炎症生物标志物、发热、心动过速。如果弄错，变量将完全不产生作用，因为“错误”方向上的偏差会被归零。体温取决于模型——在败血症和内毒素血症模型中会下降，在发热模型中会升高。数据无法替你作出这个判断：在已发表的败血症模型中，活动量确实可能高于基线，且超出的幅度大于低于基线的幅度。因此，只有从未沿声明方向变化过的变量才可被检测出来，而 build_reference() 会专门针对这种情况发出警告。

2. 参考组 — 相对于什么？ 没有参考组，RELSA 分数毫无意义。使用模型中预计负担最重的组（已发表的研究使用最高剂量组或达到终点的治疗组）。参考组太轻，所有分数都会高于 1；太重，则所有分数都会被压缩到接近 0。使用 --save-reference 保存参考组，并用 --load-reference 复用，这样后续队列就能保持在相同的量表上。

3. 基线为零的分数。 健康动物的临床评分为 0，无法进行比值归一化——0/0 无定义。改用 --score-scale score=8 来映射该评分的量表（健康 → 100%，最差状态 → 200%），这也会将其标记为升高变量。这种映射体现了一项建模选择：一个评分点相对于体重百分之一的价值。请说明这项选择。另一种做法是将该评分排除在 RELSA 之外，并将其作为独立的终点标准。

4. 全程测量哪些变量。 由于分数会对当时可用的变量取平均，某个变量在轨迹中途出现或消失，本身就会改变分数。在已发表的败血症数据中，加入仅在安乐死当天记录的体重，会使该动物的终点分数从 0.93 降至 0.83，但没有任何生物学原因。relsa_scores() 会在变量构成发生变化时发出警告；请使用全程都有记录的变量进行评分。

工作流

步骤 1 — 计算 RELSA 分数

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

系统会输出参考模型，以便审核量表：

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 包含各变量的权重及其对应分数，因此可以解释分数——例如，M01 持续恶化直至终点，而 M03 在第 3 天达到峰值后恢复：

 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

权重为 1.00 表示该变量达到参考组最大值；n_vars 表示该时间点有多少个变量纳入了分数计算。

如果需要这些对象，也可以通过 Python 执行相同操作：

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)

步骤 2 — 预测终点

使用终点前一个时间点之前的所有数据进行训练，预测终点的分数，并评估预测结果：

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

同时报告这三项指标。RMSE 表示点预测准确度，PICP 表示实际值落在区间内的百分比，MPIW 表示以 RELSA 单位计的平均区间宽度——模型可以通过把区间设得极宽来达到 PICP = 100%，但这就毫无意义了；论文中的胰腺癌一行正是如此（PICP 100%，MPIW 7.35，即 RELSA 范围的 735%）。

对于实时监测，则改为在每个时间点预测下一步：

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

在信任预测之前，有两点需要了解：

• 默认启用插值（--interpolate-step 0.1），因为每天只测量一次，对 ARIMA 来说过于稀疏。这样能进行有效的模型选择并缩窄区间，但代价是低估真实的不确定性。如果测量频率允许，请设置 --interpolate-step 0。
• **ARIMA 无法预测骤变。**它假设数据具有平稳性和线性，因此在终点前的最后几小时内突然崩溃时，模型无法根据平稳的先前轨迹预测到这一点——论文自身的失败案例就是如此。请依据区间的上限采取行动，绝不要因为预测值低，就忽视看起来不适的动物。

步骤 3 — 根据严重程度分区解读分数

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%)

阈值是分数密度的极小值——也就是分数聚类之间的稀疏谷地。请纳入终点动物、存活动物和假手术动物：这些分区旨在区分上述状态，因此必须包含所有这些动物。

**相信阈值之前，先检查带宽。**在已发表的脓毒症数据上，本实现找到的极小值位于 0.355 和 0.655（已发表结果：0.337 和 0.643），但带宽增加 10% 后，这两个极小值都会完全消失。请运行 references/thresholds-and-zones.md 中的扫描，并报告扫描结果，而不是只给出一对数字。阈值列表为空也是合理的结果：分数只形成一个聚类，没有数据驱动的切分位置。

局限：报告结果时请说明

• **RELSA 是严重程度评估的辅助工具，不是决定性指标。**如果 RELSA 分数较低的动物表现出其他不适迹象，仍必须采取相应措施。两种方法都不是经过验证的死亡预测指标。
• **KDE 分区不是监管严重程度分级。**欧盟指令 2010/63/EU 中的类别（非恢复、轻度、中度、重度）由不同流程预先确定。论文明确指出，其阈值“不得与监管严重程度分级混淆”，也不能直接转换为这些分级。
• **不同参考集或模型的分数不可比较。**RELSA 从定义上就是相对指标，且不同实验室的临床评分并未统一。报告分数时，务必同时报告参考集。
• 已发表的证据只是概念验证：涵盖 7 种模型中的 13 只动物，其中 5 行数据仅基于 1 到 2 只动物。总体 RMSE 0.069 和 PICP 96% 来自 13 次终点预测。
• 低估分数是更危险的错误，因为这会降低对动物的关注，并可能延误安乐死决策；而高估只会促使采取额外的照护措施。

报告检查清单

只有明确说明以下所有内容，严重程度分析才可复现：

1. 结局指标、单位及其方向性（哪些指标进行了反向处理，以及原因）。
2. 基线时间点或时间窗口，以及进行了归一化的变量。
3. 对有序变量应用的任何分值映射及其量表。
4. 参考集：包含哪些动物、属于哪个组、数量是多少，以及为何认为它们承受的负担最大。
5. 研究中实际采用的人道终点标准，并与 RELSA 分值分开说明。
6. 对于预测：插值步长、每只动物所选的 ARIMA 阶数，以及 RMSE、PICP 和 MPIW。
7. 对于阈值：带宽、分值数量，以及带宽敏感性扫描。
8. 软件版本，并声明阈值特定于模型，且不属于监管评级。

常见问题

1. 方向性错误 — 未列入 --turned 的上升变量会被静默地计为零，除非它从未下降过，否则不会发出警告。请自行检查参考模型表：对于下降变量，max reached 应低于 100；对于反向处理的变量，应高于 100；max delta 应与该指标的合理量级相符。
2. 百分比重复归一化 — bwc [%] 和映射后的分值已经采用百分比量表；将它们传给 --normalize 会压平数值。
3. 基线为零 — 临床评分为 0 时，比例无定义；该变量会变成全为 NaN，并发出警告。请使用 --score-scale。
4. 无法体现负担的参考集 — 如果某个变量在参考集中始终没有偏离，会报错而不是除以零；如果偏离很小，则会放大所有分值。
5. 沿轨迹改变变量组成 — 参见上文决策 4。
6. 把 MPIW 当作好事 — 区间越宽，PICP 越高，但预测的实用性也会随之降低。
7. 报告 KDE 阈值却不报告带宽 — 带宽变化 10% 就可能导致阈值消失。
8. 将预测视为可以等待的许可 — 模型无法识别突发恶化，且方案中的人道终点标准始终优先。
9. 比较不同模型的 RELSA 分值 — 只有在同一个参考框架内比较才有效。

资源

脚本

• scripts/relsa_score.py — RELSA 流程：prepare()、build_reference()、 relsa_scores()、relsa_weights()，以及可序列化为 JSON 的 ReferenceModel。 可复现 R 包已发布的示例，精确到小数点后两位。
• scripts/forecast_relsa.py — foRcast 工具：auto_arima()（Hyndman–Khandakar 逐步 AICc 选择）、forecast_animal()、predict_endpoint()、rolling_forecast()、 forecast_indirect()、summarize() 和 Figure 1 风格的图表。
• scripts/kde_thresholds.py — 严重程度分区：bw_nrd0()（R 的带宽）、density_curve()、 find_thresholds()、分区分配，以及 Figure 3 风格的密度图。
• scripts/_common.py — RELSA 格式的 I/O、验证、score_to_percent()、 percent_of_baseline() 和 forecast_metrics()（RMSE/PICP/MPIW）。

参考资料

• references/relsa-method.md — 完整介绍四个步骤、评分/基线为零的问题、变量组成陷阱、 与 R 包的对照说明，以及所有七个已发布模型的结局指标和终点标准。
• references/forecasting.md — ARIMA 选择、插值为何会造成失真、直接与间接预测、 各项指标、已发布的表 1，以及此移植版本复现的内容。
• references/thresholds-and-zones.md — KDE 方法、已发布阈值、带宽敏感性扫描、 监管界限，以及 KDE 无结果时的替代方法。

资源文件

• assets/example_cohort.csv — 合成的 6 只小鼠队列，包含体温、体重、临床评分和生物标志物；仅供示例，不是真实数据。

相关技能

• experimental-design、statistical-power — 设计研究并确定组别样本量。
• statsmodels、timesfm-forecasting — 通用时间序列建模。
• statistical-analysis、scientific-visualization — 组间比较和绘图。

关键参考文献

• Talbot, S. R. et al. (2022). RELSA — 一种多维程序，用于比较评估动物福利，并定量确定实验程序的严重程度。Front. Vet. Sci. 9:937711。R 包：https://github.com/mytalbot/RELSA
• Lutscher, S. et al. (2026). 通过时间序列预测和使用多变量严重程度评分设定阈值，改进人道终点检测。Front. Physiol. 17:1869563。
• Hyndman, R. J. & Khandakar, Y. (2008). 自动时间序列预测：用于 R 的 forecast 包。J. Stat. Softw. 27, 1–22。
• EU Commission (2010). 关于保护用于科学目的动物的 Directive 2010/63/EU。

Agent 操作流程

1. 检查环境。 明确研究目标、约束条件（预算、时间、伦理）以及用户需要作出的决策。
2. 确定输入。 根据数据或用户提供的信息，确认格式、标识符和参数。任何会影响结果的数值，都应询问，而不是猜测。
3. 先运行小规模版本。 在扩展任何选项之前，提出一份简短的候选方案清单，并明确说明假设。
4. 使用上述说明和参考资料执行完整任务。
5. 验证结果。 明确陈述假设；依据明确标准比较备选方案；统计分析计划应包含统计功效和分析选项。
6. 报告。 说明运行了什么（版本、命令、参数）、检查了什么，以及仍有哪些不确定之处。

如果发生这种情况	请执行以下操作
设计选择取决于未知的效应量或成本	给出多种情景（乐观、预期、悲观），不要只作单一猜测。
当前安装的版本中缺少这些说明所述的函数、标志或终点	查阅已安装版本自己的文档（help()、--help、官方文档），进行调整并告知用户。绝不可捏造 API。
必需的输入、标识符或参数不明确	询问用户，或在运行前明确说明假设。

诚信原则

• 绝不捏造结果、参数、标识符、引用或统计数据。如果无法运行或验证某项内容，请直言说明。
• 对推测内容作出标注；将既有证据与新假设分开。
• 将此处与特定版本相关的细节视为可能过时：请根据已安装版本的官方文档确认。
• 对于会产生费用、占用共享 GPU 或云配额、涉及个人或患者数据，或无法撤销的操作，须先征得同意。