第10回 生物検定法の故郷、2値データの用量反応に対する逆推定，平行性検定

Transcription

1 i 生物検定法の故郷 値データの 用量反応に対する逆推定, 平行性検定 第 回高橋セミナー 高橋行雄 ファイル名 :C: Documents and Settings Owner My Documents anz_seminal_9 Semi9_ 生物検定法の故郷 _ 修復 3.doc 最終保存日 :4// :9 AM 最終印刷日時 :5 年 月 3 日 /9 時 4 分

2 ii 目次 表紙裏 改訂の記録 年 月 3 日, 新規作成

3 iii 目 次. はじめに.... プロビット変換, ロジット変換 プロビット変換 ロジット変換 値反応に対する解析法 回帰曲線の 95% 信頼区間 逆推定 推定された有効量についての標準誤差の近似 フィラーの定理を用いた有効用量の信頼区間 ED5 の信頼区間の計算事例 最尤法による計算の例 吉村の例 JMP の行列計算言語による SAS の POBIT プロシジャによるプロビット法による計算 JMP の行列言語によるロジスティック回帰... 7

4 iv 目次 図表目次 表. LD 5 を求める数値例... 4 表. プロビット法による解析結果... 4 表.3 ロジスティック回帰モデルの結果... 8 表.4 プロビット法とロジット法の推定値 p の差... 8 表.5 ロジット法による逆推定... 9 表.6 ロジット法による 95% 信頼区間... 図. プロビット曲線の当てはめ... 5 図. プロビット曲線とロジット曲線の比較... 7 図.3 ロジット曲線と 95% 信頼区間の表示... 図.4 ED5 とその信頼区間の近似計算 Excel による計算... 9 図.5 ED5 とその信頼区間のフィラーの公式 Excel による計算... 図.6 ED5 とその信頼区間の近似計算... 図.7 ロジット法のフィラーの公式による 95% 信頼区間... 図.8 回反復のプロビット法のフィラーの公式による 95% 信頼区間...

5 第 回 - 高橋セミナーのご案内 年 4 月 日, 医薬安全研の午前中にセミナーを行います. 第 9 回のテーマ : 生物検定法の故郷, 値データの用量反応に対する逆推定, 平行性検定 をすべて終了することは不可能でした. 第 回目は, 第 9 回目で芳賀先生から指摘を受けた改定した ロジット曲線とプロビット曲線の比較 も含め, 値データの解析に用いられている最尤法, ニュートン ラフソン法の復習を行い, 5. 相対力値 以後をテーマとします.. プロビット変換, ロジット変換. プロビット法とロジスティック回帰によるシグモイド曲線の推定 3. 逆推定によるLD5の推定 4.LD5の95% 信頼区間 5. 相対力値 6. 自然反応 7. 非線形ロジット回帰モデル 8.complementary log-log 9. ランダム効果を含むつの用量反応データの比較このテーマは, 題 回目の課題でも有りました. この時は,JMPで取り扱うことがどこまで可能なのか, 本当にできるのか, に集中した盛りだくさんの内容でした. この内容は, biostat/ からダウンロードできます. 今回は, 実験研究者が引用できる文献を紹介しつつ,JMPでの解析の演習を行います. これにより, 実験研究者, プロウラマー, 統計家が実務でJMPあるいは他の統計パッケージの見直しができるようになります. 参加希望の方はJMPがインストールされたノートPCを準備してください. 文献は,D. Collett (99) Modeling Binary Data, Chapman and Hall/CRC の第 4 章生物検定法といくつかの適用です そして, ピンク本 5.4 節 LD5の推定です. 時間があれば,Overdispersioの問題にもふれたい. この問題は, 次のような実験データの解析に現れる. 6.. Variation between the response probabilities 実験ユニットでバッチの数が, 同じ条件で得られたとき, 応答確率がバッチとバッチの間で, それにもかかわらず異なるかもしれない. 胎児の例, 催奇形性試験として知られている. 母体の遺伝的影響により同じ実験条件でも奇形出現率は異なる. 奇形出現率の分散は, それが定数の場合よりも大きい.

6 この変動を説明する変数を実験者は知ることができない. 実験条件が十分にコントロールされていない場にも起きる. Example 6. Germination of Orobanche ( 岩波生物学辞典 : シソ花類 ; ハマウツボ ) :ageyp75 :ageyp75 :Bean :Cucumber :Bean :Cucumber y n y n y n y n インゲン豆ときゅうりの引き抜いた根の中での つのハマウツボ種種の発芽の割合を 考える.

