ロジスティック回帰による推測と多値データへの拡張)
東京理科大学工学部経営工学科
浜田知久馬
1

内容
ロジスティックモデル
最尤法による推定の原理
最尤法による検定の原理
比例オッズモデル
一般化ロジットモデル

2ロジステｨック曲線とオッズ
イベント発現確率ｐ1−p2

1−p1podds =1− pp2

p1−∞

exp(β0 + β1x)p=
1 + exp(β0 + β1x)

X1

X2

⎛ p ⎞
⎟⎟ = β0 + β1xlog⎜⎜⎝ 1− p ⎠+∞3

例と一般化薬剤−
イベント 5+イベント 95
−100計

薬剤 計+
101590185100200薬剤−
イベント a+イベント b−n･計

薬剤 計+cn-･
dn+･n･+n4

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

5説明変数が１つの場合

ｘ=0:drug−

ｘ=1:drug＋

6

説明変数が１つの場合

7

薬剤
−

薬剤+

イベント
+

5

10

イベント
−

95

90

likelihood（尤度）

尤度（Ｌ）＝モデルの下でデータが得られる確率

最尤法：β0、β１の値を動かしてＬが最も大きくなるようにする方法
MLE：Maximum Likelihood Estimator

8

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

9

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

10

尤度曲面

尤度

（-2.94,0.75)

β0

β1
11

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

β0

β1

12

絨毯爆撃
尤度

β0

β1
13

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

14

山の頂上では傾きは0∧β15対数尤度（-2.94,0.75)

β0

β1

16

薬剤
−

薬剤+

イベント
+ac

イベント
−bd

対数尤度とスコア関数17

薬剤−

薬剤+イベント
+510

イベント
−

9590

対数尤度とスコア関数

18

薬剤
−

薬剤+イベント+ac

イベント−bd

最尤推定量

19

薬剤
−

薬剤+イベント
+510

イベント
−9590

説明変数が１つの場合の最尤推定量

粗オッズ比に一致
20ロジスティック回帰のプログラム
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

10959021

ロジスティック回帰の出力
Analysis of Maximum Likelihood Estimates
Parameter

Interceptdrug

DF

Estimate StandardError∧

WaldChiSquarePr > ChiSq

β1 0 -2.9444

0.458841.1812<.0001

β1 1 0.7472

0.5671

1.73590.1877∧

Odds Ratio Estimates
Effect

Point
Estimatedrug

2.111

95% Wald
Confidence
Limits
0.695 6.416

22帰無仮説の表現

23

[MedStat:002877] Raoのスコア検定をどのように計算するのでしょうか？
ロジスティック回帰分析
で係数の決定には最尤推定法を用います。
Raoのスコア検定を使おうと考えていますが，どう算出 したらよいかわからず悩んでおります．
24

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

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

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

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

28

尤度山の頂点から帰無仮説の離れ具合を測るには？
β0∧

β（1）高度を測ってみる．
尤度比検定（２）傾斜角度を測ってみる．

（３）地図で位置をスコア検定

Wald検定 29

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

尤度比検定 孫悟空

觔斗雲でひとっ飛び，如意棒で山の高さを測る．

H0： β＝031

スコア検定 沙悟浄

水を流して勾配を測る.

H0： β＝032

Wald検定 猪八戒

地図を頼りにひたすら掘り進み距離を測る.

H0： β＝033

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

薬剤
−薬剤
+イベント
+510

イベント
−9590

3種類の検定の模式図

尤度比

スコア

β1Wald35

薬剤
−薬剤
+

イベント
+50

100

イベント
−950

900

セル度数が10倍になると

β136

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

∧β↓

37

H0
1

π+0

母数空間

H1

H0 : β1 = 0
H0 :π− =π+

β0

π−

1

β1
38

β0

∧

0, β 0 H 0∧∧β1 , β 0 H

1

β139

β0

β140

薬剤
−

薬剤
+

イベント+ac

イベント−bd

尤度比検定

41

薬剤
−

薬剤
+

イベント+

ac

イベント−bd

Wald検定

42

β0

β143

薬剤
−

薬剤
+

イベント
+ac

イベント
−bd

帰無仮説の下でのU

薬剤＋群における観測イベント数と期待イベント数の差

44

薬剤−

薬剤+イベント+ac

イベント−bd

スコア検定

45

薬剤−

薬剤+

イベント+510

イベント
−9590

FREQプロシジャの出力

統計量

自由度値

p 値

χ 2 乗値

11.8018 0.1795

尤度比χ 2 乗値

11.8341 0.1756

