ロジスティック回帰による推測 (V.9 LOGISTICプロシジャの機能拡張) 東京理科大学工学部経営工学科 浜田知久馬

内容 ロジスティックモデル 最尤法による推定の原理 最尤法による検定の原理 条件付ロジスティック回帰の数理 V9のLOGISTIC の機能拡張 　(STRATA文による　 　　条件付ロジスティック回帰)

ロジステｨック曲線とオッズ イベント発現確率ｐ 1－p2 　1－p1 p2 p1 －∞ X1 X2 +∞

例と一般化 薬剤 － + 計 イベント 5 10 15 95 90 185 100 200 薬剤 － + 計 イベント a c n-･ b d n+･ n･- n･+ n

説明変数が１つの場合 ｘ=0:drug－ ｘ=1:drug＋

説明変数が１つの場合 ｘ=0:drug－ 　　　　　ｘ=1:drug＋

説明変数が１つの場合

likelihood（尤度） 最尤法：β0、β１の値を動かしてＬが最も大きくなるようにする方法 薬剤 － + イベント 5 10 95 90 likelihood（尤度） 尤度（Ｌ）＝モデルの下でデータが得られる確率 最尤法：β0、β１の値を動かしてＬが最も大きくなるようにする方法 MLE：Maximum Likelihood Estimator

西遊記 ひたすら西を目指す．

最尤法 ひたすら尤度山の頂上を目指す．

尤度曲面 尤度 （-2.94,0.75)

対数尤度 対数尤度曲面 （-2.94,0.75)

絨毯爆撃 尤度

尤度山の頂上にいるのは？

尤 山の頂上では傾きは0

対数尤度 （-2.94,0.75)

薬剤 － + イベント a c b d 対数尤度とスコア関数

薬剤 － + イベント 5 10 95 90 対数尤度とスコア関数

薬剤 － + イベント+ a c イベント－ b d 最尤推定量

薬剤 － + イベント 5 10 95 90 説明変数が１つの場合の 最尤推定量 粗オッズ比に一致

ロジスティック回帰のプログラム data data; do drug=0 to 1; do y=0,1; input w @@; do i=1 to w;output;end;end;end; cards; 95 5 90 10 ; proc logistic descending; model y=drug; 薬剤 － + イベント 5 10 95 90

Analysis of Maximum Likelihood Estimates ロジスティック回帰の出力 Analysis of Maximum Likelihood Estimates Parameter DF Estimate Standard Error Wald Chi-Square Pr > ChiSq Intercept 1 -2.9444 0.4588 41.1812 <.0001 drug 0.7472 0.5671 1.7359 0.1877 Odds Ratio Estimates Effect Point Estimate 95% Wald Confidence Limits drug 2.111 0.695 6.416

帰無仮説の表現

[MedStat:002877] Raoのスコア検定をどのように計算するのでしょうか？

[MedStat:002878]浜田 Raoのスコア検定は説明変数が１つのときは帰無仮説の下でのUとその分散Vを 計算して, U**2/Vを カイ２乗分布と比べることで行うことができます。 しかし,説明変数が複数ある場合は行列演算が必要になりますので 手計算は困難です。 SASのPROC　LOGISTIC等の統計ソフトウエアを利用することを お勧めします。

[MedStat:002879] プログラムを作成するスキルがあるため、SASを使わなくても行列計算を 手元でさせることはできます。残念ながらスコア検定の知識をはじめ統 計的な素養が不足して行き詰っているところです。 大学ではSASを使える環境ではあるものの、センターに行かないと使えず SASを使わないでスコア検定を行うことを希望します。

[MedStat:002881]大橋先生 医療関係の研究者が統計計算で時間を費やす必要は ないと思うのですがね。 統計家からの協力を得て、もっと生産的な仕事を された方が世のためです。

[MedStat:002886] 正直先週までRaoのスコア検定というものを全く知りませんでした。 尤度比検定とWald統計量を使おうと思っていたのですが、最尤推 定値が定まらず どうしようかと思案しているときに、研究室の方からスコア検定 を教えていただきました。 最尤推 定値が求まらないのに どうして検定できるのだろうか？

尤度山の頂点から帰無仮説の 離れ具合を測るには？ β0 （２）傾斜角度を 測ってみる． （３）地図で位置を 確認する （1）高度を 測ってみる． 尤度比検定 スコア検定 Wald検定

