code{white-space: pre;}
Pius Dahinden 博士，Tillotts Pharma 分析科学经理
抽象的
自 2014 年新的欧盟 GMP 指南（EudraLex 第 4 卷）第 6 章（质量控制）生效以来，对药品稳定性数据进行趋势分析是一项额外的分析质量控制要求。欧洲合规学院 (ECA) 分析质量控制工作组 (AQCWG) 制定了评估超出预期 (OOE) 和超出趋势 (OOT) 结果的标准操作程序。对于稳定性数据的评估，建议使用回归控制图 (RCC) 程序将当前批次与历史批次的数据进行比较。但是，从历史批次数据本身估计的斜率和/或截距存在不确定性。ECA AQCWG 提出了一种简化的随机系数回归 (RCR) 方法，该方法考虑到了这一点。介绍了一种新开发的 JMP 脚本，该脚本实现了推荐的 RCC / RCR 程序。该脚本使稳定性数据中的 OOE / OOT 点的识别变得简单，并在输出中提供所有必要的信息以及适当的图形显示。
介绍
制造商采用统计过程控制 (SPC) 来确保所生产产品的质量，并防止因生产过程失控而造成的经济损失。过程控制意味着定期监控正在运行的过程，如果任何受监控的参数出现异常趋势，则采取措施使过程重新受控。SPC 是一种工具，可让制造商及早发现过程中的变化，从而让他们能够调查导致问题的问题并在早期采取纠正措施。控制图方法是一种强大的过程监控工具，由 Shewhart [6] 推广。图 1 显示了典型控制图的示意图。
与不应存在趋势的生产过程相反，许多数据管理问题都具有固有趋势。对于此类数据管理问题的分析，传统的控制图方法与回归分析相结合，这种方法被 Mandel [4] 称为“回归控制图”(RCC)。作为控制限，使用公差区间 (TI) 是合理的。为了监测在推荐的储存条件下储存的药品的稳定性数据（例如 25 °C / 60% RH），RCC 用于监测活性药物成分的降解或杂质随时间的积累，以支持产品的保质期，即确认药品在标示的保质期内的稳定性。
Figure_1图 1图 1
图1 典型控制图组成
控制图通常显示值XX或样本质量特征的多个值（例如平均值）相对于样本数量或时间的统计数据。一般来说，图表包含一条代表平均值的中心线X¯\overset{\overline{}}{X}受控过程的水平线。一个基本假设是数据呈正态分布。另外两条特征水平线是控制下限和上限（分别为 LCL 和 UCL），距离中心线 ± 3 个标准差 (SD)，因此，只要过程处于受控状态，99.7% 的数据点就会落在这些限制范围内，这是由于随机变化造成的，这被称为常见原因变化。距离中心线 ± 2 SD 的线称为警告下限和上限（分别为 LWL 和 UWL），而距离中心线 ± 1 SD 的线没有特定名称。借助这些线，可以将控制图划分为 C 到 A 区。检测特殊原因变化的规则（表示过程不受控制）使用这些区域中的观察频率作为决策基础。
程序和结果
RCC 应用于药品稳定性数据旨在回答两个问题：1）就历史批次而言，当前批次是否超出预期（OOE）或超出趋势（OOT）（即关于从历史数据得出的 TI），2）就当前批次而言，最近测量的数据点是否为 OOE（即关于其预测区间（PI））？
就问题 1 而言，要评估趋势，至少需要当前批次的三个数据点。对于问题 2 的评估，至少需要四个数据点，因为最近测量的数据点不应用于构建 PI。
由于各个历史批次的截距和斜率很可能存在差异，因此普通最小二乘模型（共同斜率 - 共同截距）可能没有意义。因此，欧洲合规学院 (ECA) 分析质量控制工作组 (AQCWG) 建议使用“随机系数回归” (RCR) 方法，该方法可以解释各个批次的截距和斜率的变化 [7]。简化的 RCR 模型（改良的 Carter 和 Yang [1] 模型）包含两个未知参数，即截距β1\beta_{1}和斜率β2\beta_{2} .该程序符合 WHO [9] 和 ICH [2] 指导原则。
英国药品和保健产品管理局 (MHRA) 将 OOT 定义为与其他稳定性批次相比或与稳定性研究期间收集的先前结果相比未遵循预期趋势的稳定性测试结果 [8]。MHRA 将 OOS 定义为不符合预定验收标准的测试结果或超出官方药典和/或公司文件中确定的验收标准的测试结果 [8]。目前尚无 OOE 一词的正式定义。ECA AQCWG 建议（针对定性和定量测试）将 OOE 结果定义为不能归类为超出规格 (OOS) 或 OOT 的异常、意外或不寻常的发现 [7]。
RCR 模型
简化的 RCR 模型可以写成
和伊=𝛃1i+𝛃2i•𝐭伊+𝛆伊\mathbf{y}_{\text{ij}} = \mathbf{\beta}_{1i} + \mathbf{\beta}_{2i} \bullet \mathbf{t}_{\text{ij}} + \mathbf{\varepsilon}_{\text{ij}}
在哪里和i\mathbf{y}_{i}是ni×1n_{i} \times 1向量表示i日第i个批次nin_{i}重复测量（其中j=1,⋯,nij = 1,\ \cdots,\ n_{i} ， 和nin_{i}表示批次重复测量的总次数ii ）对响应变量（例如含量[％]，杂质总和[％]或溶出度[％]）， 𝛃1i\mathbf{\beta}_{1i}和𝛃2i\mathbf{\beta}_{2i}分别是批次特定截距和斜率的向量， 𝐭i\mathbf{t}_{i}是i日第i个批次的时间向量（例如 [Month]），以及𝛆i\mathbf{\varepsilon}_{i}是ni×1n_{i} \times 1向量nin_{i}随机误差。该 RCR 模型的先决条件是响应和解释变量之间存在线性关系。如果关系不是线性的，则必须在执行 RCR 建模之前应用适当的转换。
假设每个ε伊\varepsilon_{\text{ij}}来自均值为零、方差为σ2\sigma^{2}以及所有ε伊\varepsilon_{\text{ij}}彼此独立（条件独立），并且与模型系数无关β1i\beta_{1i}和β2i\beta_{2i} 。这一假设是合理的，因为分析方法的变异性不太可能取决于被测试的批次或执行测试的时间点。进一步假设每个批次的回归系数向量𝛃i\mathbf{\beta}_{i} ，与它们和所有其他批次的随机误差向量无关， 𝛆k\mathbf{\varepsilon}_{k} 。
RCR 建模包括三个步骤。
步骤a）对每个批次的数据拟合一个简单的线性回归模型
与批次相关的数据ii表示为(t伊,y伊)\left( t_{\text{ij}},y_{\text{ij}} \right) ， 在哪里y伊y_{\text{ij}}代表i日第i个批次在时间点的测量响应tjt_{j} ， 在哪里j=1,⋯,nij = 1,\ \cdots,\ n_{i} . 各个时间点tjt_{j}不同批次之间可能存在差异，时间点总数也可能存在差异nin_{i} . 设计矩阵𝐗i\mathbf{X}_{i}用于批量ii可以写成矩阵形式如下。
𝐗i=[1ti1⋮⋮1tini]\mathbf{X}_{i} = \begin{bmatrix} 1 & t_{i1} \\ \vdots & \vdots \\ 1 & t_{in_{i}} \\ \end{bmatrix}
该矩阵被归一化并指定为𝐌i\mathbf{M}_{i} 。
𝐌i=(𝐗i′𝐗i)−1\mathbf{M}_{i} = \left( \mathbf{X}_{i}^{'}\mathbf{X}_{i} \right)^{- 1}
截距的估计值βi1\beta_{i1} ，斜坡βi2\beta_{i2} ，回归的“均方误差” 微分方程i=σiy^2\text{MSE}_{i} = \sigma_{i\hat{y}}^{2}以及自由度自由度i=(ni−2)\text{df}_{i} = \left(n_{i} - 2 \right)与...相关微分方程i\text{MSE}_{i}是通过对每个批次的数据拟合简单的线性回归模型获得的。截距和斜率(𝛃1,𝛃2)′\left( \mathbf{\beta}_{1},\mathbf{\beta}_{2} \right)^{'}估计为
𝐛i=(𝐗i′𝐗i)−1𝐗i′和i\mathbf{b}_{i} = \left( \mathbf{X}_{i}^{'}\mathbf{X}_{i} \right)^{- 1}\mathbf{X}_{i}^{'}\mathbf{y}_{i}.
步骤 b）估计误差方差-协方差矩阵
估计随机系数的误差方差-协方差矩阵首先要估计合并的微分方程\text{均方误差}所有批次，即微分方程水池\text{MSE}_{\text{池}} ，如下所示。
微分方程水池=σ^2=∑i=1n自由度i•微分方程i∑i=1B自由度i\text{MSE}_{\text{pool}} = {\hat{\sigma}}^{2} = \frac{\sum_{i = 1}^{n}{\text{df}_{i} \bullet \text{MSE}_{i}}}{\sum_{i = 1}^{B}\text{df}_{i}}
该等式中的分子可以表示如下。
自由度i•微分方程i=上交所i=𝛆^i′𝛆^i=(和i−𝐗i𝐛i)′(和i−𝐗i𝐛i)\text{df}_{i} \bullet \text{MSE}_{i} = \text{SSE}_{i} = {\hat{\mathbf{\varepsilon}}}_{i}^{'}{\hat{\mathbf{\varepsilon}}}_{i} = \left( \mathbf{y}_{i} - \mathbf{X}_{i}\mathbf{b}_{i} \right)^{'}\left( \mathbf{y}_{i} - \mathbf{X}_{i}\mathbf{b}_{i} \right)
样本方差-协方差矩阵𝐒bb\mathbf{S}_{\text{bb}}估计截距和斜率(𝛃1,𝛃2)′\left( \mathbf{\beta}_{1},\mathbf{\beta}_{2} \right)^{'}计算如下。
𝐒bb=1n−1∑i=1n(𝐛i−β^超低密度脂蛋白)(𝐛i−β^超低密度脂蛋白)′\mathbf{S}_{\text{bb}} = \frac{1}{n - 1}\sum_{i = 1}^{n}\left( \mathbf{b}_{i} - {\hat{\beta}}_{\text{ULS}} \right)\left( \mathbf{b}_{i} - {\hat{\beta}}_{\text{ULS}} \right)^{'}
在哪里β^超低密度脂蛋白=1n∑i=1n𝐛i{\hat{\beta}}_{\text{ULS}} = \frac{1}{n}\sum_{i = 1}^{n}\mathbf{b}_{i} ，即系数的平均向量，即截距和斜率（ULS 代表“非加权最小二乘法”）。现在，方差-协方差矩阵𝚺^BB{\hat{\mathbf{\Sigma}}}_{\text{ββ}}可以如下估算。
𝚺^BB=𝐒bb−σ^2n•∑i=1n(𝐗i′𝐗i)−1=𝐒bb−σ^2•𝐌¯=[σβ12σβ1β2σβ1β2σβ22]{\hat{\mathbf{\Sigma}}}_{\text{ββ}} = \mathbf{S}_{\text{bb}} - \frac{{\hat{\sigma}}^{2}}{n} \bullet \sum_{i = 1}^{n}\left( \mathbf{X}_{i}^{'}\mathbf{X}_{i} \right)^{- 1} = \mathbf{S}_{\text{bb}} - {\hat{\sigma}}^{2} \bullet \overset{\overline{}}{\mathbf{M}} = \begin{bmatrix} \sigma_{\beta_{1}}^{2} & \sigma_{\beta_{1}\beta_{2}} \\ \sigma_{\beta_{1}\beta_{2}} & \sigma_{\beta_{2}}^{2} \\ \end{bmatrix}
在哪里𝐌¯=1n•∑i=1n𝐌i\overset{\overline{}}{\mathbf{M}} = \frac{1}{n} \bullet \sum_{i = 1}^{n}\mathbf{M}_{i} ，即均值模型方差-协方差矩阵。误差方差σ^2{\hat{\sigma}}^{2}估计分析方法的方差。截距方差分量σβ12\sigma_{\beta_{1}}^{2}表示产品放行时观察到的响应的批次内方差，这完全是由于制造过程造成的，即排除了由于分析方法造成的方差。斜率方差分量σβ22\sigma_{\beta_{2}}^{2}表示降解率的批次间差异。
由于方差-协方差矩阵𝚺^BB{\hat{\mathbf{\Sigma}}}_{\text{ββ}}是通过计算两个矩阵的差值来估计的，结果可能不是正定的，即截距和/或斜率方差分量可能为负。一种可能的补救方法是用 0 替换负值。这是负方差估计的一种简单急救方法，相当于将随机效应转换为固定效应。如果两个方差都为零，则 RCR 模型对应于普通最小二乘 (OLS) 模型。
步骤c）通过加权最小二乘（WLS）估计量估计系数均值向量
系数均值向量的估计首先计算𝐖i\mathbf{W}_{i} ，
𝐖i=(𝚺^BB+σ^2•𝐌i)−1\mathbf{W}_{i} = \left( {\hat{\mathbf{\Sigma}}}_{\text{ββ}} + {\hat{\sigma}}^{2} \bullet \mathbf{M}_{i} \right)^{- 1}
它对应于回归参数估计的方差的倒数i日第i个批次。基于此方程， 2×i2 \乘以 i截距和斜率系数估计的维矩阵，即加权最小二乘估计，现由以下方程得出：
𝛃^加权最小二乘法=(∑i=1n𝐖i)−1(∑i=1n𝐖i𝐛i)=𝛀(∑i=1n𝐖i𝐛i){\hat{\mathbf{\beta}}}_{\text{WLS}} = \left( \sum_{i = 1}^{n}\mathbf{W}_{i} \right)^{- 1}\left( \sum_{i = 1}^{n}\mathbf{W}_{i}\mathbf{b}_{i} \right) = \mathbf{\Omega}\left( \sum_{i = 1}^{n}\mathbf{W}_{i}\mathbf{b}_{i} \right)
在哪里𝛀\mathbf{\欧米茄}是加权最小二乘估计量的方差估计值。截距估计值的平均值表示在制造过程中产品发布时的响应平均值，斜率估计值的平均值表示共同斜率。
RCR 模型的公差区间 (TI)
RCR 模型系数β1\beta_{1}和β2\beta_{2}描述总体趋势，即历史批次的共同趋势。RCR 模型的 TI 可借助估计值计算y^j{\hat{y}}_{j}回归线及其标准误差σy^,RCR\sigma_{\hat{y},RCR}按照 Scholz 和 Vangel [5] 提出的方法，其方法利用了有效样本量的概念。
y^j±k•σy^,相对反应率{\hat{y}}_{j} \pm k \bullet \sigma_{\hat{y},\text{RCR}}
误差项σy^,RCR\sigma_{\hat{y},RCR}计算如下
σy^,RCR=D(𝐭(𝚺^BB+1n𝛀)𝐭′+σ^2)\sigma_{\hat{y},RCR} = D\left( \sqrt{\mathbf{t}\left( {\hat{\mathbf{\Sigma}}}_{\text{ββ}} + \frac{1}{n}\mathbf{\Omega} \right)\mathbf{t}^{\mathbf{'}}\mathbf{+}{\hat{\sigma}}^{2}} \right)
在哪里𝐭=[1t1⋮⋮1tk],kϵ N\mathbf{t} = \begin{bmatrix}
1 & t_{1} \\
\vdots & \vdots \\
1 & t_{k} \\
\end{bmatrix},\ k\ \epsilon\mathbb{\text{\ N}} ， 和DD是参数的对角线。因子kk给出的是
k=n−1n•1n*−1•tn*−1,zp•n*,γk = \sqrt{\frac{n - 1}{n}} \bullet \frac{1}{\sqrt{n^{*} - 1}} \bullet t_{n^{*} - 1,\ z_{p} \bullet \sqrt{n^{*}},\gamma}
其中因素tn*−1,zp•n*,γt_{n^{*} - 1,z_{p} \bullet \sqrt{n^{*}},\gamma}是γ\伽马 t分布的分位数n*−1n^{*} - 1自由度和非中心参数zp•n*z_{p} \bullet \sqrt{n^{*}} ， 在哪里zpz_ {p}是pp标准正态分布的分位数。参数n*n^{*}有效样本量由以下公式给出
n*=[σb2σb2+σe2•∑i=1nb(nin)2+1n•σe2σb2+σe2]−1=[ρ•1f+1+(1−ρ)•1n]−1{n^{*} = \left\lbrack \frac{\sigma_{b}^{2}}{\sigma_{b}^{2} + \sigma_{e}^{2}} \bullet \sum_{i = 1}^{n_{b}}\left( \frac{n_{i}}{n} \right)^{2} + \frac{1}{n} \bullet \frac{\sigma_{e}^{2}}{\sigma_{b}^{2} + \sigma_{e}^{2}} \right\rbrack}^{- 1} = \left\lbrack \rho \bullet \frac{1}{f + 1} + \left( 1 - \rho \right) \bullet \frac{1}{n} \right\rbrack^{- 1}
其中参数ff和方差比ρ\rho定义为
| f=1∑i=1nb(nin)2−1f = \frac{1}{\sum_{i = 1}^{n_{b}}\left( \frac{n_{i}}{n} \right)^{2}} - 1 |
ρ=σb2σb2+σe2\rho = \frac{\sigma_{b}^{2}}{\sigma_{b}^{2} + \sigma_{e}^{2}} |
在哪里σb2\sigma_{b}^{2}是批次间方差， σe2\sigma_{e}^{2}是批次内方差（或方法方差）， nin_{i}是每批测量次数ii ， IE i=1,⋯,nbi = 1,\ \cdots,\ n_{b} ， nbn_{b}是历史批次的数量， nn是样本总数，即n=∑i=1nbnin = \sum_{i = 1}^{n_{b}}n_{i} 。
TI 涵盖p(100%)p\left( 100\% \right)人口的概率γ(100%)\gamma\left( 100\% \right)有效样本量概念考虑到批次内数据比批次间数据更具相关性，从而减少了独立信息量。因此，以下情况始终正确： n*≤nn^{*} \leq n . 理想情况n*=nn^{*} = n如果批次内数据完全不相关，则可以实现。
JMP 中的 RCC/RCR 实现
实现 RCC/RCR 程序的 JSL 脚本分为几个部分。
首先必须定义参考和测试批次名称以及响应、时间、批次和存储条件的各种变量名称。
RCR 模型的构建（简化的 Carter 和 Yang [1]）如上所述。
RCR 模型的 TI 是使用 Scholz 和 Vangel [5] 提出的方法构建的。
生成一个报告窗口，以表格和图形形式总结计算出的参数。
例子
为了评估脚本功能，我们为 9 个历史批次和一个当前批次生成了各种数据集。以下代码片段执行 RCR 模型拟合的步骤a )，即截距 (a、 𝛃1\mathbf{\beta}_{1} ），斜率（b， 𝛃2\mathbf{\beta}_{2} )，(根)均方误差（(RMSE ) MSE， 微分方程i\text{MSE}_{i} )，自由度（df， 自由度i\text{df}_{i} ) 和观测值数量 (n, nin_{i} ） 的ii收集历史批次并将其存储到长度为nbn_{b} 。
For( i = 1, i <= n.refBatches, i++,
fm = Fit Model(
Y( :Eval( responseVariable ) ),
Effects( :Month ),
Personality( "Standard Least Squares" ),
Emphasis( "Minimal Report" ),
Where( :Batch == refBatches[i] ),
Invisible
);
fit = fm << Run;
fitr = fit << Report;
a.i[i, 1] = fitr["Parameter Estimates"][Number Col Box( 1 )] << Get( 1 );
b.i[i, 1] = fitr["Parameter Estimates"][Number Col Box( 1 )] << Get( 2 );
RMSE.i[i, 1] = fitr["Summary of Fit"][Number Col Box( 1 )] << Get( 3 );
df.i[i, 1] = fitr["Analysis of Variance"][Number Col Box( 1 )] << Get( 2 );
n.i[i, 1] = fitr["Summary of Fit"][Number Col Box( 1 )] << Get( 5 );
);
For( i = 1, i <= n.refBatches, i++,
MSE.i[i] = RMSE.i[i] ^ 2
);
步骤b ) 从估计池化 MSE (MSE.pool， 微分方程水池\text{MSE}_{\text{池}} ）随后，平均模型方差-协方差矩阵（M.bar， 𝐌¯\overset{\overline{}}{\mathbf{M}} ) 进行估计，然后将其与样本方差-协方差矩阵 (S, 𝐒bb\mathbf{S}_{\text{bb}} ）得到方差-协方差矩阵（Sigma.hat， 𝚺^BB{\hat{\mathbf{\Sigma}}}_{\text{ββ}} ）。
ab.i = a.i || b.i; // Matrix of the intercept and slope vectors
S = Covariance( ab.i ); // Sample variance-covariance matrix of the intercepts and slopes
MSE.pool = Sum( E Mult( df.i, MSE.i ) ) / Sum( df.i ); // Pooled MSE
RMSE.pool = Sqrt( MSE.pool );
M.i = J( 2, 2, 0 );
For( i = 1, i <= n.refBatches, i++,
dat << Select All Rows << Invert Row Selection;
dat << Select Where( :Batch == refBatches[i] );
x = dat << Get As Matrix( {"Month"} );
x.length = N Row( x ); // Number of observations
colReprIntcpt = J( x.length, 1, 1 ); // Column of 1 representing the intercept
x = colReprIntcpt || x; // Design matrix of the OLS model
M.i = M.i + Inverse( x` * x );
);
M.bar = M.i / n.refBatches; // Mean model variance-covariance matrix
Sigma.hat = S - MSE.pool * M.bar;
Sigma.hat[Loc( Diag( Sigma.hat ) < 0 )] = 0; // Replace negative values on the diagonal of Sigma.hat by 0.
In step c) the RCR model coefficients ab.hat (𝛃^加权最小二乘法{\hat{\mathbf{\beta}}}_{\text{WLS}}) are estimated by aid of the following code snippet.
M.i = J( 2, 2, 0 );
W.i = J( 2, 2, 0 );
For( i = 1, i <= n.refBatches, i++,
dat << Select All Rows << Invert Row Selection;
dat << Select Where( :Batch == refBatches[i] );
x = dat << Get As Matrix( {"Month"} );
colReprIntcpt = J( n.i[i, 1], 1, 1 ); // Column of 1 representing the intercepts
x = colReprIntcpt || x; // Design matrix of the linear model
M.i = M.i + Inverse( x` * x );
W.i = W.i + Inverse( Sigma.hat + MSE.pool * M.i );
);
Omega = Inverse( W.i );
ab.hat = Omega * (W.i * ab.i`);
ab.hat = ab.hat`;
ab.hat.bar = V Mean( ab.hat );
下一个代码片段显示了根据 Scholz und Vangel [5] 构建 RCR 模型的 TI。为了构建区间，使用了存储在 ab.hat.bar 向量中的平均截距和斜率系数估计值。
utp = Associative Array( dat:Month ); // Unique time points
time.points = Matrix( utp << Get Keys ); // Vector with (new) time points for which the response is sought
n.new = N Row( time.points ); // Number of data points of the new data set
m.x = J( n.new, 1, 1 ) || time.points; // Design matrix
n = Sum( n.i ); // Total number of data points
df.th = n - 2; // Number of degrees of freedom if only fixed effects are allowed, i.e. if a "common slope, common intercept" model is appropriate.
df = df.th;
variance.ratio = Sigma.hat[1, 1] / (Sigma.hat[1, 1] + MSE.pool);
t.f = 1 / Sum( (df.i / n) ^ 2 ) - 1;
If( Sigma.hat[1, 1] == 0,
df.eff = df,
df.eff = 1 / (variance.ratio / (t.f + 1) + (1 - variance.ratio) / n)
);
// Calculation of the k.eff factor on the basis of the effective number of degrees of freedom.
k.eff = Sqrt( (n - 1) / n ) * 1 / Sqrt( df.eff - 1 ) * t Quantile(
gamma,
(df.eff - 1),
(Normal Quantile( (1 - alpha) ) * Sqrt( df.eff ))
);
t.quantil = t Quantile( (1 - alpha / 2), df.eff );
y.hat = ab.hat.bar[1] + ab.hat.bar[2] * time.points;
LTL = y.hat - k.eff * Vec Diag( Sqrt( m.x * (Sigma.hat + Omega / n.refBatches) * m.x` + MSE.pool ) );
UTL = y.hat + k.eff * Vec Diag( Sqrt( m.x * (Sigma.hat + Omega / n.refBatches) * m.x` + MSE.pool ) );
在脚本的最后一部分，结果准备显示在报告窗口中。示例报告如图 2 所示。
报告窗口总结的结果可以回答第 3 章开头提出的以下问题：1) 由于当前批次的 18 个月数据点（即最近的测量值）位于与历史批次数据相适应的 TI 之外，因此必须得出结论，当前批次为 OOE/OOT，以及 2) 由于当前批次的 18 个月数据点位于当前批次的 PI 之外，因此必须得出结论，就当前批次而言，当前批次也是 OOE/OOT。尽管 12 个月时间点位于当前批次的 PI 范围内，但在前一个时间站也应该得出相同的结论。但是，由于该点位于与历史批次数据相适应的 TI 之外，因此 12 个月时间站的总体结论一定是这是一个 OOE/OOT 点。
如第 3 章开头所述，各个历史批次的截距和斜率很可能存在差异，因此 RCR 模型优于普通最小二乘模型。可以区分五种情况：情况 I：无截距变化，无斜率变化；情况 II：截距变化，无斜率变化；情况 III：无截距变化，斜率变化；情况 IV：截距变化，斜率变化；情况 V：非线性。图 3 总结了所有五种情况的结果。图 3 中的图表显示，在截距和/或斜率变化的情况下，RCR 模型可以更好地表示数据，而在截距和斜率没有变化的情况下，RCR 模型与基本方法非常相似。
Figure_2图 2图 2
图 2 实现 RCC/RCR 方法的脚本输出
窗口分为“历史批次评估”、“当前批次评估（单个模型参数）”和汇总历史和当前批次数据的图形显示。历史批次的数据用蓝点（）表示，当前批次的数据用红叉（）表示。历史批次数据的拟合 TI 由蓝色虚线划定的带表示，当前批次数据的拟合 PI（但不包括最新观察值）由红色虚线划定的带表示。
图3 不同情况下RCR模型区间与基本线性回归模型区间对比
RCR 模型拟合度用蓝色实线回归线和相应公差区间（95% 置信度和 95% 覆盖度）的蓝色虚线表示。基本模型拟合度用紫色虚线回归线和相应预测区间（95% 置信度）的紫色阴影表示，即单个预测值的置信限度。面板 Va 中的灰色虚线表示基于平方响应计算的非线性拟合度。六个面板比较了以下情况下的 RCR 模型区间和基本线性回归模型区间：面板 I：“无截距变化，无斜率变化”情况；面板 II：“截距变化，无斜率变化”情况；面板 III：“无截距变化，斜率变化”情况；面板 IV：“截距变化，斜率变化”情况；面板 Va：“非线性”情况；面板 Vb：“非线性”情况，即响应变量变换（即平方）后。
结论
本文介绍的 JMP 脚本是一种合适且强大的工具，可用于评估药品当前批次的 OOE 和 OOT 结果的稳定性数据。前提条件是至少有三批药品的历史稳定性数据，每批药品至少有四个时间点。
参考书目
Carter, RL 和 Yang, MCK 随机系数回归模型中的大样本推断。 《统计通讯 - 理论与方法》 （1986 年） 15 (8): 2507-2525。 链接
Conover, WJ (1980)，《实用非参数统计学》，纽约：John Wiley and Sons, Inc.
ICH 指南 Q1A(R2)，新药物物质和产品的稳定性测试，第 4 步版本，2003 年 2 月 6 日。 链接
Mandel, BJ 回归控制图。J Qual Technol (1969) 1 (1): 1-9。 链接
Scholz, F. 和 Vangel, M. 批次效应下的公差界限和 Cpk 置信界限。发表于：可靠性、质量和安全随机模型高级版。由 Kahle, W.、von Collani, E.、Franz, J. 和 Jensen, U. Birkhäuser 编辑，波士顿，1998 年，第 361-379 页。 链接
Shewhart, WA (1939)。从质量控制角度看统计方法。（Deming，WE）。华盛顿特区，农业部研究生院。
标准操作程序 – 实验室数据管理指南 – 超出预期 (OOE) 和超出趋势 (OOT) 结果。欧洲合规学院 (ECA) 分析质量控制工作组，Christopher Burgess，草案 V 0.3，2015 年。
英国药品和保健产品管理局 (MHRA)，《不合格调查》，演讲笔记和代表资料包，2013 年 8 月 28 日。 链接
世卫组织技术报告系列第 953 号附件 2 活性药物成分和成品药物的稳定性测试，2009 年。 链接
