code{white-space: pre;}
Pius Dahinden 博士，Tillotts Pharma 分析科学经理
抽象的
自新的欧盟 GMP 指南（EudraLex 第 4 卷）第 6 章（质量控制）于 2014 年生效以来，对药品稳定性数据进行趋势分析是一项额外的分析质量控制要求。欧洲合规学院 (ECA) 分析质量控制工作组 (AQCWG) 制定了用于评估超出预期 (OOE) 和超出趋势 (OOT) 结果的标准操作程序。为了评估稳定性数据，建议使用回归控制图（RCC）程序将当前批次与历史批次的数据进行比较。然而，根据历史批次数据本身估计的斜率和/或截距具有不确定性。 ECA AQCWG 提出了一种考虑到这一点的简化随机系数回归 (RCR) 方法。提出了一个新开发的 JMP 脚本，该脚本实现了推荐的 RCC / RCR 程序。该脚本使稳定性数据中 OOE / OOT 点的识别变得简单，并在输出中提供所有必要的信息以及适当的图形显示。
介绍
制造商应用统计过程控制（SPC）来确保生产产品的质量并防止因生产过程失控而造成的财务损失。过程控制意味着对正在运行的过程进行例行监控，如果任何监控参数出现异常趋势，则采取措施使过程恢复受控。 SPC 是一种工具，使制造商能够及早检测流程中的变化，从而使他们能够调查引起问题的问题并在早期阶段采取纠正措施。过程监控的一个强大工具是由 Shewhart [6] 推广的控制图方法。图 1 显示了典型控制图的示意图。
与不应该有趋势的生产过程相反，存在许多具有固有趋势的数据管理问题。为了分析此类数据管理问题，传统的控制图方法与回归分析相结合，这种方法被 Mandel 称为“回归控制图”（RCC）[4]。由于控制限制了公差区间 (TI) 的使用，因此是合理的。为了监测在推荐储存条件（例如 25 °C / 60% RH）下储存的药品的稳定性数据，RCC 用于监测活性药物成分随时间的降解或杂质的积累，以支持产品的有效期，即确认药品在标签保质期内的稳定性。
Figure_1图1
图1 典型控制图组成
控制图通常显示值X X或样本质量特性的几个值（例如平均值）相对于样本数量或时间的统计。一般来说，图表包含一条代表平均值的中心线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 的评估，至少需要四个数据点，因为最近测量的数据点不应用于评估问题 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模型可以写为
𝐲ij=𝛃1i+𝛃2i•𝐭ij+𝛆ij\mathbf{y}_{\text{ij}} = \mathbf{\beta}_{1i} + \mathbf{\beta}_{2i} \bullet \mathbf{t}_{\text{ij}} + \mathbf{\varepsilon}_{\text{ij}}
在哪里 𝐲我\mathbf{y}_{i}是一个n我× 1 n_{i} \乘以1向量表示我th我^{\text{第}}批次的n我你}重复测量（其中j = 1 , ⋯ , n我j = 1,\ \cdots,\ n_{i} ， 和n我你}代表批次重复测量的总数我我) 响应变量（例如含量 [%]、杂质总和 [%] 或溶出度 [%]）， 𝛃 1我\mathbf{\beta}_{1i}和 𝛃 2我\mathbf{\beta}_{2i}分别是批次特定截距和斜率的向量， 𝐭我\mathbf{t}_{i}是个我th我^{\text{第}}批次的时间向量（例如[月]），以及 𝛆我\mathbf{\varepsilon}_{i}是一个n我× 1 n_{i} \乘以1的向量n我你}随机错误。该 RCR 模型的先决条件是响应与解释变量之间存在线性关系。如果关系不是线性的，则必须在执行 RCR 建模之前应用适当的变换。
假设每个ε ij \varepsilon_{\text{ij}}来自均值为零、方差为零的单变量正态总体σ 2 \西格玛^{2} ， 和所有ε ij \varepsilon_{\text{ij}}彼此独立（条件独立）并且与模型系数独立β 1我\beta_{1i}和β 2我\beta_{2i} 。这一假设是合理的，因为分析方法的变异性极不可能取决于正在测试的批次或执行测试的时间点。进一步假设每个批次的回归系数向量， 𝛃我\mathbf{\beta}_{i} ，独立于它们和所有其他批次的随机误差向量， 𝛆 k \mathbf{\varepsilon}_{k} 。
RCR 建模包括三个步骤。
步骤 a) 对每个批次的数据拟合一个简单的线性回归模型
与批次相关的数据我我表示为 ( t ij , y ij ） \left( t_{\text{ij}},y_{\text{ij}} \right) ， 在哪里y ij y_{\text{ij}}代表我th我^{\text{第}}批次在时间点测量的响应t j t_{j} ， 在哪里j = 1 , ⋯ , n我j = 1,\ \cdots,\ n_{i} 。个别时间点t j t_{j}批次之间可能有所不同，时间点的总数也可能有所不同n我你} 。设计矩阵 𝐗我\mathbf{X}_{i}对于批次我我可以写成矩阵形式如下。
𝐗i=[1ti1⋮⋮1tini]\mathbf{X}_{i} = \begin{bmatrix}1&t_{i1}\\\vdots&\vdots\\1&t_{in_{i}}\\\end{bmatrix}
该矩阵被归一化并指定为 𝐌我\mathbf{M}_{i} 。
𝐌i=(𝐗i′𝐗i)−1\mathbf{M}_{i} = \left( \mathbf{X}_{i}^{'}\mathbf{X}_{i} \right)^{- 1}
截距的估计β我1 \beta_{i1} ，斜坡β我2 \beta_{i2} ，回归的“均方误差” 均方误差我= σ我y ^ 2 \text{MSE}_{i} = \sigma_{i\hat{y}}^{2} ，以及自由度 df我= ( n我- 2 ） \text{df}_{i} = \left( n_{i} - 2 \right)有关联均方误差我\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{MSE}所有批次中，即均方误差水池\text{MSE}_{\text{池}} ， 如下。
MSEpool=σ^2=∑i=1ndfi•MSEi∑i=1Bdfi\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}}
该方程中的分子可以表示如下。
dfi•MSEi=上证所i=𝛆^i′𝛆^i=(𝐲i−𝐗i𝐛i)′(𝐲i−𝐗i𝐛i)\text{df}_{i} \bullet \text{MSE}_{i} = \text{SSE}_{i} = {\hat{\mathbf{\valuepsilon}}}_{i}^{' }{\hat{\mathbf{\varepsilon}}}_{i} = \left( \mathbf{y}_{i} - \mathbf{X}_{i}\mathbf{b}_{i}\右)^{'}\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} - {\帽子{\beta}}_{\text{ULS}} \right)\left( \mathbf{b}_{i} - {\hat{\beta}}_{\text{ULS}} \right)^ {'}
where β^超大规模激光雷达=1n∑i=1n𝐛i{\hat{\beta}}_{\text{ULS}} = \frac{1}{n}\sum_{i = 1}^{n}\mathbf{b}_{i}, i.e. the mean vector of the coefficients, i.e. intercepts and slopes (ULS stands for "Unweighted Least-Squares"). Now, the variance-covariance matrix 𝚺^BB{\hat{\mathbf{\Sigma}}}_{\text{ββ}} can be estimated as follows.
𝚺^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 \translated{\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}
在哪里 𝐌 ￣ = 1 n • Σ我= 1 n 𝐌我\overset{\overline{}}{\mathbf{M}} = \frac{1}{n} \bullet \sum_{i = 1}^{n}\mathbf{M}_{i} ，即平均模型方差-协方差矩阵。误差方差σ ^ 2 {\帽子{\西格玛}}^{2}估计分析方法的方差。截距方差分量σ β 1 2 \sigma_{\beta_{1}}^{2}表示产品放行时观察到的响应的批次内差异，该差异仅由制造过程引起，即排除由分析方法引起的差异。斜率方差分量σ β 2 2 \sigma_{\beta_{2}}^{2}表示降解率的批次间方差。
Since the variance-covariance matrix 𝚺^BB{\hat{\mathbf{\Sigma}}}_{\text{ββ}} is estimated calculating the difference of two matrices, the result may not be positive-definite, i.e. the intercept and / or the slope variance component(s) may be negative. One possible remedy to this is replacing the negative value(s) with 0. This is a simple first-aid approach for negative variance estimates and is equivalent to converting a random effect into a fixed effect. In case that both variances are zero the RCR model corresponds to the ordinary least squares (OLS) model.
步骤 c) 通过加权最小二乘 (WLS) 估计器估计系数均值向量
系数均值向量的估计从计算估计值开始 𝐖我\mathbf{W}_{i} ,
𝐖i=(𝚺^BB+σ^2•𝐌i)−1\mathbf{W}_{i} = \left( {\hat{\mathbf{\sigma}}}_{\text{ββ}} + {\hat{\sigma}}^{2} \bullet \mathbf {M}_{i}\右)^{-1}
它对应于回归参数估计的方差的倒数我th我^{\text{第}}批。基于这个方程2 ×我2 \乘以我截距和斜率系数估计的维矩阵，即加权最小二乘估计，现在通过以下等式获得：
𝛃^世界LS=(∑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 {\帽子{y}}_{j}回归线及其标准误差σ y ^ ,右C右\sigma_{\hat{y},RCR}根据Scholz和Vangel[5]提出的方法如下。他们的方法利用了有效样本量的概念。
y^j±k•σy^,RCR{\hat{y}}_{j} \pm k \bullet \sigma_{\hat{y},\text{RCR}}
误差项σ y ^ ,右C右\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)
在哪里 H = [ 1 t 1 ⋮ ⋮ 1 t k ] ， 1999 。 k ε 经过\mathbf{t} = \begin{bmatrix}1 & t_{1}\\\vdots&\vdots\\1 & t_{k}\\end{bmatrix},\k\ \epsilon\mathbb{\text {\经过}} ， 和D D是参数的对角线。因素k k是（谁）给的
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}
其中因子t n * - 1 , z p • n * , γ t_{n^{*} - 1,z_{p} \bullet \sqrt{n^{*}},\gamma}是个γ \伽玛t分布的分位数n * - 1 n^{*} - 1自由度和非中心参数z p • n * z_{p} \bullet \sqrt{n^{*}} ， 在哪里z p z_{p}是个p p标准正态分布的分位数。参数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}
其中参数F F和方差比ρ \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}} |
在哪里σ乙2 \sigma_{b}^{2}是批次间方差， σ e 2 \sigma_{e}^{2}是批次内方差（或方法方差）， n我你}是每批次的测量次数我我， IE我= 1 , ⋯ , n乙i = 1,\ \cdots,\ n_{b} , n乙n_{b}是历史批次的数量，并且n n是样本总数，即n = Σ我= 1 n乙n我n = \sum_{i = 1}^{n_{b}}n_{i} 。
TI 涵盖p ( 100 % ） p\左( 100\% \右)人口的概率为γ ( 100 % ） \伽玛\左(100\%\右) 。有效样本量概念考虑到批次内数据比批次间数据更具相关性，从而减少了独立信息量。因此，以下内容始终正确： n * ≤ n n^{*} \leq n 。理想的情况n * = n n^{*} = n如果批次内数据完全不相关，则可以实现。
在 JMP 中实施 RCC / RCR
实现 RCC/RCR 过程的 JSL 脚本分为几个部分。
首先，必须定义参考和测试批次名称以及响应的各种变量名称、时间、批次和存储条件。
RCR 模型的构建（简化的 Carter 和 Yang [1]）如上所述。
RCR 模型的 TI 是使用 Scholz 和 Vangel [5] 提出的方法构建的。
生成一个报告窗口，以表格和图形形式总结计算出的参数。
例子
为了评估脚本功能，为九个历史批次和一个当前批次生成了各种数据集。下面的代码片段执行RCR模型拟合的步骤a )，即截距(a, 𝛃 1 \mathbf{\beta}_{1} ），斜坡（b， 𝛃 2 \mathbf{\beta}_{2} ), (根) 均方误差 ((RMSE ) MSE, 均方误差我\text{MSE}_{i} ), 自由度 (df, df我\text{df}_{i} ) 和观测值数量 (n, n我你} ） 的我我收集历史批次并将其存储到长度向量中n乙n_{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
);
Step b) starts with estimating the pooled MSE (MSE.pool, MSEpool\text{MSE}_{\text{pool}}). Subsequently, the mean model variance-covariance matrix (M.bar, 𝐌¯\overset{\overline{}}{\mathbf{M}}) is estimated, which is then combined with the sample variance-covariance matrix (S, 𝐒bb\mathbf{S}_{\text{bb}}) to obtain the variance-covariance matrix (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 (𝛃^世界LS{\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 模型优于普通最小二乘模型。可以区分五种情况： 情况一：没有截距变化，没有斜率变化；情况二：截距变化，无斜率变化；情况III：没有截距变化、斜率变化；情况IV：截距变化、斜率变化；情况V：非线性。图 3 总结了所有五种情况的结果。图 3 中的图表表明，RCR 模型更好地表示了截距和/或斜率变化的情况下的数据，并且与截距和斜率没有变化的情况下的基本方法非常相似。
Figure_2图_2
图2 实现RCC/RCR方法的脚本输出
该窗口分为“历史批次评估”、“当前批次评估（单个模型参数）”部分以及总结历史和当前批次数据的图形显示。历史批次的数据用蓝点 () 表示，当前批次的数据用红叉 () 表示。适合历史批次数据的 TI 由蓝色虚线界定的带表示，适合当前批次数据（但不包括最近观察）的 PI 由红色虚线界定的带表示。
图3 不同情况下RCR模型区间与基本线性回归模型区间对比
RCR 模型拟合由蓝色实线回归线和相应公差区间的蓝色虚线表示（95% 置信度和 95% 覆盖率）。基本模型拟合由紫色虚线回归线和相应预测区间（95% 置信度）的紫色阴影表示，即单个预测值的置信限。 Va 面板中的灰色虚线表示基于平方响应计算的非线性拟合。六个面板比较了以下情况下的 RCR 模型区间和基本线性回归模型区间： 面板 I：“无截距变化，无斜率变化”情况；第二组：“截距变化，无斜率变化”情况；第三组：“无截距变化，斜率变化”案例；第四组：“截距变化、斜率变化”案例；面板 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 回归控制图。质量技术杂志（1969） 1 （1）：1-9。 关联
Scholz, F. 和 Vangel, M. 批次效应下的容差界限和 Cpk 置信界限。发表于：可靠性、质量和安全随机模型高级。由 Kahle, W.、von Collani, E.、Franz, J. 和 Jensen, U. Birkhäuser 编辑，波士顿，1998 年，361-379。 关联
华盛顿州休哈特 (1939)。从质量控制的角度来看的统计方法。 （戴明，我们）。华盛顿特区，农业部研究生院。
标准操作程序 – 实验室数据管理指南 – 超出预期 (OOE) 和超出趋势 (OOT) 结果。欧洲合规学院 (ECA) 分析质量控制工作组，Christopher Burgess，草案 V 0.3，2015 年。
英国药品和保健品监管局 (MHRA)，不合格调查、演示说明和代表包，2013 年 8 月 28 日。链接
世界卫生组织技术报告系列，第 953 号，附件 2 活性药物成分和成品药物的稳定性测试，2009 年。 链接
