code{white-space: pre;}
ピウス・ダヒンデン博士、ティロッツ・ファーマ、分析科学マネージャー
抽象的な
新しい EU GMP ガイドライン (EudraLex 第 4 巻) の第 6 章 (品質管理) が 2014 年に発効して以来、医薬品の安定性データに対する傾向分析の実行は、分析品質管理の追加要件となっています。欧州コンプライアンス アカデミー (ECA) の分析品質管理ワーキング グループ (AQCWG) は、予想外 (OOE) および傾向外 (OOT) の結果を評価するための標準操作手順を策定しました。安定性データの評価には、現在のバッチと過去のバッチのデータの比較に回帰管理図 (RCC) 手順を使用することが推奨されます。ただし、過去のバッチ データ自体から推定される傾きや切片には不確実性が伴います。ECA AQCWG は、この点を考慮した簡略化されたランダム係数回帰 (RCR) 法を提案しています。推奨される RCC / RCR 手順を実装する、新しく開発された JMP スクリプトを紹介します。このスクリプトにより、安定性データ内の OOE/OOT ポイントの識別が簡単になり、適切なグラフィック表示とともに必要なすべての情報が出力に提供されます。
導入
製造業者は、生産品の品質を確保し、生産プロセスの制御不能によって生じる経済的損失を防止するため、統計的工程管理（SPC）を適用しています。工程管理とは、実行中の工程を日常的に監視し、監視対象のパラメータのいずれかに異常な傾向が見られる場合には、工程を再び制御下に戻すための措置を講じることを意味します。SPCは、製造業者が工程の変化を早期に検出し、問題の原因を調査して早期に是正措置を講じることを可能にするツールです。工程監視の強力なツールとして、シューハート[6]によって普及された管理図法があります。図1は、典型的な管理図の模式図を示しています。
傾向が存在しないはずの製造プロセスとは対照的に、データ管理の問題には固有の傾向が伴うものが数多く存在します。こうしたデータ管理の問題の分析には、従来の管理図法と回帰分析が組み合わされており、この手法はMandel [4]によって「回帰管理図」（RCC）と呼ばれています。管理限界として、許容区間（TI）の利用は合理的です。推奨保管条件（例：25℃ / 60% RH）で保管された医薬品の安定性データをモニタリングするために、RCCは、例えば、有効成分の劣化や不純物の経時的蓄積をモニタリングし、製品の有効期限設定をサポートするために、すなわち、ラベルに記載されている有効期間にわたって医薬品の安定性を確認するために使用されます。
Figure_1図_1図_1
図1 典型的な管理図の構成
管理図は通常、値を表示しますXXまたは、サンプルの品質特性の複数の値（例えば平均値）とサンプル数または時間との関係を示す統計。一般的に、チャートには平均値を表す中心線が含まれる。 X¯\overset{\overline{}}{X}管理下の工程の。基本的な前提は、データが正規分布しているということである。追加の 2 つの特徴的な水平線は、中心線から標準偏差 (SD) ± 3 の距離にある下側管理限界 (LCL) と上側管理限界 (UCL) であり、工程が管理されている限り、共通原因変動と呼ばれるランダム変動により、データ ポイントの 99.7% がこれらの限界内に収まる。中心線から SD ± 2 の距離にある線は、下側警告限界 (LWL) と上側警告限界 (UWL) と呼ばれ、中心線から SD ± 1 の距離にある線には特定の名前はありません。これらの線を利用して、管理図をゾーン C から A に分割できます。工程が管理されていないことを示す特殊原因変動の検出ルールでは、これらのゾーンでの観測頻度を判断基準として使用します。
手順と結果
医薬品の安定性データに適用される RCC は、次の 2 つの質問に答えることを目的としています。1) 現在のバッチは、過去のバッチと比較して期待外 (OOE) または傾向外 (OOT) ですか (つまり、過去のデータから得られた TI に関して)。2) 現在のバッチの最も最近測定されたデータ ポイントは、現在のバッチと比較して OOE ですか (つまり、その予測区間 (PI) に関して)。
質問 1 の観点で傾向を評価するには、現在のバッチの少なくとも 3 つのデータ ポイントが必要です。質問 2 の評価には、最後に測定されたデータ ポイントを PI の構築に使用しない必要があるため、少なくとも 4 つのデータ ポイントが必要です。
個々のバッチの切片と傾きは変動する可能性が高いため、通常の最小二乗モデル（共通傾き-共通切片）は意味をなさない可能性があります。このため、欧州コンプライアンスアカデミー（ECA）分析品質管理ワーキンググループ（AQCWG）は、個々のバッチの切片と傾きの変動を考慮した「ランダム係数回帰」（RCR）法の使用を推奨しています[7]。簡略化されたRCRモデル（Carter and Yang [1]の修正モデル）には、切片と傾きの2つの未知パラメータが含まれます。 β1\beta_{1}そして斜面β2\beta_{2}この手順はWHO [9]およびICH [2]のガイドラインに準拠しています。
OOTという用語は、英国医薬品医療製品規制庁（MHRA）によって、他の安定性バッチと比較しても、または安定性研究中に収集された以前の結果に関しても、予想される傾向に従わない安定性試験結果と定義されています[8]。OOSという用語は、MHRAによって、事前に決定された許容基準に準拠しない試験結果、または公式の公定書および/または企業の文書によって確立された許容基準から外れた試験結果と定義されています[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}}
どこそしてi\mathbf{y}_{i}はni×1n_{i} \times 1ベクトルを表すi番目i^{\text{th}}バッチの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^{\text{th}}バッチの時間ベクトル（例：[月]）、および𝛆i\mathbf{\varepsilon}_{i}はni×1n_{i} \times 1ベクトルnin_{i}ランダム誤差。このRCRモデルの前提条件は、応答変数と説明変数の間に線形関係があることです。関係が線形でない場合は、RCRモデリングを実行する前に適切な変換を適用する必要があります。
それぞれεij\varepsilon_{\text{ij}}平均0、分散1の単変量正規分布から得られるσ2\シグマ^{2} 、そしてすべてεij\varepsilon_{\text{ij}}は互いに独立しており（条件付き独立）、モデル係数からも独立している。 β1i\beta_{1i}そしてβ2i\beta_{2i}この仮定は妥当である。なぜなら、分析法の変動が試験対象のバッチや試験実施時点に依存する可能性は極めて低いからである。さらに、各バッチの回帰係数ベクトルは、 𝛃i\mathbf{\beta}_{i}は、それらと他のすべてのバッチのランダムエラーのベクトルから独立しており、 𝛆k\mathbf{\varepsilon}_{k} 。
RCR モデリングは 3 つのステップで構成されます。
ステップa) 各バッチのデータに単純な線形回帰モデルを当てはめる
バッチに関連付けられたデータiiは次のように表される(tij,yij)\left( t_{\text{ij}},y_{\text{ij}} \right) 、 どこyijy_{\text{ij}}を表すi番目i^{\text{th}}バッチの測定された応答の時点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}回帰の「平均二乗誤差」 MSEi=σiy^2\text{MSE}_{i} = \sigma_{i\hat{y}}^{2} 、そして自由度DFi=(ni−2)\text{df}_{i} = \left(n_{i} - 2 \right)関連するMSEi\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) 誤差分散共分散行列の推定
ランダム係数の誤差分散共分散行列の推定は、プールされた誤差共分散行列の推定から始まる。 MSE\text{MSE}すべてのバッチ、すなわちMSEプール\text{MSE}_{\text{プール}} 、 次のように。
MSEプール=σ^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{\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−β^ULS)(𝐛i−β^ULS)′\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)^{'}
どこβ^ULS=1n∑i=1n𝐛i{\hat{\beta}}_{\text{ULS}} = \frac{1}{n}\sum_{i = 1}^{n}\mathbf{b}_{i} 、つまり係数の平均ベクトル、つまり切片と傾きです（ULSは「Unweighted Least-Squares（重み付けなし最小二乗法）」の略です）。さて、分散共分散行列は𝚺^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{\ハット{\シグマ}}^{2}分析方法の分散を推定する。切片分散成分はσβ12\sigma_{\beta_{1}}^{2}は、製品出荷時に観測されるバッチ内変動のうち、製造工程のみに起因するもの、すなわち分析方法に起因する変動を除いたものを表す。傾きの変動成分はσβ22\sigma_{\beta_{2}}^{2}バッチ間の分解速度の変動を表します。
分散共分散行列は𝚺^BB{\hat{\mathbf{\Sigma}}}_{\text{ββ}} 2つの行列の差を計算して推定する場合、結果は正定値にならない可能性があります。つまり、切片および／または傾きの分散成分が負になる可能性があります。この問題への対策の一つとして、負の値を0に置き換えることが挙げられます。これは、負の分散推定値に対する簡単な応急処置であり、ランダム効果を固定効果に変換することと同等です。両方の分散が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^{\text{th}}バッチ。この式に基づいて、 2×i2 \times i切片と傾きの係数の推定値の次元行列、つまり加重最小二乗推定値は、次の式で得られます。
𝛃^WLS=(∑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^,RCR{\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つまり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つの過去のバッチと1つの現在のバッチについて様々なデータセットが生成されました。以下のコードスニペットは、RCRモデルフィッティングのステップa )、すなわち切片(a, 𝛃1\mathbf{\beta}_{1} ）、傾き（b、 𝛃2\mathbf{\beta}_{2} ）、（ルート）平均二乗誤差（（RMSE ） MSE、 MSEi\text{MSE}_{i} ）、自由度（df、 DFi\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、 MSEプール\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.
ステップc )では、RCRモデルの係数ab.hat( 𝛃^WLS{\hat{\mathbf{\beta}}}_{\text{WLS}} ) は次のコード スニペットを使用して推定されます。
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 モデルは通常の最小二乗モデルよりも優れています。5 つのケースを区別できます。ケース I: 切片の変化なし、傾きの変化なし。ケース II: 切片の変化あり、傾きの変化なし。ケース III: 切片の変化なし、傾きの変化。ケース IV: 切片の変化あり、傾きの変化。ケース V: 非線形。5 つのケースすべてについて、図 3 に結果がまとめられています。図 3 のグラフは、切片や傾きが変化するケースで RCR モデルがデータをより適切に表現し、切片および傾きの変動がない場合には基本手法と非常に似ていることを示しています。
Figure_2図_2図_2
図2 RCC / RCR方式を実装したスクリプトの出力
ウィンドウは、「過去のバッチ評価」、「現在のバッチ評価（個々のモデルパラメータ）」、そして過去のバッチデータと現在のバッチデータを要約したグラフ表示の3つの部分に分かれています。過去のバッチデータは青い点（）で、現在のバッチデータは赤い十字（）で表されます。過去のバッチデータに適合したTIは青い点線で区切られた帯で表され、現在のバッチデータ（ただし最新の観測値は含まない）に適合したPIは赤い破線で区切られた帯で表されます。
図3 さまざまな状況におけるRCRモデルの区間と基本線形回帰モデルの区間の比較
RCRモデルの適合度は、青色の実線回帰線と、対応する許容区間（95%信頼度および95%カバレッジ）の青色の点線で表されます。基本モデルの適合度は、紫色の破線回帰線と、対応する予測区間（95%信頼度）、すなわち個々の予測値の信頼限界の紫色の網掛けで表されます。パネルVaの灰色の破線は、応答変数の二乗に基づいて計算された非線形適合度を表しています。6つのパネルは、以下の状況におけるRCRモデルの区間と基本線形回帰モデルの区間を比較しています。パネルI：「切片の変動なし、傾きの変動なし」の場合。パネルII：「切片の変動あり、傾きの変動なし」の場合。パネルIII：「切片の変動なし、傾きの変動あり」の場合。パネルIV：「切片の変動あり、傾きの変動あり」の場合。パネルVa：「非線形」の場合。パネルVb：応答変数の変換（二乗）後の「非線形」の場合。
結論
本稿で紹介するJMPスクリプトは、医薬品の現在のバッチの安定性データをOOEおよびOOT結果に関して評価するための適切かつ堅牢なツールです。前提条件として、少なくとも3バッチの過去の安定性データ（各バッチについて少なくとも4つの時点）が利用可能であることが必要です。
参考文献
Carter, RL, Yang, MCK「ランダム係数回帰モデルにおける大規模サンプル推論」 統計におけるコミュニケーション – 理論と方法 (1986) 15 (8): 2507-2525. リンク
Conover, WJ (1980)、「実用ノンパラメトリック統計」、ニューヨーク:John Wiley and Sons、Inc.
ICHガイダンスQ1A(R2)、新医薬品原薬および製剤の安定性試験、ステップ4版、2003年2月6日。 リンク
マンデル、BJ「回帰管理図」 J クアルテクノロジー (1969) 1 (1): 1-9. リンク
Scholz, F., Vangel, M. バッチ効果下における許容範囲とCpk信頼限界。掲載誌：『信頼性、品質、安全性のための確率モデルにおける高度化』。編集者：Kahle, W., von Collani, E., Franz, J., Jensen, U. Birkhäuser, Boston, 1998, 361-379。 リンク
シューハート, WA (1939). 品質管理の観点からの統計的手法. (デミング, WE). ワシントンD.C. 農学部大学院.
標準操作手順 - 実験データ管理ガイダンス - 期待外（OOE）およびトレンド外（OOT）の結果。欧州コンプライアンスアカデミー（ECA）分析品質管理ワーキンググループ、クリストファー・バージェス、草案V 0.3、2015年。
英国医薬品・医療製品規制庁（MHRA）、規格外調査、プレゼンテーションノートおよび参加者用資料、2013年8月28日。 リンク
WHO技術報告書シリーズ、第953号、付録2「有効医薬品成分および完成医薬品の安定性試験、2009年」。 リンク
