code{white-space: pre;}
Pius Dahinden 氏、PhD、分析科学マネージャー、Tillotts Pharma
抽象的な
新しい EU GMP ガイドライン (EudraLex Volume 4) 第 6 章 (品質管理) が 2014 年に発効して以来、医薬品の安定性データの傾向分析の実行が追加の分析品質管理要件となっています。欧州コンプライアンス アカデミー (ECA) の分析品質管理ワーキング グループ (AQCWG) は、期待外れ (OOE) および傾向外れ (OOT) の結果を評価するための標準作業手順を策定しました。安定性データの評価には、現在のバッチと過去のバッチのデータを比較する回帰管理図 (RCC) 手順を使用することをお勧めします。ただし、過去のバッチ データ自体から推定される傾きや切片には不確実性が伴います。 ECA AQCWG は、これを考慮した単純化されたランダム係数回帰 (RCR) 方法を提案しています。推奨される RCC / RCR プロシージャを実装する、新しく開発された JMP スクリプトが紹介されています。このスクリプトにより、安定性データ内の OOE / OOT ポイントの識別が簡単になり、必要なすべての情報が適切なグラフィック表示とともに出力に提供されます。
導入
メーカーは統計的工程管理 (SPC) を適用して、生産品目の品質を確保し、生産プロセスの制御不能による経済的損失を防ぎます。プロセス制御とは、実行中のプロセスが定期的に監視され、監視されたパラメータのいずれかが異常な傾向を示す場合には、プロセスを制御下に戻すための措置が講じられることを意味します。 SPC は、メーカーがプロセスの変更を早期に検出し、トラブルの原因となっている問題を調査し、早期に是正措置を講じることができるツールです。プロセス監視のための強力なツールは、シューハートによって普及された管理図法です [6]。図 1 は、典型的な管理図の概略図を示しています。
傾向が存在しないはずの生産プロセスとは対照的に、固有の傾向を伴うデータ管理の問題が数多くあります。この種のデータ管理問題の分析には、従来の管理図法と回帰分析が組み合わされており、このアプローチはマンデルによって「回帰管理図」(RCC) と呼ばれています [4]。管理限界があるため、許容範囲 (TI) を使用するのが合理的です。推奨される保管条件 (例: 25 °C / 60% RH) で保管された医薬品の安定性データを監視するために、RCC を使用して、たとえば、時間の経過に伴う医薬品有効成分の分解や不純物の蓄積を監視し、安全性をサポートします。製品の有効期限、つまり、ラベルに表示された保存期間中の医薬品の安定性を確認するため。
Figure_1図1
図 1 典型的な管理図の構成
管理図には通常、値が表示されます。バツバツまたは、サンプル数または時間に対するサンプルの品質特性のいくつかの値 (平均など) の統計。一般に、グラフには平均値を表す中心線が含まれています。バツ￣ \overset{\overline{}}{X}制御内プロセスの。基本的な前提は、データが正規分布しているということです。追加の 2 本の特徴的な水平線は、中心線から ± 3 標準偏差 (SD) の距離にある下限管理限界と上限管理限界 (それぞれ LCL と UCL) であり、プロセスが制御されている限り、99.7% の管理限界が維持されます。データ ポイントは、共通原因変動と呼ばれるランダムな変動により、これらの制限内に収まります。中心線から ± 2 SD の距離にある線は警告下限値および上限警告限界 (それぞれ LWL および UWL) と呼ばれますが、中心線から ± 1 SD の距離にある線には特定の名前がありません。これらの線を使用して、管理図をゾーン C から A に分割できます。制御されていないプロセスを特徴付ける特殊原因変動を検出するためのルールは、これらのゾーンでの観察頻度を決定基準として使用します。
手順と結果
医薬品の安定性データに適用される RCC は、次の 2 つの質問に答えることを目的としています。 1) 現在のバッチは、過去のバッチに関して予想外 (OOE) か傾向外 (OOT) ですか (つまり、履歴データから導出された TI に関して) 、および 2) 現在のバッチに関して (つまり、その予測間隔 (PI) に関して) 現在のバッチ OOE の最も最近測定されたデータ ポイントですか?
質問 1 の傾向を評価するには、現在のバッチの少なくとも 3 つのデータ ポイントが必要です。質問 2 の評価には、最後に測定されたデータ ポイントを使用する必要がないため、少なくとも 4 つのデータ ポイントが必要です。 PIの構築。
個々の履歴バッチは切片と傾きに関して異なる可能性が非常に高いため、通常の最小二乗モデル (共通の傾き - 共通の切片) は意味がない可能性があります。このため、欧州コンプライアンス アカデミー (ECA) 分析品質管理ワーキング グループ (AQCWG) は、個々のバッチの切片と傾きの変動を考慮した「ランダム係数回帰」 (RCR) 手法を利用することを推奨しています [7]。単純化された RCR モデル (修正された Carter および Yang [1] モデル) には、2 つの未知のパラメーター、つまり切片が含まれています。 β 1 \ベータ_{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}}
どこ 𝐲私\mathbf{y}_{i}ですn私× 1 n_{i} \times 1を表すベクトル私番目i^{\text{th}}バッチのn私n_{i}繰り返し測定（ここで、 j = 1 、 ⋯ 、 n私j = 1,\ \cdots,\ n_{i} 、 とn私n_{i}バッチの繰り返し測定の総数を表します。私私) 応答変数 (例: アッセイ [%]、不純物の合計 [%]、または溶解 [%])、 𝛃 1私\mathbf{\beta}_{1i}そして 𝛃 2私\mathbf{\beta}_{2i}はそれぞれバッチ固有の切片と傾きのベクトルです。 𝐭私\mathbf{t}_{i}それは私番目i^{\text{th}}バッチの時間ベクトル (例: [月])、および 𝛆私\mathbf{\varepsilon}_{i}ですn私× 1 n_{i} \times 1のベクトルn私n_{i}ランダムなエラー。この 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 モデリングは 3 つのステップで構成されます。
ステップ a) 単純な線形回帰モデルを各バッチのデータに当てはめる
バッチに関連付けられたデータ私私で表されます ( t ij 、 y ij ) \left( t_{\text{ij}},y_{\text{ij}} \right) 、 どこy ij y_{\text{ij}}を表します私番目i^{\text{th}}ある時点でのバッチの測定された応答t j t_{j} 、 どこj = 1 、 ⋯ 、 n私j = 1,\ \cdots,\ n_{i} 。個々の時点t j t_{j}バッチ間で異なる場合があり、時点の合計数も異なる場合があります。 n私n_{i} 。デザインマトリックス 𝐗私\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} 、回帰の「平均二乗誤差」 MSE私= σ私y ^ 2 \text{MSE}_{i} = \sigma_{i\hat{y}}^{2} 、および自由度 DF私= ( n私− 2 ) \text{df}_{i} = \left( n_{i} - 2 \right)と関連した MSE私\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{プール}} 、 次のように。
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=SSEi=𝛆^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−β^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 = 1 n ∑私= 1 n 𝐛私{\hat{\beta}}_{\text{ULS}} = \frac{1}{n}\sum_{i = 1}^{n}\mathbf{b}_{i} 、つまり係数の平均ベクトル、つまり切片と傾きです (ULS は「重み付けされていない最小二乗法」の略です)。さて、分散共分散行列 𝚺 ^ ββ {\hat{\mathbf{\Sigma}}}_{\text{ββ}}は以下のように推定できる。
𝚺^ββ=𝐒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 \transrated{\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 {\hat{\sigma}}^{2}分析手法の分散を推定します。切片分散成分σ β 1 2 \sigma_{\beta_{1}}^{2}は、製造プロセスのみに起因する、製品リリース時に観察された応答のバッチ内分散を表します。つまり、分析方法に起因する分散は除外されます。傾き分散成分σ β 2 2 \sigma_{\beta_{2}}^{2}は、分解率のバッチ間の分散を表します。
分散共分散行列なので、 𝚺 ^ ββ {\hat{\mathbf{\Sigma}}}_{\text{ββ}}は 2 つの行列の差を計算して推定されますが、結果は正定値ではない可能性があります。つまり、切片および/または傾き分散成分が負になる可能性があります。これに対する考えられる解決策の 1 つは、負の値を 0 に置き換えることです。これは、負の分散推定に対する簡単な応急処置のアプローチであり、変量効果を固定効果に変換するのと同じです。両方の分散がゼロの場合、RCR モデルは通常の最小二乗 (OLS) モデルに対応します。
ステップ c) 重み付き最小二乗 (WLS) 推定器による係数平均ベクトルの推定
係数平均ベクトルの推定は、次の推定値の計算から始まります。 𝐖私\mathbf{W}_{i} 、
𝐖i=(𝚺^ββ+σ^2•𝐌i)−1\mathbf{W}_{i} = \left( {\hat{\mathbf{\sigma}}}_{\text{ββ}} + {\hat{\sigma}}^{2} \bullet \mathbf {M}_{i}\右)^{-1}
これは、回帰パラメータ推定値の分散の逆数に対応します。私番目i^{\text{th}}バッチ。この方程式に基づいて、 2 ×私2 \倍私切片と傾き係数の推定値の次元行列、つまり重み付き最小二乗推定値は、次の方程式によって取得されます。
𝛃^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 \ベータ_{1}そしてβ 2 \beta_{2}全体的な傾向、つまり過去のバッチの共通の傾向を説明します。 RCR モデルの TI は、推定値を利用して計算できます。 y ^ j {\hat{y}}_{j}回帰直線とその標準誤差についてσ y ^ 、 R C R \sigma_{\hat{y},RCR} Scholz と Vangel によって提案された方法 [5] に従って、次のようになります。彼らの方法では、有効サンプルサイズの概念が利用されています。
y^j±k•σy^,RCR{\hat{y}}_{j} \pm k \bullet \sigma_{\hat{y},\text{RCR}}
エラー用語σ y ^ 、 R C R \sigma_{\hat{y},RCR}によって計算されます
σy^,RCR=D(𝐭(𝚺^ββ+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}} |
どこσ b 2 \sigma_{b}^{2}はバッチ間の分散、 σ e 2 \sigma_{e}^{2}はバッチ内分散 (またはメソッド分散)、 n私n_{i}バッチごとの測定数です私私、つまり私= 1 、 ⋯ 、 n b i = 1,\ \cdots,\ n_{b} 、 n b n_{b}は過去のバッチの数であり、 n nはサンプルの総数です。つまり、 n = ∑私= 1 n b n私n = \sum_{i = 1}^{n_{b}}n_{i} 。
TI がカバーしますp ( 100 % ) p\left( 100\% \right) ～の確率で母集団のγ ( 100 % ) \gamma\left( 100\% \right) 。有効サンプル サイズの概念では、バッチ内のデータはバッチ間データよりも相関性が高く、独立した情報の量が減少することが考慮されています。このため、次のことが常に当てはまります。 n * ≤ n n^{*} \leq n 。理想的なケースn * = n n^{*} = 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、 MSE私\text{MSE}_{i} )、自由度 (df、 DF私\text{df}_{i} )、および観測値の数 (n, n私n_{i} )の私私履歴バッチが収集され、長さのベクトルに保存されますn b 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
);
ステップb ) は、プールされた MSE (MSE.pool, MSEプール\text{MSE}_{\text{プール}} ）。続いて、平均モデルの分散共分散行列 (M.bar, 𝐌 ￣ \overset{\overline{}}{\mathbf{M}} ) が推定され、サンプルの分散共分散行列 (S, 𝐒 bb \mathbf{S}_{\text{bb}} ) 分散共分散行列を取得します (Sigma.hat, 𝚺 ^ ββ {\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 RCC / RCR メソッドを実装するスクリプトの出力
このウィンドウは、「履歴バッチ評価」、「現在のバッチ評価 (個々のモデル パラメーター)」、および履歴および現在のバッチ データを要約するグラフィカル表示の部分に分かれています。過去のバッチのデータは青い点 () で表され、現在のバッチのデータは赤い十字 () で表されます。過去のバッチのデータに当てはめた 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 Qual Technol (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 年)。品質管理の観点からの統計的手法。 （デミング、我々）。ワシントン DC、農学部大学院。
標準操作手順 – 検査データ管理ガイダンス – 期待外れ (OOE) および傾向外れ (OOT) の結果。欧州コンプライアンス アカデミー (ECA) 分析品質管理ワーキング グループ、Christopher Burgess、ドラフト V 0.3、2015 年。
英国医薬品・ヘルスケア製品規制庁 (MHRA)、規格外調査、プレゼンテーションノートおよび参加者パック、2013 年 8 月 28 日。リンク
WHO 技術報告書シリーズ、No. 953、付録 2、医薬品有効成分および最終医薬品の安定性試験、2009 年。 リンク