連続性補正χ 2 乗値

11.1532 0.2829

Mantel-Haenszel のχ 2 乗値

11.7928 0.1806

φ係数0.0949

一致係数0.0945

Cramer の V 統計量0.0949
46

LOGISTICの3種類の検定の出力
薬剤−

薬剤+イベント+510

イベント
−9590

Testing Global Null Hypothesis: BETA=0Test
Chi-Square DF Pr > ChiSq
Likelihood Ratio 1.8341
10.1756
Score
1.801810.1795
Wald
1.735910.1877
47

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

ロジスティックモデルの多値データへの拡張
・名義尺度
一般化ロジットモデル(generalized logit model)
・順序尺度
比例オッズモデル(累積オッズモデル)(cumulative logit model)49一般化ロジットモデル

50

一般化ロジットモデル

51

比例オッズモデル

52

比例オッズモデル

53

比例オッズモデル

54

比例オッズモデル

55

ワニの体長と餌の選択
data gator;
input length choice $ @@;cards;
1.24 I 1.30 I 1.30 I 1.32 F
1.45 I 1.45 O 1.47 I 1.47 F1.63 I 1.65 O 1.65 I 1.65 F
1.78 I 1.78 I 1.78 O 1.80 I
1.98 I 2.03 F 2.03 F 2.16 F
2.36 F 2.39 F 2.41 F 2.44 F
2.79 F 2.84 F 3.25 O 3.28 O
3.68 O 3.71 F 3.89 F;

1.32 F1.50 I
1.65 F1.80 F2.26 F2.46 F3.33 F1.40 F1.52 I
1.68 F1.85 F2.31 F
2.56 O
3.56 F1.42 I
1.55 I
1.70 I
1.88 I
2.31 F2.67 F3.58 F1.42 F1.60 I
1.73 O
1.93 I
2.36 F2.72 I
3.66 F56

一般化ロジットモデル
proc logistic;
model choice = length / L=GLOGIT;
output out=out p=p;proc gplot;
plot p*length=̲level̲;
symbol1 i=spline c=red w=3;
symbol2 i=spline c=green w=3;
symbol3 i=spline c=blue w=3;57

一般化ロジットモデル
Type 3 Analysis of Effects
EffectlengthDF2Wald
Chi-Square
8.9360Pr > ChiSq
0.0115

Analysis of Maximum Likelihood Estimates

Parameter
Intercept
Intercept
lengthlengthchoice
FIFIF

Estimate
1.6177
5.6974-0.1101-2.4654

IStandard
WaldError
Chi-Square
Pr > ChiSq
1.3073
1.53140.21591.793810.08810.00150.51710.04530.83140.89977.51010.0061

O
58

59

比例オッズモデル
data gator;set gator;
if choice='I' then y=1;
if choice='F' then y=2;
if choice='O' then y=3;
proc logistic;
model y= length/l=logit;
output out=out p=p;proc gplot;
plot p*length=̲level̲;
symbol1 i=spline c=red w=3;
symbol2 i=spline c=green w=3;
symbol3 i=spline c=blue w=3;
run;
60

比例オッズモデル
Score Test for the Proportional Odds Assumption
Chi-Square
4.9162

DF1

Analysis of Maximum Likelihood Estimates
StandardParameter
DFEstimate
Error
Intercept 111.88680.8247
Intercept 2
14.8765
1.0929length1-1.2420
0.3953Pr > ChiSq
0.0266Wald
Chi-Square
5.234519.9098
9.8698Pr > ChiSq
0.0221<.00010.0017

61

Others

FishInsect62

比例オッズ性の検定

63

まとめスコア

尤度比Wald∧β

H0： β＝064

スコア検定の利点
・Wald，尤度比検定はMLEが求まらないとできない．

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

65

次のうちスコア検定に相当するのはどれでしょう．
1)Z検定(分散既知のｔ検定)
2)Pearsonのカイ2乗検定
3)McNemar検定4)Cochran-Armitage検定5)Mantel-Haenzel検定
6)ログランク検定66[MedStat:002887]浜田
スコア検定は確かに 最尤推定値が求まらなくても 帰無仮説が検定できるのが 利点です. しかし最尤推定値が求まらないのは モデルが破綻しているということですし 帰無仮説の検定だけでは推測としては不十分です. 根本的に最尤推定値が求まらない原因(0セルがある等）を追究しとく 必要があるかと思います. やはり統計の専門家に相談した方がよいと思います.
67

参考文献
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 CaseControl 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 LogisticRegression. 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 forOrdered 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
68