7 . はじめに 生物検定法の故郷, 値データの用量反応に対する逆推定, 平行性検定について基礎にもどって講義と実習を行う このテーマは, 第 回目セミナーの課題であった この時は,JMP で生物検定法の代表的な課題に対してどこまで対応可能なのか, 本当にできるのか, などについての解析結果を主体にした ( からダウンロードできる ) 今回は, 実験研究者が 引用できる文献を紹介しつつ,JMP での解析の演習を行う. これにより, 実験研究者, プログラマー, 統計家が実務で JMP あるいは他の統計パッケージを利用した生物検定法の見直しができるように意図した. 今回の講義と実習は D. Collett (99) Modeling Binary Data をベースにする この本は 統計モデルを用いた豊富な薬理試験 毒性試験などの生物系の実験データの解析事例を用いつつ 値データの解析の理論的な背景が丁寧に解説されている 生物検定法は 4 章 Bioassay and some applications で取り上げられている 参加者になじみのある吉村編著 (987) 毒性 薬効データの統計解析 5.4 節 LD5 の推定法 の事例を取り上げる ロジスティック回帰については 丹後 (996) ロジスティック回帰 SAS を利用した統計解析の実際 がロジスティック回帰全般について丁寧に述べられている この本にも 3.9 節で LD5 ED5 の推定が取り上げられている 生物検定法 (Biological Assay) とは生物を用いて未知の化合物の生物活性を既知の化合物の生物活性に対して相対的に比較するために体系化された応用統計学の一つの分野である代表的な生物検定法 : 5 パーセント致死量の推定 生物検定法は回帰分析の応用直線の当てはまりの欠如 (Lack of Fit) 非平行性 (Lack of Parallelism) 逆推定とその信頼区間用量反応関係を論ずるために欠かせない他の応用統計の分野では軽視

8 SAS とJMPによる生物検定法 SAS: proc PROBIT 逆推定を取り扱えるのは 群の場合プロビット変換よりむしろロジット変換ロジスティック回帰分析として定式化 SAS: proc LOGISTIC proc CATMOD proc GENMOD JMP: Fit Y by X Fit Model 逆推定 Y となる X は 生物検定法で常用される逆推定 95% 信頼区間 SAS では標準的には求められない JMP では Inverse Prediction の問題として対応生物検定法のための統計パッケージとして JMP が適している 今回は次の課題を 基礎概念も含めて再度 取り上げ 非線形混合モデルの応用の可能性についても触れてみたい. プロビット変換, ロジット変換. プロビット法とロジスティック回帰によるシグモイド曲線の推定 3. 逆推定による LD5 の推定 4.LD5 の 95% 信頼区間 5. 相対力値 6. 自然反応 7. 非線形ロジット回帰モデル 8. 重対数変換 9. ランダム効果を含む つの用量反応データの比較

9 プロビット変換, ロジット変換 3. プロビット変換, ロジット変換.. プロビット変換毒性 薬効データの統計解析の 5.4 節 LD5 の推定法生物検定法では プロビット法による LD5 の推定とその 95% 信頼区間の計算法が示されている を毒性の目安とするときには 暗黙のうちに 用量と反応の関係が単調増加 LD 5 であることを前提にしている その単調増加曲線としては ほとんどの場合プロビット曲線を想定している プロビット曲線 probit curveとは 図. に示すように反応率が 用量の対数に対して 正規分布の累積確率 すなわち正規分布関数の関係をもつ曲線のことである 見方を変えれば 用量を横軸にとり 反応率を縦軸に取って データを対数正規確率紙にプロットしたときに 関係が正の勾配を持った直線になることである 式で書けば用量 d と反応率 p ( 原文では y ) の関係が p = log d ( x μ) exp{ } σ dx πσ (.) となるものである このとき μ = log LD5 である probit は probability unit の省略である 以下に プロビット曲線を前提にして 実験データから LD 5 を推定する手法を紹介 するが 計算では次の式で定義するプロビット関数を使う y 5 x exp{ } probit( p) p = = dx となるyの値 (.) π y は基準正規分布の累積確率が p となるところの横軸 ( 正規偏差 ) に 5 を加えたものである 数値例として表. が示され 本的な LD 5 の求め方として重み付き最小 乗最の繰 返し計算による最尤法が示されている 計算手順は煩雑なので JMP の行列計算言語に よるプログラミング例と結果を付録に示した ここでは 得られた結果の解釈を通じて プロビット法の考え方を示す 式 (.) のプロビットは 標準正規分布の逆関数 Φ ( p) に 5 を加える JMP では Normal Quartile 関数を用いることにより直接計算できる 死亡 率が % あるいは % は マイナス無限大とプラスの無限大となる

10 4 プロビット変換, ロジット変換 表. を求める数値例 LD 5 群 投与量 mg/kg 群の大きさ 死亡数 死亡率 プロビット i ( 公比.35) p i probit( p i ) 投与量の常用対数をとり 変換されたプロビットに 重み付き最小 乗法の繰り返し により尤度を最大化するような次の回帰直線をあてはめる 3. 節に示した JMP の行列 計算プログラムにより 次の回帰係数が得られる η = probit( ˆ ) = log ( ) i pi dosei 表. に この式によって計算されたプロビットの推定値を示す 死亡率の推定値は 標準正規分布 Φ ( η i 5) よって計算される 表. プロビット法による解析結果 群 log ( dose i ) probit( p ˆ i ) p ˆi 死亡率 p が.5 となる 投与量は式 (.) の μ とσ は η = ˆ β + ˆ β x i i が x ( ( ˆ β ˆ 5)/ β) ηi = / ˆ β と書き換えられるので ˆ μ = log LD = (( ˆ β 5) / ˆ β ) = (.44 5) / 6.63) =.78 5 となり

