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 年。 鏈接