三蔵法師が尤度山の最高天竺にいる．弟子たちは どれくらい離れているか．

尤度比検定　孫悟空　 觔斗雲でひとっ飛び，如意棒で山の高さを測る． H0： β＝0

　スコア検定　沙悟浄 水を流して 勾配を測る. H0： β＝0

Wald検定　猪八戒 地図を頼りにひたすら掘り進み距離を測る. H0： β＝0

尤度原理に基づく3種類の検定 尤度比検定，Wald検定，スコア検定 例 H0： β＝0 の検定 1)尤度比検定 山の高さの違い 1)尤度比検定　　 　山の高さの違い 2）スコア検定　　 　 β＝ 0における傾きが0に近いか 3) Wald検定　　 　　最尤推定量からの隔たり

薬剤 － + イベント 5 10 95 90 3種類の検定の模式図 スコア 尤度比 Wald

薬剤 － + イベント 50 100 950 900 セル度数が10倍になると

セル度数がk倍になると β0 ↓

母数空間 H1 H0 1 1

薬剤 － + イベント+ a c イベント－ b d 尤度比検定

薬剤 － + イベント+ a c イベント－ b d Wald検定

薬剤 － + イベント a c b d 帰無仮説の下でのU 薬剤＋群における観測イベント数と期待イベント数の差

薬剤 － + イベント+ a c イベント－ b d スコア検定

FREQプロシジャの出力 統計量 自由度 値 p 値 χ 2 乗値 1 1.8018 0.1795 尤度比χ 2 乗値 1.8341 薬剤 － + イベント 5 10 95 90 FREQプロシジャの出力 統計量 自由度 値 p 値 χ 2 乗値 1 1.8018 0.1795 尤度比χ 2 乗値 1.8341 0.1756 連続性補正χ 2 乗値 1.1532 0.2829 Mantel-Haenszel のχ 2 乗値 1.7928 0.1806 φ係数   0.0949 一致係数 0.0945 Cramer の V 統計量

Testing Global Null Hypothesis: BETA=0 LOGISTICの3種類の検定の出力 薬剤 － + イベント 5 10 95 90 Testing Global Null Hypothesis: BETA=0 Test Chi-Square DF Pr > ChiSq Likelihood Ratio 1.8341 1 0.1756 Score 1.8018 0.1795 Wald 1.7359 0.1877

LOGISTICのSTRATA文 V.9からSTRATA文が追加 層，マッチングを行った場合の条件付の推測 （条件付ロジスティック回帰) PHREGのDISCRETEオプションによる解析と 等価 EXACT文と組み合わせて正確な解析も可能

1:1マッチングを行った ケース・コントロール研究 曝露 非曝露 Case Cont Cont E- E+ 計 Case 48 4 52 12 16 28 60 20 80 Cont E- E+ 計 Case a b c d N

McNemar検定 曝露 E+ 非曝露 E- Cont E- E+ 計 Case a b c d N

片側P値=Pr(12)+Pr(13)+Pr(14)+Pr(15)+Pr(16)=0.0384 2項分布　Bin(n=16,p=0.5) 片側P値=Pr(12)+Pr(13)+Pr(14)+Pr(15)+Pr(16)=0.0384 両側P値= 0.0384×2=0.0768

Bin(16,0.5)の正規近似 片側P値=0.0227 両側P値=0.0455

Cont E- E+ Case　E- a b Case　E+ c d McNemar検定

FREQプロシジャによるMcNemar検定 proc freq data=cc; tables case*control/agree; exact agree; Cont E- E+ Case 48 4 12 16

FREQプロシジャによる McNemar検定の出力 統計量 (S) 4.0000 自由度 1 Pr > S ( 漸近 ) 0.0455 Pr >= S ( 正確 ) 0.0768

ケース・コントロール研究データの2×2の分割表(曝露×疾患)での集計 Cont E- E+ 計 Case 48 4 52 12 16 28 60 20 80 E － + 計 Case 52 28 80 Cont 60 20 112 48 160

条件付きでない解析のプログラム proc logistic data=cc3 descending; class strata; model response=exposure strata; E－ E+ 計 Case 52 28 80 Cont 60 20 112 48 160 strata 1,･･･,80 マッチング した層