11 プロビット変換, ロジット変換 5.78 LD = = 9. mg/kg 5 が得られる σ は ˆ σ = / ˆ β = / 6.63 =.54 と計算される これらの μ =.78 σ =.54 となる正規分布 および.4 節で述べる方法でもとめた回帰直線の 95% 信頼区間を逆プロビット変換して 当てはめたのが図. のシグモイド曲線である Y log(dose) 図. プロビット曲線の当てはめ LD 5 の 95% 信頼区間は 図. のY 軸の.5 を通る水平線上の 95% 信頼区間を横切る ときの X 軸の目盛りを読むことにより得られる 解析的には フィラーの式により 常 用対数で (.6,.347) となり 投与用量に変換して (6, 3) が得られる フィラー の式は.7 節で詳しく述べるので ここでは結果のみを示した ここで示した計算結果は 重みつき最小 乗法を 回繰り返した場合であり 数値計 算的には収束していない プロビット法で十分に収束した結果を求めるためにSASのproc PROBITを用いた結果 (3. 節に計算プログラムと結果 ) を示す これは JMPでは プロビット法の計算はサポートされていないためである なお SASでは proc LOGISTIC proc GENMODでもプロビット法での計算は行えるが LD5 の 95% 信頼区間の計算がサポートされていない

12 6 プロビット変換, ロジット変換 SAS の proc PROBIT の反復計算は 初期の回帰係数 β = β = からスタートして 6 回のニュートン ラフソン法による反復の結果 β = 5.73 β = 6.65 が得られ た 切片が 回の重み付き最小 乗法で得られた -.44 に比べてかなり小さくなって いるが これは プロビット変換時に定数 5 を加えたためであり -5 を加えると となりほほ同じである LD5 の 95% 信頼区間は フィラーの式により 常用対数で (.4,.348) と計算され 投与用量に変換して (6, 3) と JMP の重み付最小 乗法の 回の反復と計算誤差範囲内で一致している.. ロジット変換シグモイド曲線となる 値反応の出現率に 正規分布をあてはめることにより さまざまな生物検定法における応用がされてきた シグモイド曲線を得るために正規分布の数値計算は煩雑であることから 数値計算が簡単なロジスティック分布をシグモイド曲線に用いる方法が多用されるようになってきた f( x) = x μ exp τ x μ τ + exp τ, < x < (.3) ここで < μ < τ < であり 平均と分散は それぞれ μ とπτ /3である 確率 密度関数 f ( x) は 正規分布に比べ簡潔とはいえないが 分布関数は p i x μ log d μ log d i i exp exp τ τ dx (.4) x μ log di μ + exp = = τ + exp τ τ と簡潔な指数関数となる ここで β = μ/ τ β = / τ とおけば ロジスティック分布関数は ( β + β logdi ) ( β β d ) exp pi = (.5) + exp + log となる 簡単な式の変形により ロジット i p i logit( pi) = ln = β + β log pi d i

13 プロビット変換, ロジット変換 7 が得られる 式 (.5) が正規分布を用いた場合とほぼ同様のシグモイド曲線を与える このシグモイド曲線を用いて 5% 致死量 LD5 を推定するのがロジスティック回帰モデ ルである ロジスティック法で推定された回帰係数は プロビット法のそれとは異なる プロ ビット法では / β が標準偏差となるが ロジスティック法での標準偏差は πτ /3であり プロビット法で推定された推定された標準偏差を ˆ σ とすれば となる 3 ˆ σ 3.54 ˆ τ = = =.83 π 3.4 JMPによるLD5 の推定は プロビット法ではなく ロジット法 による推定となっている プロビット法によって得られた LD 5 =.78 と ˆ σ =.54 から推定されるプロビット曲線 LD 5 =.78 と ˆ τ =.83から推定されるロジスティック曲線を比較した結果を図. に示す 違いは ロジスティック回帰が裾広がりとなるが 図にしてみるとごくわずかである..8.6 Y x 図. プロビット曲線とロジット曲線の比較 プロビット ロジット

14 8 プロビット変換, ロジット変換 JMPでロジスティック回帰を行った結果を表.3 に示す 回帰係数 β =.59 と勾 配が急になり 裾広がりを解消している 表.3 ロジスティック回帰モデルの結果 パラメータ推定値項推定値切片 log(dose).5999 標準誤差 カイ 乗 p 値 (Prob>ChiSq) <. <. 5% 致死量は 求められた回帰係数から ˆ μ = log LD = ˆ β / ˆ β = ( 6.5) /.59 =.75 5 となり.75 LD = = 88. mg/kg 5 と推定できる プロビット法の 9.mg/kg と比べて約 % の違いである 分散は πτ 3.4 (/.5) ˆ σ = = = とプロビット法の.5 とほぼ同じである プロビット法とロジット法による推定値 p の差を表.4 に示す その結果でもロジット法の結果は 裾広がりとなっているが その差は.5% 程度の違いとわずかである 表.4 プロビット法とロジット法の推定値 p の差 log ( dose i ) Probit logit logit-probit JMP では 逆推定の機能により 任意の死亡率について計算できるので 9% 5%

