代码
library(ivdtools)
library(readr)4 稳定性研究
4.1 分析目的
稳定性研究用于评价试剂、校准品或质控品在规定储存、运输或使用条件下，检测结果随时间 变化的程度。本文演示如何使用 ivdtools 包建立稳定性分析计划，完成节点偏差、条件间及 相对首节点的回归趋势、 达到允许限值的时间预测、阿伦纽斯（Arrhenius）加速稳定性分析及平均动力学温度（mean kinetic temperature，MKT）计算。
本文数据均为确定性的教学示例，仅用于验证分析流程。示例中的 5% 允许漂移、10% 加速 降解限值、目标温度和活化能不是通用接受标准；正式研究必须使用方案、风险评估或适用标准 预先规定的参数。选择浓度值或信号值作为响应值应根据方案确定。
4.2 函数概述
| 函数 | 主要用途 | 关键输入或输出 |
|---|---|---|
stability_plan() |
规划回归系列的最少时间点 | 允许漂移、变异、重复数、预期漂移和效能 |
stability_bias() |
计算各条件相对首节点的绝对/相对偏差 | 节点汇总、条件内偏差、条件间偏差和限值状态 |
stability_regression() |
单条件趋势或试验/参考条件配对偏差回归 | 斜率、单侧置信限、节点结果和诊断信息 |
stability_time() |
预测点估计或单侧置信限达到允许偏差的时间 | 预测时间、是否外推及搜索范围 |
arrhenius() |
根据多个升高温度的降解速率外推目标温度稳定期 | 动力学级数、速率、活化能和预测稳定期 |
mkt() |
计算等间隔或时间加权的平均动力学温度 | 各探头及合并 MKT |
stability_regression() 支持两种模式：
mode = "self"：对单一条件的测量结果随时间回归；mode = "compare"：在匹配时间点计算试验条件相对参考条件的偏差，再对偏差随时间回归。
分析前应明确条件映射、时间单位、偏差尺度、预期变化方向和允许限。只有在方案提供接受 标准时，才能据此作出是否可接受的结论。
4.3 示例一：建立分析计划
假设研究方案规定允许相对漂移为 5%，预期漂移为 2.5%，重复性变异为 1.5%，目标效能为 90%。比较每个时间点测量 1、2 或 3 次时所需的最少时间点：
代码
plan <- ivdtools::stability_plan(
allowable_drift = 5,
variability = 1.5,
replicates = c(1, 2, 3),
expected_drift = 2.5,
power = 0.9,
mode = "simple",
bias_type = "relative"
)
print(plan)
Stability Study Time-Point Plan
Allowable drift: 5.0000; expected drift: 2.5000
Success target: 90%; method: CLSI EP25 Appendix A Table A2
Replicates Residual Allow/Resid Min nodes
------------------------------------------------
1 1.5000 3.3333 n/a
2 1.0607 4.5455 28
3 0.8660 5.7471 18
n/a -- No feasible entry in EP25 table.4.4 示例二：偏差分析
4.4.1 读取和核验原始数据
代码
bias_data <- read_csv(
"./data/stability-bias-example.csv",
show_col_types = FALSE
)
bias_data代码
dim(bias_data)[1] 30 4代码
str(bias_data)spc_tbl_ [30 × 4] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
$ Condition: chr [1:30] "2-8C" "2-8C" "2-8C" "2-8C" ...
$ Day : num [1:30] 0 0 0 7 7 7 14 14 14 21 ...
$ Replicate: num [1:30] 1 2 3 1 2 3 1 2 3 1 ...
$ Result : num [1:30] 100.2 99.8 100 99.9 100.1 ...
- attr(*, "spec")=
.. cols(
.. Condition = col_character(),
.. Day = col_double(),
.. Replicate = col_double(),
.. Result = col_double()
.. )
- attr(*, "problems")=<externalptr> 4.4.2 节点偏差计算
以 2-8C 为参考条件，按各条件自身首节点归一化，评价相对偏差是否在示例性的 ±5% 限值内：
代码
bias_result <- ivdtools::stability_bias(
data = bias_data,
condition = "Condition",
time = "Day",
value = "Result",
reference_condition = "2-8C",
limit = 5,
bias_type = "relative",
time_unit = "d"
)
print(bias_result)
Stability Bias Analysis
Conditions: 2 (2-8C, 25C)
Reference: 2-8C
Time unit: d
Node comparisons:20
Acceptance: +/- 5.0000 (relative); 10/10 nodes passed
Within-condition bias (vs first node)
Cond Time Value Bias(%) Status
----------------------------------------------
2-8C 0.0000 100.0000 0.0000 Pass
2-8C 7.0000 99.9333 -0.0667 Pass
2-8C 14.0000 99.7333 -0.2667 Pass
2-8C 21.0000 99.5000 -0.5000 Pass
2-8C 28.0000 99.3333 -0.6667 Pass
25C 0.0000 100.0000 0.0000 Pass
25C 7.0000 99.2333 -0.7667 Pass
25C 14.0000 98.3333 -1.6667 Pass
25C 21.0000 97.2333 -2.7667 Pass
25C 28.0000 95.8333 -4.1667 Pass
Between-condition bias (reference: 2-8C)
Time Condition TestVal Bias(%)
--------------------------------------
0.0000 25C 100.0000 0.0000
7.0000 25C 99.2333 -0.7005
14.0000 25C 98.3333 -1.4037
21.0000 25C 97.2333 -2.2781
28.0000 25C 95.8333 -3.5235条件内偏差以各条件自己的首节点为基线；条件间偏差则在相同时间点比较试验条件和参考条件。
代码
plot(bias_result)
4.5 示例三：不同保存条件比较的回归分析
4.5.1 读取和核验原始数据
代码
regression_data <- read_csv(
"./data/stability-regression-example.csv",
show_col_types = FALSE
)
regression_data4.5.2 试验条件与参考条件比较
在每个匹配时间点计算 25C 相对 2-8C 的相对偏差，并对偏差随时间进行回归。示例 预先指定结果随时间降低，因此使用 direction = "decrease"：
代码
regression_result <- ivdtools::stability_regression(
data = regression_data,
time = "Day",
value = "Result",
condition = "Condition",
mode = "compare",
test_condition = "25C",
reference_condition = "2-8C",
bias_type = "relative",
direction = "decrease",
conf.level = 0.95,
limit = 5,
time_unit = "d"
)
print(regression_result)
Stability Regression
Mode: compare
Test condition: 25C
Reference: 2-8C
Fit: WLS (1/SE^2); intercept = 0.0799, slope = -0.1025
Direction: decrease
One-sided CI: 95.0%
R-squared: 0.9956
Cook distance: row(s) 6 > 1; not removed
Time ObsBias Estimate CI_limit Status
--------------------------------------------------
0.0000 0.0000 0.0799 -0.0750
7.0000 -0.6669 -0.6376 -0.7566
14.0000 -1.3013 -1.3552 -1.4522
21.0000 -2.0040 -2.0728 -2.1717
28.0000 -2.6747 -2.7904 -2.9140
42.0000 -4.3842 -4.2256 -4.4296 代码
plot(regression_result)
模型使用单侧置信限评价最不利方向。本示例识别出一个 Cook 距离较大的节点，但函数不会 自动删除；应回查原始记录、实验过程和方案规定，并在有充分理由时进行包含/不包含该节点 的敏感性分析。
4.6 示例四：与首节点比较的回归分析
4.6.1 单个条件相对首节点分析
除示例三的试验条件/参考条件配对回归外，mode = "self" 可分析单个条件自身随时间的变化。 以下仍使用回归示例数据，但仅选择 25C 条件；各节点的 observed_bias 相对于该条件首节点 均值计算：
代码
self_regression_result <- ivdtools::stability_regression(
data = regression_data,
time = "Day",
value = "Result",
condition = "Condition",
mode = "self",
test_condition = "25C",
bias_type = "relative",
direction = "decrease",
conf.level = 0.95,
limit = 5,
time_unit = "d"
)
print(self_regression_result)
Stability Regression
Mode: self
Test condition: 25C
Fit: WLS (1/SE^2); intercept = 100.0849, slope = -0.1115
Direction: decrease
One-sided CI: 95.0%
R-squared: 0.9966
Cook distance: row(s) 1, 6 > 1; not removed
Time ObsBias Estimate CI_limit Status
--------------------------------------------------
0.0000 0.0000 -0.0000 -0.0000
7.0000 -0.7000 -0.7797 -0.8278
14.0000 -1.4000 -1.5594 -1.6555
21.0000 -2.2000 -2.3391 -2.4833
28.0000 -2.9667 -3.1188 -3.3110
42.0000 -4.7667 -4.6782 -4.9665 代码
plot(self_regression_result)
首节点观测偏差固定为 0；模型估计值则由全部时间点共同拟合，不要求回归截距恰好通过首节点。 单侧置信限用于评价预先指定的最不利方向，示例性的 ±5% 限值仍不能替代研究方案规定。
4.7 示例五：预测达到限值的时间
stability_time() 使用示例三的回归对象，分别计算点估计和单侧置信限首次达到 -5% 偏差的时间：
代码
time_result <- ivdtools::stability_time(
regression_result,
limit = 5,
max_time = 90
)
print(time_result)
Stability Limit Time
Limit: 5.0000 (relative); direction: decrease
estimate 49.5548 d [EXTRAPOLATED]
one_sided_ci 47.2316 d [EXTRAPOLATED]代码
plot(time_result)
本例最长实际观察时间为 42 天，而预测达到限值的时间超过了观察范围，因此输出标记为 EXTRAPOLATED。外推结果依赖线性趋势在未观察区间继续成立的假设，不能等同于实测稳定性 证据；应显著报告外推距离和不确定性，并优先用更长的实际观察数据确认。
4.8 示例六：阿伦纽斯（Arrhenius）加速稳定性分析
4.8.1 读取和核验原始数据
代码
arrhenius_data <- read_csv(
"./data/stability-arrhenius-example.csv",
show_col_types = FALSE
)
arrhenius_data4.8.2 动力学和阿伦纽斯拟合
示例使用 30、37 和 45 °C 的加速数据，外推至 5 °C；允许相对变化 10%，并预先指定 降解方向为降低：
代码
arrhenius_result <- ivdtools::arrhenius(
data = arrhenius_data,
temperature = "Temperature",
time = "Day",
value = "Result",
target_temp = 5,
order = "auto",
direction = "decrease",
limit = 10,
bias_type = "relative",
temp_unit = "C",
time_unit = "d",
conf.level = 0.95
)
print(arrhenius_result)
Arrhenius Stability Analysis
Kinetic order: zero (auto-selected)
Direction: decrease
Activation Ea: 97.159 kJ/mol
Target rate: 0.002346 /d
Limit time: 4263.7878 d [1585.3589, 11467.3633]代码
plot(arrhenius_result)
order = "auto" 会在零级和一级动力学候选中按实现规则选择模型。正式分析应在实验设计和 机理知识支持下确定候选级数，并核对每个温度的降解方向、速率显著性、降解程度和模型残差。 从升高温度外推至目标储存温度可能跨越很大范围，预测稳定期及其区间必须与模型假设、实际 时间数据和验证研究一起解释。
4.9 示例七：MKT 计算
4.9.1 读取和核验原始温度数据
代码
mkt_data <- read_csv(
"./data/stability-mkt-example.csv",
show_col_types = FALSE
)
mkt_data$Timestamp <- as.POSIXct(
mkt_data$Timestamp,
format = "%Y-%m-%d %H:%M:%S",
tz = "UTC"
)
mkt_data4.9.2 时间加权 MKT
使用两个温度探头和实际时间间隔，温度单位为摄氏度，活化能采用示例值 83.144 kJ/mol：
代码
mkt_result <- ivdtools::mkt(
data = mkt_data,
temp_cols = c("Probe_A", "Probe_B"),
time = "Timestamp",
temp_unit = "C",
ea = 83.144
)
print(mkt_result)
Mean Kinetic Temperature
Ea: 83.144 kJ/mol
Method: Time-weighted trapezoidal MKT
Probe_A 299.538 K 26.388 C
Probe_B 299.194 K 26.044 C
Combined 299.367 K 26.217 CMKT 不是温度的简单算术平均值，高温暴露会因阿伦纽斯关系获得更大权重。正式计算应确认 探头校准状态、时间顺序、缺失区间、温度单位、曝光持续时间以及活化能假设。MKT 只能概括 温度暴露，不能替代产品特异的稳定性研究。
4.10 多样本数据拆分
当不同浓度水平、样本类型或批次需要分别评价时，应保留样本标识，先拆分数据，再对每个 样本使用相同且预先规定的参数分析。读取并核验多样本原始数据：
代码
multiple_samples_data <- read_csv(
"./data/stability-multiple-samples-example.csv",
show_col_types = FALSE
)
multiple_samples_data使用 split() 建立按样本命名的数据列表，并先查看其中一个样本：
代码
sample_data_list <- split(
multiple_samples_data,
multiple_samples_data$Sample
)
names(sample_data_list)[1] "Sample_A" "Sample_B"代码
sample_a_data <- sample_data_list[["Sample_A"]]
sample_a_data逐个样本执行分析。
4.11 结果判读要点
- 方案优先：分析前明确条件、节点、重复数、时间单位、方向、偏差尺度和接受限。
- 原始数据只读：报告缺失、重复、非数值记录和每一步排除行，不自动删除异常记录。
- 配对正确：条件比较必须在匹配时间点进行，并明确试验条件和参考条件。
- 模型诊断：报告斜率、置信限、拟合方法、残差、影响点及是否发生 WLS 回退。
- 外推醒目：超过观察范围的稳定时间预测应明确标注，不作为实测证据。
- 温度假设明确：阿伦纽斯分析和 MKT 必须记录摄氏/开尔文单位、时间单位和活化能。
- 统计与接受性分开：统计显著不等于产品可接受，只能依据预设允许限作出判断。
- 独立复核：临床、监管或放行决策应由合适的专业人员独立复核。