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}控制過程。一個基本假設是資料呈常態分佈。另外兩條特徵水平線是距離中心線 ± 3 個標準差 (SD) 的下限和上限控制限 (分別為 LCL 和 UCL)，只要過程處於控制之中，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{MSE}_{i}}}{\sum_{i = 1}^{B}}
此等式中的分子可以表示如下。
自由度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)
樣本變異數-協方差矩陣𝐒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)\if^\math}\p;
在哪裡β^超低密度脂蛋白=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( \right} {S}_{\text{bb}} - {\hat{\sigma}}^{2} \bullet \overset{\overline{}}{\mathbf{M}} = \begin{bmatrix} \sigma_{\beta_{1}}^{2} & \sigma_{\beta_{1}\beta_{\beta_{1 \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}\mathf{Wright}\ = 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}\math」{\Omega} \right)\mathbf; 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( ma_{e}^{2}}{\sigma_{b}^{2} + \sigma_{e}^{2}} \right\rbrack}^{- 1} = \left\lbrack \rho \bullet \frac{1}{f + 1} + \iftrack 1 - \rho \bullet \frac{1}{f + 1} + \rbC 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] 提出的方法建構的。
產生一個報告窗口，以表格和圖形形式總結計算出的參數。
例子
為了評估腳本的功能，我們針對九個歷史批次和一個目前批次產生了各種資料集。以下程式碼片段執行 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；2) 由於當前批次的 18 個月數據點位於當前批次的 PI 批次為 OOE/OOT；2) 由於當前批次的 18 個月數據點位於當前批次的 PI 批次為 OOE 批次；雖然 12 個月的時間點落在目前批次的 PI 範圍內，但在前一個時間站也應該已經得出了同樣的結論。然而，由於該點位於與歷史批次資料相適應的 TI 之外，因此 12 個月時間站的總體結論必須是這是一個 OOE/OOT 點。
如第 3 章開頭所提到的，各個歷史批次的截距和斜率很可能有所不同，因此 RCR 模型優於一般最小平方法模型。可區分出五種情況：情況 I：截距無變化，斜率無變化；情況 II：截距有變化，斜率無變化；情況三：無截距變化，斜率變化；案例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：「無截距變化，斜率變化」情況；第四面板：「截距變化、斜率變化」案例；面板 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。 關聯
休哈特，WA（1939 年）。從品質控制角度的統計方法。 （戴明，WE）。華盛頓特區，農業部研究生院。
標準作業程序 – 實驗室資料管理指南 – 超出預期（OOE）和超出趨勢（OOT）的結果。歐洲合規學院 (ECA) 分析品質控制工作小組，Christopher Burgess，草案 V 0.3，2015 年。
英國藥品和保健產品管理局 (MHRA) ， 《不合格調查》，演講筆記和代表資料包，2013 年 8 月 28 日。
世衛組織技術報告系列第953號附件 2 活性藥物成分和成品藥物的穩定性測試，2009 年。