15 プロビット変換, ロジット変換 9 および % のそれぞれについて推定した結果を表.5 に示す LD 5 は フィラーの式により 常用対数で (.96,.35) と推定される の 95% 信頼区間 表.5 ロジット法による逆推定 逆推定確率.9.5. 予測値 log(dose) 下限 上限 Alpha.95 JMP のロジスティック回帰モデルの計算手順は ニュートン ラフソン法による反復 計算による尤度が最大となるような解が求められている どのような計算手順なのかを 示すために 3.3 節 JMPの行列計算プログラムと実行結を示す 反復計算は 初期の回帰係数 β = β = からスタートして 6 回の反復の結果 β = 6. β =.53が得 られ 結果は一致する.3. 値反応に対する解析法 これまでに多くの 法 を使ってきたので これらの関連を整理してみよう 値データの変換法の反復プロビット変換 (5 を加える場合も有る ) ロジット変換シグモイド曲線をあてはめて 5% 推定値を求める方法の総称プロビット変換による場合 : プロビット法ロジット変換による場合 : ロジット法 または ロジスティック回帰モデル最尤解を求める反復計算の方法 どちらの方法でも解は一致する重みつき最小 乗法ニュートン ラフソン法 方法の組合せプロビット変換 重み付き最小 乗法の反復吉村らの本 SAS/LOGISTIC プロビット変換 ニュートン ラフソン法

16 プロビット変換, ロジット変換 SAS/PROBIT( 逆推定 ) SAS/LOGISTIC SAS/GENMOD ロジット変換 重み付き最小 乗法の反復 SAS/LOGISTIC ロジット変換 ニュートン ラフソン法 JMP( 逆推定 ) SAS/LOGISTIC SAS/GENMOD.4. 回帰曲線の 95% 信頼区間 プロビット曲線またはロジスティック曲線の 95% 信頼区間を求めよう η = x β ˆ の 推定値の分散は 分散共分散行列を Σ としたときに σ = x Σ x x i i となる 単回帰分析の場合は η ˆ ˆ i = β + log ( dosei) β なので ˆ β の分散を v ˆβ の分散を v ˆ β と ˆβ の共分散を で計算される σ = x Σ x = v + log ( dose ) v + log ( dose ) v x i i i i v とすると JMPのロジスティック回帰で求められた回帰係数は ˆ β = 6. ˆ β =.5 JMP の 行列計算言語で計算した分散共分散行列は Σ= である η の場合につい計算してみる η = ˆ β + log ( dose ) ˆ β = 6.+ log ().5 = 3.59 i e p = + e i 3.5 ˆ 3.59 σ = x Σ x = v + log ( dose ) v + log ( dose ) v x i i i i = log () ( 7.63) + log () 7.73=.739 η ±.96σ x = 3.59 ± = ( 4.8,.43) e e pˆ ( L95) = =.8, pˆ 4.8 ( U95) = = e + e i i 他の用量についても同様に計算した結果を表.6 に示す

17 プロビット変換, ロジット変換 表.6 ロジット法による 95% 信頼区間 i dose n Y p log (dose) eta p^ Var(eta) s.e.(eta) L95% U95% 投与量を連続的に変えて推定値とその 95% 信頼区間計算した結果を図.3 に示す 図の右は拡大した図であり 反応が.5 となるような 95% の投与量が逆推定されている Y.6.4 Y log(dose) log(dose) 図.3 ロジット曲線と 95% 信頼区間の表示.5. 逆推定応答確率が 単一の説明変数 x の関数としてモデル化される場合 応答確率に対応する x の値を推定することにしばしば興味がある 例えば 生物検定法では 化学薬品に曝露された個体の 5% に反応すると期待される濃度に しばしば興味がもたれる この用量はメディアン用量と名付けられ 一般に ED5 値と呼ばれる 試験の反応が死である場合 これはメディアン致死量 あるいは LD5 値と言い換えられる 他の反応率を期待する濃度は 同様の方法出で表される 例えば 各個体の 9% に反応させる用量は ED9 値などである これらの量での要約は 分析しようとしている化学薬品の力価 および その後の異なる物質の間の比較の基礎となる 生物検定法により応答確率 p および用量 d の関連を記述するために 次に示す

18 プロビット変換, ロジット変換 線形ロジスティックス モデルを使用するとしよう p logit( p) = log = β + β d p p =.5 が ED5 値である用量は logit (.5) = log () = であるので ED5 値は 式 β + β = ED5 となり その結果 ED5 = β / β となる 線形ロジスティックス モデルをあてはめた 後に 未知のパラメータの推定値 ˆ β ˆ β が得られ ED5 値は ˆ β ED5 = ˆ β によって推定される ED9 値を得るために p =.9 を 線形ロジスティックス モ デル.9 log = β+ β ED9. で用い その結果 ED9 値は (.97 β) / β であるので ED9 = (.97 ˆ β) / βˆ と推定される 用量ではなく log(dose) が モデルの中で説明変数として使用される場合 これらの 式は 修正を少し必要とする モデル logit( p)= β + β log( d) があてはめられる場合 ED5 値は β + β log( ED5) = となるので ED5 = exp( β/ β) となり ED5 = exp( ˆ β ˆ / β) で推定される 同様に ED9 値は ED9.97 ˆ β exp ˆ β = によって推定される 他の対数の底が モデルの中で使用される場合 この底は この