Analysis of Maximum Likelihood Estimates 95% Wald Confidence Limits 条件付きでない解析の出力 Analysis of Maximum Likelihood Estimates Parameter   DF Estimate Standard Error WaldChi-Square Pr > ChiSq Intercept 1 -0.6592 0.2944 5.0135 0.0252 exposure 2.1972 0.8165 7.2417 0.0071 strata 0.6592 1.4271 0.2133 0.6442 ・・・ 2 Odds Ratio Estimates Effect Point Estimate 95% Wald Confidence Limits exposure 9.000 1.817 44.591 strata 1 vs 80 3.000 0.040 223.087 ・・・ 正しいオッズ比：3 観測値：160 母数：81

条件付きの解析のモデル 非曝露 曝露 曝露なし 1 1 1 0 0 1 0 0 曝露あり 0 0 0 1 1 0 1 1 pix：疾患を発症する確率 i：マッチした層(i=1,･･･,80) x：曝露の有無 　　　 case cont　　 曝露なし 1　1 1 0 0 1 0 0 曝露あり 0 0 0 1 1 0 1 1 　　 48(a) 4(b) 12(c) 16(d)

のパターンが得られる条件付確率 caseのみ曝露をうける確率 1 １ ２ 1 １ ２ case E+　 cont E- case E-　 cont E+

条件付ロジスティック回帰の尤度と最尤推定量

条件付きロジスティック回帰のプログラム proc logistic descending data=cc3; class strata; model response=exposure; strata strata; exact exposure /estimate=both outdist=out;

条件付きロジスティック回帰の結果 Testing Global Null Hypothesis: BETA=0 Test Chi-Square DF Pr > ChiSq Likelihood Ratio 4.1860 1 0.0408 Score 4.0000 0.0455 Wald 3.6208 0.0571 Analysis of Maximum Likelihood Estimates Parameter DF Estimate Standard Error Wald Chi-Square Pr > ChiSq exposure 1 1.0986 0.5774 3.6208 0.0571

95% Wald Confidence Limits 条件付きロジスティック回帰の結果 Odds Ratio Estimates Effect Point Estimate 95% Wald Confidence Limits exposure 3.000 0.968 9.302 Exact Odds Ratios Parameter Estimate 95% Confidence Limits p-Value exposure 3.000 0.909 12.762 0.0768

正確な推測 48 4 12 16 層を固定したものでの正確な条件付分布 ・・・ ー 16 ー 1 15 ー 2 14 Cont E- E+ 1 １ ２ 1 １ ２ 層を固定したものでの正確な条件付分布 ー 16 ー 1 15 ー 2 14 Cont E- E+ Case 48 4 12 16 ・・・ ー 14 2 ー 15 1 ー 16

OBS Ｂ Ｃ Score Prob 1 16 16.00 0.00002 2 15 12.25 0.00024 3 14 9.00 0.00183 4 13 6.25 0.00854 5 12 4.00 0.02777 6 11 2.25 0.06665 7 10 1.00 0.12219 8 9 0.25 0.17456 0.00 0.19638 17

スコアカイ2乗の 正確な分布とカイ２乗近似 P値=Pr(4)+Pr(6.25)+Pr(9)+Pr(12.25)+Pr(16)=0.0768 確率 P値=Pr(4)+Pr(6.25)+Pr(9)+Pr(12.25)+Pr(16)=0.0768

mid-p型の信頼区間の計算プログラム proc logistic descending data=cc3; class strata; model response=exposure; strata strata; exact exposure /estimate=both outdist=out cltype=midp;

Exact Parameter Estimates mid-p型の信頼区間の出力 Exact Parameter Estimates Parameter Estimate 95% Confidence Limits p-Value Type exposure 1.0986 0.004279 2.3773 0.0490 MidP(0.5) Exact Odds Ratios Parameter Estimate 95% Confidence Limits p-Value Type exposure 3.000 1.004 10.776 0.0490 MidP(0.5) exactP値=Pr(4)+Pr(6.25)+Pr(9)+Pr(12.25)+Pr(16)=0.0768 midP値=0.5×Pr(4)+Pr(6.25)+Pr(9)+Pr(12.25)+Pr(16)=0.0490

モンテカルロシミュレーションによる近似検定のプログラム proc logistic descending data=cc3 exactoptions (method=networkmc seed=4989 n=20000); class strata; model response=exposure; strata strata; exact exposure /estimate=both outdist=out cltype=exact;

モンテカルロシミュレーションによる近似検定の出力 Exact Parameter Estimates Parameter Estimate 95% Confidence Limits p-Value Type exposure 1.1151 -0.0885 2.7066 0.0745 Exact Exact Odds Ratios Parameter Estimate 95% Confidence Limits p-Value Type exposure 3.050 0.915 14.978 0.0745 Exact

