code{white-space: pre;}
ピウス・ダヒンデン博士、ティロッツ・ファーマ、分析科学マネージャー
抽象的な
新しい EU GMP ガイドライン (EudraLex Volume 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 °C / 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モデルは次のように記述できる。
そして私j=𝛃1i+𝛃2i•𝐭私j+𝛆私j\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}}バッチのniニ_{i}繰り返し測定（ここでj=1,⋯,nij = 1,\ \cdots,\ n_{i} 、 とniニ_{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ベクトルのniニ_{i}ランダム エラー。この RCR モデルの前提条件は、応答と説明変数の間に線形関係があることです。関係が線形でない場合は、RCR モデリングを実行する前に適切な変換を適用する必要があります。
それぞれε私j\varepsilon_{\text{ij}}平均0、分散1の単変量正規分布から得られるσ2\シグマ^{2} 、そしてすべてε私j\varepsilon_{\text{ij}}は互いに独立しており（条件付き独立）、モデル係数からも独立している。 β1i\beta_{1i}そしてβ2i\beta_{2i}この仮定は妥当である。なぜなら、分析法の変動が試験されるバッチや試験が実行される時点に依存する可能性は非常に低いからである。さらに、各バッチの回帰係数ベクトルは、 𝛃i\mathbf{\beta}_{i}は、それらと他のすべてのバッチのランダムエラーのベクトルとは独立しており、 𝛆k\mathbf{\varepsilon}_{k} 。
RCR モデリングは 3 つのステップで構成されます。
ステップa) 各バッチのデータに単純な線形回帰モデルを当てはめる
バッチに関連付けられたデータiiは次のように表される(t私j,y私j)\left(t_{\text{ij}},y_{\text{ij}} \right) より 、 どこy私jy_{\text{ij}}を表すi番目i^{\text{th}}バッチの測定された応答tjt_{j} 、 どこj=1,⋯,nij = 1,\ \cdots,\ n_{i}個々の時点tjt_{j}バッチ間で異なる場合があり、時間点の合計数も異なる場合がある。 niニ_{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−β^UL)(𝐛i−β^UL)′\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)^{'}
どこβ^UL=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{\hat{\sigma}}^{2}分析方法の分散を推定する。切片分散成分σβ12\sigma_{\beta_{1}}^{2}は、製品出荷時に観察されるバッチ内変動を表し、製造プロセスのみに起因するものであり、分析方法に起因する変動は除外されている。傾き変動成分σβ22\sigma_{\beta_{2}}^{2}バッチ間の分解速度の変動を表します。
分散共分散行列𝚺^BB{\hat{\mathbf{\Sigma}}}_{\text{ββ}} 2 つの行列の差を計算して推定する場合、結果は正定値ではない可能性があります。つまり、切片および/または傾きの分散成分が負になる可能性があります。これに対する 1 つの解決策は、負の値を 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切片と傾き係数の推定値の次元行列、すなわち加重最小二乗推定値は、次の式によって得られます。
𝛃^ウェブ=(∑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^,RCRC{\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*−11 は 1 です。自由度と非心度パラメータzp•n*z_{p} \bullet \sqrt{n^{*}} 、 どこzpz_{p}はpp標準正規分布の分位数。パラメータ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}バッチ内分散（またはメソッド分散）である。 niニ_{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、 niニ_{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( 𝛃^ウェブ{\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方式を実装したスクリプトの出力
ウィンドウは、「履歴バッチ評価」、「現在のバッチ評価 (個別のモデル パラメータ)」、および履歴と現在のバッチ データを要約したグラフィック表示の部分に分かれています。履歴バッチのデータは青い点 () で表され、現在のバッチのデータは赤い十字 () で表されます。履歴バッチのデータに適合した TI は青い点線で区切られたバンドで表され、現在のバッチのデータに適合した PI (最新の観測値は含まない) は赤い破線で区切られたバンドで表されます。
図3 さまざまな状況におけるRCRモデル区間と基本線形回帰モデル区間の比較
RCR モデルの適合は、青色の実線回帰線と、対応する許容区間の青色の点線 (95% の信頼度と 95% の適用範囲) で表されます。基本モデルの適合は、紫色の破線回帰線と、対応する予測区間 (95% の信頼度) の紫色の陰影 (つまり、個々の予測値の信頼限界) で表されます。パネル Va の灰色の破線は、応答の二乗に基づいて計算された非線形適合を表します。6 つのパネルは、次の状況で RCR モデル区間と基本線形回帰モデル区間を比較します。パネル I:「切片の変化なし、傾きの変化なし」の場合。パネル II:「切片の変化あり、傾きの変化なし」の場合。パネル III:「切片の変化なし、傾きの変化」の場合。パネル IV:「切片の変化あり、傾きの変化」の場合。パネル Va:「非線形」の場合。パネル Vb: 応答変数の変換 (つまり、二乗) 後の「非線形」の場合。
結論
この原稿で紹介されている JMP スクリプトは、医薬品の現在のバッチの安定性データを OOE および OOT 結果に関して評価するための適切で堅牢なツールです。前提条件は、各バッチの少なくとも 4 つの時点を含む、少なくとも 3 つのバッチの過去の安定性データが利用可能であることです。
文献
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 Qual Technol (1969) 1 (1): 1-9。 リンク
Scholz, F.、および Vangel, M. バッチ効果における許容範囲と Cpk 信頼範囲。掲載先: Advanced in Stochastic Models for Reliability, Quality and Safety。編集者: Kahle, W.、von Collani, E.、Franz, J.、および Jensen, U. Birkhäuser、ボストン、1998、361-379。 リンク
シューハート、WA (1939)。品質管理の観点から見た統計的手法。(デミング、WE)。ワシントン DC、農学部大学院。
標準操作手順 - ラボデータ管理ガイダンス - 予想外 (OOE) および傾向外 (OOT) の結果。欧州コンプライアンス アカデミー (ECA) 分析品質管理ワーキング グループ、Christopher Burgess、ドラフト V 0.3、2015 年。
英国医薬品・医療製品規制庁 (MHRA)、規格外調査、プレゼンテーション ノートおよび代表者パック、2013 年 8 月 28 日。 リンク
WHO技術報告書シリーズ、第953号、付録2「医薬品有効成分および医薬品完成品の安定性試験、2009年」。 リンク