19 プロビット変換, ロジット変換 3 式の e を代えなければならない コンピューターによる解析では 出力によって さら なる計算を正確に実行するために 対数の底 ( または e) を知ることが必要である ED5 と ED9 値の推定は プロビット回帰モデルのでも同様に得ることができる 例 えば log(dose) が説明変数として使用される場合 モデルは - probit( p)= Φ ( p) = β + β log( d) となる P =.5 とした場合 probit (.5) = であり ロジスティックス モデルと同様に ED5 値は ED5 = exp( ˆ β ˆ / β) によって推定される ここで ˆ β ˆ β はプロビット モデルでのパラメータ推定値である P =.9 のとき probit (.9) =.86 であり ED9 値は exp[(.86 ˆ β ˆ ) / β] によって推定される ED5 のような推定量を得たときに 推定精度を判断することができるように 標準誤差 あるいは信頼区間を与えることが一般に望まれる パラメータ推定関数の標準誤差 対応する信頼区間を得るため用いられる式の展開は 重要であり単純ではないので ここで詳しく示す 推定されたED5 値の近似の標準誤差は 最初にセクション.6 で誘導され この標準誤差に基づいた真のED5 値の近似の信頼区間が与えられる 真の ED5 値のより正確な区間は セクション.7 節で記述されるアプローチ フィラーの定理を使用して求められる.6. 推定された有効量についての標準誤差の近似 ED5 が例えば つのパラメータの関数 g = ( ˆ β ˆ, β) であるばあい 関数の分散の近似の標準的な解法 テイラー級数近似法 または デルタ法により 推定値の標準誤差を得るために使用できる この結果は g = ( ˆ β ˆ, β) の分散がほぼ式 (.6) であることを示している g ˆ Var( ˆ ) g Var( ) g g β β Cov( ˆ, ˆ ) ˆ + β ˆ + ˆ ˆ β β β (.6) β β モデルが説明変数として用量を含んでいる場合 g = ( ˆ β ˆ, β) は ˆ β ˆ / β であることから得られる ˆ ρ = ˆ β ˆ / β v ˆ = Var( β ) v ˆ = Var( β ) および v = Cov( ˆ β ˆ, β ) とおくと 式 (.6) は

20 4 プロビット変換, ロジット変換 v Var( ED5) ˆ ρv + ˆ ρ v ˆ β となる 実際には この式での分散および共分散は 正確には わからないのであるが ˆv ˆv および ˆv の近似値が 式の中で用いられるであろう 推定されたED5 値と s.e.( ED5 ) の標準誤差は その後 式 (.7) によって与えられる v ˆ ρv + ˆ ρ v s.e.( ED5) ˆ β / (.7) このような方法による推定値の標準誤差は 対応するパラメータについての近似的な信頼区間を得ることである 例えば ED5 ±.96 s.e.( ED5) を使用して ED5 値の 95% の信頼区間を得ることができる この方法で計算された信頼区間の低い限界は 負になることが起きるかもしれない ゼロと単にこれを取り替えることができるかもしれないが この問題を回避する つの方法はすべて ED5 値の対数を使うことである log( ED5 ) の標準誤差は 標準的な計算法により s.e.{log( ED5)} ( ED5) s.e.( ED5) を得る log( ED5 ) 値の近似的な 95% 信頼区間は log( ED5) ±.96 s.e.{log( ED5)} であり ED5 値の信頼区間は これらの つの範囲の指数をとることにより得られる しかしながら この方法は 信頼区間のより低い限界が負でないことを保証するが 一般的には勧められない これは推定された log (ED5) 値は 対称的に ED5 値自体が分布しにくいからであり 従って 信頼区間を構築するために用いられる正規分布の仮定は合理的とは思われない ロジスティックスのモデルの中で用いられたとき 説明変数がlog (dose) であるとき ED5 値は exp( ˆ β ˆ / β) によって推定され そして 式 (.7) の関係から log( ED5 ) の標準誤差は 式 (4.3) によって与えられる vˆ ˆ ρvˆ + ˆ ρ vˆ s.e.{log( ED5)} ˆ β / (.8) ED5 自体の標準誤差は 式 (.9) を用いて得られる