結果のまとめ 条件無 条件付 正確法 カイ2乗 8.000 4.000 ｐ値　 0.003 0.046 0.077 オッズ比 9.000 3.000 信頼下限 1.817 0.968 0.909 信頼上限 44.591 9.302 12.762

尤 まとめ スコア 尤度比 Wald H0： β＝0

スコア検定の利点 ・ MLEを求めるためには反復計算が必要 ・スコア検定はＨ0の下でのＵがわかれば計算可能 ・Wald，尤度比検定はMLEが求まらないとできない． ・ MLEを求めるためには反復計算が必要 ・スコア検定はＨ0の下でのＵがわかれば計算可能 ・多くのモデルについての計算が必要な総当り法では,スコア検定が行われる． ・単純な問題については，よく知られた検定に一致 ・スコア検定では収束しない場合でも,Ｈ0の検定が可能

次のうちスコア検定に相当するのはどれでしょう． 1)Z検定(分散既知のｔ検定) 2)Pearsonのカイ2乗検定 3)McNemar検定 4)Cochran-Armitage検定 5)Mantel-Haenzel検定 6)ログランク検定

[MedStat:002887]浜田 スコア検定は確かに 最尤推定値が求まらなくても 帰無仮説が検定できるのが 利点です. しかし最尤推定値が求まらないのは モデルが破綻しているということですし 帰無仮説の検定だけでは推測としては不十分です. 根本的に最尤推定値が求まらない原因(0セルがある等）を追究しとく 必要があるかと思います. やはり統計の専門家に相談した方がよいと思います.

参考文献 Derr, R.E.(2000) Performing exact logistic　regression with the SAS System. SUGI'2000 Proceedings, Paper 254 Gail, M.H., Lubin, J.H., and Rubinstein, L.V. (1981)　Likelihood Calculations for Matched Case-Control Studies and Survival Studies with Tied Death Times. Biometrika, 68, 703-07. Hirji, K.F., Mehta, C.R., and Patel, N.R. (1987) Computing Distributions for Exact Logistic Regression. Journal of the American Statistical Association, 82, 1110 - 1117. Hosmer, D.W, Jr. and Lemeshow, S. (2000), Applied Logistic Regression, Second Edition, New York: John Wiley & Sons, Inc. Mehta, C.R., Patel, N. and Senchaudhuri, P. (1992), Exact Stratified Linear Rank Tests for Ordered Categorical and Binary Data. Journal of Computational and Graphical Statistics, 1, 21 - 40. Mehta, C.R., Patel, N. and Senchaudhuri, P. (2000) Efficient Monte Carlo Methods for Conditional Logistic Regression. Journal of the American Statistical Association, 95, 99 - 108. Truett，J., Cornfield, J. and Kannel, W.(1967)　A Multivariate Analysis of the Risk of Coronary Heart Disease in Framingham. J.Chron.Dis. 20, 511-524 浜田知久馬(1994)ＳＡＳによる条件付きロジスティック回帰. 日本ＳＡＳユーザー会94論文集,527-540 浜田知久馬(2000)LOGISTICのV. 8の機能拡張. 日本ＳＡＳユーザー会2000論文集,13-38 浜田知久馬(2001)SAS V. 8における正確な推測とシミュレーションによる近似法. 日本ＳＡＳユーザー会2001論文集,165-187

data likelihood; do b0=-2.9444; do b1=-0.4 to 1.2 by 0.05; p0=1/(1+exp(-(b0))); p1=1/(1+exp(-(b0+b1))); l=p0**10*(1-p0)**190*p1**20*(1-p1)**180; logl=log(l); output; end;end; proc gplot; plot l*b1; symbol1 i=spline ;run; plot logl*b1;

理想的な比較

data data;phi=0.50;n=16; do y=0 to 16; p=pdf('binomial',y,phi,n); fn=pdf('normal',y,n*phi,(n*0.25)**.5); output;end; proc gplot;plot p*y fn*y/vzero overlay; symbol1 i=needle c=red; symbol2 i=spline c=green; run;

data data; b=4;c=12; do beta=0 to 2 by 0.1; logl=c*beta-(b+c)*log(1+exp(beta)); output; end; proc gplot; plot logL*beta; symbol1 i=spline; run;

UNIVARIATEプロシジャによるMcNemar検定 data cc; input case control w @@; do i=1 to w; dif=case-control;output;end; cards; 0 0 48 1 1 16 0 1 4 1 0 12 ; proc univariate data=cc;var dif; case control dif 0 0 0 1 0 1 0 1 -1 1 1 0

UNIVARIATEプロシジャによる 対応のある検定の出力 位置の検定 : μ 0=0 検定 統計量 p 値 Student の t 統計量 t 2.039 Pr > |t| 0.0448 符号検定 M 4 Pr >= |M| 0.0768 符号付順位検定 S 34 Pr >= |S|

マッチングを無視した解析のプログラム data cc2; input response exposure w @@; 　　do i=1 to w;output;end; 　cards; 　1 0 60 1 1 20 0 0 52 0 1 28 　; 　proc logistic data=cc2 descending; 　model response=exposure;

Analysis of Maximum Likelihood Estimates 95% Wald Confidence Limits マッチングを無視した解析の結果 Analysis of Maximum Likelihood Estimates Parameter DF Estimate Standard Error Wald Chi-Square Pr > ChiSq Intercept 1 -0.1431 0.1895 0.5705 0.4501 exposure 0.4796 0.3487 1.8912 0.1691 Odds Ratio Estimates Effect Point Estimate 95% Wald Confidence Limits exposure 1.615 0.816 3.200 正しいオッズ比：3 層を無視するため誤差的なバラツキが増大

正確な推測 層を固定したものでの正確な条件付分布 Cont E- E+ Case a B C d

結果のまとめ ロジ(1) 層無視 ロジ(2) 層条件無 条件付 正確法 カイ2乗 1.905 8.000 4.000 ｐ値　 0.122 0.003 0.046 0.077 オッズ比 1.615 9.000 3.000 信頼下限 0.816 1.817 0.968 0.909 信頼上限 3.200 44.591 9.302 12.762

良性乳癌に対する1:3マッチングを 行ったケースコントロール研究 ケース: 50人　　コントロール:150人 STR：(層)　１～50　年齢と施設でマッチング AGMT(面接時年齢)　FNDX(乳癌の有無) CHK(定期的な診断の有無) AGMN(初経年齢) HIGD(就学期間)　DEG(学歴) NLV(死産児の数) LIV(生誕時の数) WT(体重ポンド) AGLP(閉経時の年齢) MST2(結婚歴)　1:婚姻歴有　2:婚姻歴無

変数減少法による変数選択 proc logistic descending; class str mst2/param=ref ref=last; model fndx= chk|deg|higd|agmn|aglp|wt|mst2@2 /selection=backward; strata str;

変数減少法による変数選択 定期的な診断 初経年齢 体重 結婚歴 Summary of Backward Elimination Step Effect Removed DF Number In Wald Chi-Square Pr > ChiSq 1 HIGD 6 2.3863 0.1224 2 DEG 5 1.7836 0.1817 3 AGLP 4 3.0579 0.0803 Type 3 Analysis of Effects Effect DF Wald Chi-Square Pr > ChiSq CHK 1 6.7503 0.0094 AGMN 7.8913 0.0050 WT 8.0069 0.0047 MST2 4.6865 0.0304 定期的な診断 初経年齢 体重 結婚歴

パラメータ推定値 Analysis of Maximum Likelihood Estimates Parameter DF   DF Estimate Standard Error Wald Chi-Square Pr > ChiSq CHK 1 -1.1613 0.4470 6.7503 0.0094 AGMN 0.3592 0.1279 7.8913 0.0050 WT -0.0282 0.00998 8.0069 0.0047 MST2 -1.5934 0.7360 4.6865 0.0304 CHK:定期的な診断　AGEN:初経年齢　WT:体重　MST2:結婚歴

条件付の解析 Odds Ratio Estimates Effect Point Estimate 95% Confidence Limits CHK 0.313 0.130 0.752 AGMN 1.432 1.115 1.840 WT 0.972 0.953 0.991 MST2 1 vs 2 0.203 0.048 0.860 定期的な診断 初経年齢 体重 結婚歴 条件無しの解析　Odds Ratio Estimates Effect Point Estimate 95% Confidence Limits CHK 0.185 0.064 0.536 AGMN 1.741 1.272 2.382 WT 0.959 0.935 0.983 MST2 1 vs 2 0.0840 0.014 0.512