☆
もう4年以上前に、”アレニウス式加速試験40℃~60℃は25℃の何ヶ月かExcelで計算”というタイトルで、加速試験結果から長期保存安定性を推定する記事を書きました。
今回それに関連して、統計解析により医薬品の有効期間を推定する方法について記事にしようかと思います。
というのも、色々なところでその解析方法の有料セミナーが開催されているようで、KenUはセミナーに参加したことはありませんがMicrosoft Excelを使って計算してみたので。
*
さて、長期保存安定性試験結果の範囲を超えて有効期間を外挿(※1)する方法が、ICH(日米EU医薬品規制調和国際会議)ガイドライン “ EVALUATION FOR STABILITY DATA Q1E ” に示されています。
厚生労働省は、これを翻訳した「安定性データの評価に関するガイドラインについて」を平成15年6月3日 医薬審発第0603004号 厚生労働省医薬局審査管理課長通知として発出しています。
しかし、具体的にどのような統計モデルに基づいて外挿したらよいのか理解が難しいことから、製薬協(日本製薬工業協会)(←HPリンク)が、統計解析の進め方について解説した統計・DM部会資料(以下「資料」という。)(※2)を作成しています。
————————————-
※1:外挿(がいそう、英: Extrapolation)とは、ある既知の数値データを基にして、そのデータの範囲の外側で予想される数値を求めること。出典Wikipedia
※2:医薬品評価委員会、統計・DM部会資料、新有効成分含有医薬品の安定性試験データの評価-リテスト期間及び有効期間設定のための試験デザインと統計的推定法-
————————————-
この ” 資料 ”(←参照リンク・PDFファイル)は、とっても親切で分かりやすいと思います。
で、そこには「読者には、本章の数値例を表計算ソフトウェア(例えば MS-EXCEL)や電卓を用いて、自分で確認することを薦める。」と書かれています。
なので、資料の数値をそのまま使って、試しにエクセルの関数や数式を入力してトレースしてみたんです。
*
まず1ロット目のデータ解析。
下図のエクセルシートは、“単一ロットの場合の統計的推定”と“複数ロットの場合の併合誤差分散を用いた統計的推定”を一つにまとめています。
1 ロットの各測定時点で少なくとも3回の繰り返し数の測定が行われますが、測定値の分散は全ての時点で等しいことを仮定し、平均値を解析に用います。
数字もグラフも資料を再現できました。
なお、Excel分析ツール、マクロ、VBAなどは使っていません。
各セルに入力した数式はこちら↓(図をクリックすると別タブで開きます。)
F列6行(t 分布の片側逆関数)の値は、t分布表から拾ってきましたが、エクセルの関数 TINV(確率,自由度)を使っても良いです。
ただし、TINVは両側t分布の関数なので片側t分布の場合には確率×2にし、つまり片側0.05の場合は0.1を入力し、=tinv(0.1,3)とすると2.353363434801820になります。
また、併合誤差分散を求める数式の値は、3つのロットのSSEの数値を直接入力しましたが、各ロットのシートへリンクしても良いです。
ちなみに、ロット1のシートをコピーして含量データだけ変えた他2ロットのシートの図はこちら↓
3ロット分の経時変化の推定を図にまとめると、こんな図↓になります。
併合誤差分散を使って計算すると、推定の有効期限が長くなることが分かります。
*
その他、資料には上記3ロットの ” 一括評価に関する検定と統計的推定 ” についても解説が載っています。
ですが、ここまで分かればそれも難なく計算できると思うので、このブログでの説明は省略します。
資料で少し分かり難かったところだけ捕捉すると・・・
傾きの一様性の検定に関する解説に「F値=0.96となり、p=0.418となる」と書いてあって、そのp値の求め方で少し考え込んだのですが、EXCELのFDIST関数を使うと・・・
=FDIST(F値,C-1,n^total-2c)
=FDIST(F値,ロット数-1,全観測数-2×ロット数)
=FDIST(F値,3-1,15-2×3)
=FDIST(0.9625,2,9)
=0.4180
となり、
切片の一様性の検定の「F値=7.53となり、p=0.009となる」というところも同様にFDIST関数を使うと・・・
=FDIST(F値,C-1,n^total-c-1)
=FDIST(F値,ロット数-1,全観測数-ロット数-1)
=FDIST(F値,3-1,15-3-1)
=FDIST(7.5328,2,11)
=0.0087
となりました。
資料では有効期間は29.6 箇月と推定されていて、KenUがエクセルで計算してグラフを書かせたのが下図。はい、含量が95%になるのは、資料とほぼ同じ29.7か月でした(※小数点以下を丸めずにそのまま計算しているので少しのずれがあります)。
*
ということで、意外と簡単にできました?
観測数nが変わる場合にはX値,Y値の入力範囲、自由度、t値などが変わってくるので式の見直しは必要ですが、簡単に修正できるのでは?と思います。
頑張ってください。
☆