21 プロビット変換, ロジット変換 5 s.e.( ED5) ED5 s.e.{log( ED5)} (.9) log (dose) が説明変数として使用される場合 ED5 値の対数の限界の指数をとって ED5 値の信頼区間を計算するほうがよく 式 (.9) に基づいた対称な間隔を使用する のではなく 式 (.8) を使用することが望ましい 同様の手順は ED9 値のような他の多量の興味がある量の標準誤差および対応する 信頼区間を得るために用いられる.7. フィラーの定理を用いた有効用量の信頼区間 フィラーの定理は つの正規分布の確率変数の比率の信頼区間によって得ることが できる一般的な計算結果である この結果は ED5 値の信頼区間の構成を適用する前 の一般的な用語で最初に与えられる ρ = β / βを考えよう ここで β と β は ˆ β と ˆβ によって推定され その平均が β と β 分散が v と v 共分散が v の正規分布となると仮定される関数 ψ = ˆ β ρβˆ を考えよう そのとき ˆ β と ˆβ が β と β の不偏推定量であるので E( ψ ) = β ρβ = となり ψ の分散は V = Var( ψ ) = v + ρ v ρv (.) で与えられる ˆ β と ˆβ は 正規分布に従うと仮定されるので ψ は 同様に正規分布 に従い ˆ β ρβˆ V は 標準正規分布となる 従って z α が 標準正規分布の上側 α /点であるとしたときに ρ の (-α ) % / 信頼区間は ˆ β ρβˆ z V α / 値のセットとなる

22 6 プロビット変換, ロジット変換 両辺を 乗し 等式とし ˆ β + ρ ˆ β ρβˆ ˆ β z V = α / を与える 式 (.) により V を代入したのちに 式の整理と 次のように ρ に関する 二次方程式が得られる ˆ ˆ ˆ ˆ ( β z vˆ ) (vˆ z ) ( vˆ / ρ + / ββ ρ + β z ) α α α / = (.) この二次方程式の つの根は ρ のための信頼限界を構成する これが フィラーの 結果である この結果を ED5 = β / β の信頼区間を得るために用いるために 式 (.) の ρ を ED5 と置き換える さらに パラメータの推定が 線形のロジスティックス モデルで得られたものならば 大きなサンプルについてのみの近似であるので 分散共分散 および ˆv の近 似を および v の代わりに使用しなければならない v v ˆv ˆv ED5 による二次方程式を書き換えると ( ˆ β z vˆ ) ED5 (vˆ z ˆ β ˆ β ) ED5 + ( ˆ β vˆ z ) α/ α / α / = がえられ この 方程式を標準的な手順により解き ED5 値の (-α ) % の信頼限界 のために次の式を得る vˆ z / ˆ α v ˆ ρ g vˆ ˆˆ ˆ ˆ ˆ ρv ρ v g v ˆ vˆ ± + ˆ β v ED5 = g / (.) ここで ˆ ˆ ˆ ρ = β/ β g z vˆ / ˆ = α β である / 強い用量反応関係があるとき ˆβ は にたいして高度に有意にとなり また ˆ / vˆ β は より極めて大きくなるであろう これらの場合に g は 小さくなる すなわち z α / より有意となるような関連の場合 g はより無視できるようになる gが式 (.) でゼロである場合 ED5 値の信頼限界は 式 (.7) の中で与えられたE5 値の標準誤差の近似に基づくものと一致する

23 プロビット変換, ロジット変換 7 log(dose) が説明変数として使用されている場合 ED5 値の信頼区間は フィラーの 定理を用い log(ed5)= β / β について信頼限界を最初に得ることにより計算でき 次に その値についての間隔の推定限界の結果の指数をとればよい

24 8 プロビット変換, ロジット変換.8. ED5 の信頼区間の計算事例 表. の SAS/PROBIT によるプロビット法の 6 回のニュートン ラフソン法の反復 の結果 ˆ β = 5.73 ˆ β = 6.65 が得られ 5% 致死量は log( ED5 ) =.784 そ の 95% 信頼区間は (.4.348) が出力されている.6 節および.7 節で示したデル タ法による近似 フィラーの公式による 95% 信頼区間を計算例として示す 信頼区間の計算のためには ˆ β および ˆβ の分散と共分散が必要であり SAS/PROBIT の model ステートメントのオプションで covb を指定しておく SASの proc PROBIT の反復計算は 初期の回帰係数 β = β = からスタートして 回帰係数として 分散共分散は β = 5.73 β = Σ= が出力される デルタ法 ( 近似法 ) による信頼区間の計算をしてみよう ˆ β 5.73 ED5 = = = ˆ β 6.65 log( ).78 ˆ β 5.73 ˆ β 6.65 ˆ = = =.78 ρ vˆ ˆ ˆ ρvˆ ˆ + ρ v s.e.{log( ED5)} ˆ β /.534 (.78) ( ) + (.78). = 6.65 =.3389 log( ED5) ±.96 s.e.(log( ED5)) =.78±.3389 = (.3,.343) / これらの数値計算を行うためには図.4 に示したように Excel での使用が便利である この結果は フィラーの式での結果 (.4.348) に対して小さめに推定されている

25 プロビット変換, ロジット変換 9 図.4 ED5 とその信頼区間の近似計算 Excel による計算 次に フィラーの式での計算をしてみよう 対数の ED5 について 95% 信頼区間は 次の 次式の根を求めることにより得られる ( ˆ β z vˆ )log( ED5) (vˆ z ˆ β ˆ β )log( ED5) + ( ˆ β vˆ z ) = α/ α/ α/ それぞれの係数は a= ˆ β z v ˆ = = α / b= (vˆ z ˆ ˆ α ββ) = ( ( 4.383).96 ( 5.44) 6.63) = / c= ˆ β vˆ z α = ( 5.44) = / となる これを 次式の公式に代入すると 95% 信頼区間が得られる b± b 4ac log( ED5) ±.96 s.e.(log( ED5)) = = (.4,.348) a これらの数値計算は煩雑であるので Excel による計算シートを図.5 に示す

26 プロビット変換, ロジット変換 プロビット 回帰係数からEDxx% の推定 信頼区間の計算 ln(dose) Var^Cov^ 切片 beta^= 傾き beta^= y% = 5 ln(ed(y%))^=.783 ED(y%)^= Q_a= ln_l95%^=.465 L95%^= Q_b= ln_u95%^= U95%^=.7735 Q_c= 計算はフィラーの式 次式の根の公式により計算 図.5 ED5 とその信頼区間のフィラーの公式 Excel による計算 ED5 とその信頼区間の計算公式は プロビット法の場合のみならずロジット法の場合もそのまま適用できる JMP のロジスティック回帰係数は ˆ β = 6.5 ˆ β =.59 分散共分散は JMP の行列計算により Σ= が計算されているので デルタ法 ( 近似法 ) による信頼区間の計算は Excel の計算シートに数値を入れ替えることにより log( ED5) ±.96 s.e.(log( ED5)) =.748 ±.343 = (.75,.34) 求められる ロジットの回帰係数からEDxx% の推定 信頼区間の計算 ln(dose) Var^Cov^ 切片 beta^= 傾き beta^= y% = 5 ln(ed(y%))^= s.e.(ln(ed(y%)))^=.349 ED(y%)^= ln_l95%^=.7546 L95%^= 6.67 ln_u95%^= U95%^= 図.6 ED5 とその信頼区間の近似計算 フィラーの公式も同様にロジスティック回帰の場合に適用できる JMP による信頼区 間の計算は フィラーの公式によって計算されている Excel の計算シートに数値を入

27 プロビット変換, ロジット変換 れ替えた結果を図.7 に示す ロジットの回帰係数からEDxx% の推定 信頼区間の計算 ln(dose) Var^Cov^ 切片 beta^= 傾き beta^= y% = 5 ln(ed(y%))^= ED(y%)^= Q_a= 3.8 ln_l95%^=.979 L95%^= Q_b= ln_u95%^= U95%^=.483 Q_c= 図.7 ロジット法のフィラーの公式による 95% 信頼区間 注 ) 計算精度が悪いので 分散の有効数字を増やす 吉村の例では 投与量の常用対数をとり 変換されたプロビットに 重み付き最小 乗法を 回繰り返しにより 次の回帰係数が得られた η probit( ˆ i = pi) = log ( dosei) 死亡率の推定値は ˆ μ = log LD = (( ˆ β 5) / ˆ β ) = (.44 5) / 6.63) =.78 5 と標準正規分布 Φ ( η i 5) によって計算された JMP の行列計算により 分散共分散は ˆ Σ= と推定されるので これらを Excel のフィラーの計算シートから 切るときの X 軸の目盛りを読むことにより得られる 解析的には フィラーの式により 常用対数で (.6,.347) となり 投与用量に変換して (6, 3) が得られる フィ ラーの式は.7 節で詳しく述べるので ここでは結果のみを示した

28 プロビット変換, ロジット変換 プロビット 回帰係数からEDxx% の推定 信頼区間の計算 ln(dose) Var^Cov^ 切片 beta^= 傾き beta^= y% = 5 ln(ed(y%))^= ED(y%)^= Q_a= ln_l95%^=.89 L95%^= Q_b= ln_u95%^= U95%^=.6959 Q_c= 計算はフィラーの式 次式の根の公式により計算 図.8 回反復のプロビット法のフィラーの公式による 95% 信頼区間 小数点 4 桁目が吉村のテキストと一致しない 計算誤差の問題であろう

29 3 3. 最尤法による計算の例 3.. 吉村の例 JMP の行列計算言語による プロビット変換 重み付最小 乗法の繰り返し 回 吉村本の表 5-3 の結果に合わせて計算結果を出しているので 対応を取ることにより 最小 乗法の計算原理を学習してもらいたい 分散共分散行列の計算は 信頼区間の計 算のために追加している /* Probit, yoshimura(988), p6, Weighted linear Regression */ d = [,36,83,47,333,45] ; x = log(d) ; y = [,, 5, 8, 9, ] ; n = [,,,,, ] ; p = y :/ n ; X = J(nrow(X),) X; beta = [-7.68, ] ; show(round(beta,4)) ; eta = X*beta ; show(round(eta,4)) ; A = d x y n p eta ; show(round(a,4)) ; /*** step ***/ pi = Normal Distribution(eta - 5) ; show(round(pi, 4)) ; z = Normal Density(eta - 5) ; show(round(z, 4)) ; w = z^ :/ (pi :* (-pi)) ; show(round(w,4)); nw = n :* w ; show(round(nw,4)); y = eta + (p - pi) :/ z ; show(round(y,4)) ; W = diag(w) ; Beta = Inv(X`*W*X)*X`*W*y ; show(round(beta,4)) ; eta = X*beta ; show(round(eta,4)) ; diff = abs(eta - eta) ; show(round(diff,4)) ; B = X pi z nw y eta diff ; show(round(b,4)) ; eta=eta ; /*** step ***/ pi = Normal Distribution(eta - 5) ; show(round(pi, 4)) ; z = Normal Density(eta - 5) ; show(round(z, 4)) ; w = z^ :/ (pi :* (-pi)) ; show(round(w,4)); nw = n :* w ; show(round(nw,4)); y = eta + (p - pi) :/ z ; show(round(y,4)) ; W = diag(w) ; Beta = Inv(X`*W*X)*X`*W*y ; show(round(beta,4)) ; eta = X*beta ; show(round(eta,4)) ; diff = abs(eta - eta) ; show(round(diff,4)) ; B = X pi z nw y eta diff ; show(round(b,4)) ; res=(y-eta)`*w*(y-eta) ; covb = Inv(X`*W*X)*res ;

30 4 Round(beta, 4):[-7.68,5.5594] Round(A, 4): /* d x y n p eta */ [ , , , , , ] Round(beta, 4):[ ,6.43] Round(B, 4): /* X pi z nw y() eta() diff() */ [ , , , , , ] Round(beta, 4):[-.44,6.63] < ================ 回帰係数 Round(B, 4): /* X pi() z() nw() y() eta() diff() */ [ , , , , , ] covb = Inv(X`*W*X)*res < ================ 分散共分散行列 [ , ]

31 5 3.. SAS の POBIT プロシジャによるプロビット法による計算 プロビット変換 ニュートン ラフソン法 Title 'ANZ9a.sas --5 Y.Takahahsi ' ; data d ; input i dose n y p eta p_hat ; datalines ; ; proc probit data=d log inversecl ; model y/n = dose / dist=normal itprint covb ; output out=out p=p std=std xbeta=xbeta ; run ; proc print data=out ; run ; Probit Procedure Iteration History for Parameter Estimates Iter Ridge Loglikelihood Intercept Log(dose) : Model Information Number of Observations 6 Number of Events 34 Number of Trials 6 Name of Distribution Normal Log Likelihood ================ 回帰係数 ================ Analysis of Parameter Estimates Standard 95% Confidence Chi- Parameter DF Estimate Error Limits Square Pr > ChiSq Intercept <. Log(dose) <. ================ 分散共分散行列 ================ Estimated Covariance Matrix Intercept Log(dose) Intercept Log(dose) Probit Model in Terms of Tolerance Distribution MU SIGMA

32 6 Probit Analysis on Log(dose) Probability Log(dose) 95% Fiducial Limits Probit Analysis on dose Probability dose 95% Fiducial Limits OBS i dose n y p eta p_hat p xbeta std

33 JMP の行列言語によるロジスティック回帰 ロジット変換 ニュートン ラフソン法 /* The linear logistic model, Newton-Raphson --3 Y.Takahashi */ d = [,36,83,47,333,45] ; x = log(d) ; y = [,, 5, 8, 9, ] ; n = [,,,,, ] ; X = J(nrow(x),) x ; p = y :/ n ; pi = (n:/) :/ n ; A = X y n p pi ; show(round(a,3)) ; beta=[, ] ; L = [-E] ; epsilon= ; /*** step << i >> ** */ for (i=, epsilon>e-8, i++, show("iteration"); show(round(i, )) ; u = X`* (y - n :* pi) ; show(round(u,3)) ; v = pi :* (-pi) :* n ; show(round(v,3)) ; V = diag(v) ; Im = X`*V*X ; show(round(im,3)); beta= beta + inv(im)*u ; show(round(beta,3)); inv_im=inv(im) ; show(round(inv_im,3)); eta = X*beta ; pi = exp(eta) :/ ( + exp(eta)) ; B = X eta pi ; show(round(b, 3)) ; L = y`* ln(pi) + (n - y)` * ln (-pi) ; epsilon=abs(( L-L)/L) ; C =L L epsilon ; show(c) ; L=L ; ) ;

34 8 Round(A, 3): [.4.5,.34..5, , , ,.653.5] "iteration" Round(i, ): Round(u, 3):[4,4.] Round(v, 3):[.5,.5,.5,.5,.5,.5] Round(Im, 3): [ , ] Round(beta, 3):[-4.9,6.55] Round(inv_Im, 3): [ , ] Round(B, 3): [ , , , , , ] C:[ ] "iteration" Round(i, ): Round(beta, 3):[-.999, 9.653] "iteration" Round(i, ):3 Round(beta, 3):[-5.499,.6] "iteration" Round(i, ):4 Round(beta, 3):[-6.89,.53] "iteration" Round(i, ):5 Round(beta, 3):[-6.,.53] "iteration" Round(i, ):6 Round(u, 3):[,] Round(v, 3):[.47,.373,.488,.66,.55,.4] Round(Im, 3): [ , ] Round(beta, 3):[-6.,.53] < ================= 回帰係数 Round(inv_Im, 3): [ , ] < ================= 分散共分散行列 Round(B, 3): [ , , , , , ] C:[ e-3]