玄関雑学の部屋雑学コーナー統計学入門

10.3 ロジスティック回帰分析の計算方法

(1) 最小2乗法を利用する方法

ロジスティック回帰分析ではロジットの回帰誤差が特殊な分布になり、普通はその分布を理論的に確定することができません。 そのため回帰誤差が近似的に正規分布すると見なして、重回帰分析と同じように最小2乗法を利用して回帰分析を行う方法と、最尤法を利用した繰り返し近似計算によって回帰分析を行う方法の2種類があります。

最初に最小2乗法を利用する方法について説明しましょう。 この場合の線形ロジスティックモデルは次のようになります。 このモデルでは説明変数xjの値を研究者が任意に指定した時、ロジットの回帰誤差εが近似的に正規分布すると仮定します。 ただしロジスティック回帰式を計算するには回帰誤差の正規性は必要ではなく、検定を行う時だけ回帰誤差の正規性が必要になります。 そして説明変数は研究者が任意の値を指定するため誤差がなく、正規性という概念そのものが当てはまりません。 重回帰分析と同様に、これはよく誤解されていることです。

線形ロジスティックモデル:
η:出現率πのロジット  β0:定数  βj:偏回帰係数(j=1,…,p)  ε:誤差、近似的に正規分布N(0,σ2)に従う

このモデルは目的変数が少し特殊なだけで、重回帰モデルとほとんど同じです。 そのため計算方法も重回帰分析とほとんど同じです。 例として表10.3.1のデータにこのモデルを当てはめ、最小2乗法を利用して解を求めると次のようになります。 (注1)

表10.3.1 ロジスティック回帰分析用
テストデータ
説明変数x反応無反応有出現率p
10282300.067
20314350.114
303314470.298
40813210.619
50639450.867
合計106721780.404
線形ロジスティックモデル:
ロジスティック回帰式:l=-4.378 + 0.122x   ロジスティック曲線:
寄与率:r2=0.983(98.3%)   回帰係数の95%信頼区間:下限 βL=0.088 上限 βU=0.155
直線性(回帰)の検定:χβ2=50.789(p=1.0280×10-12)>χ2(1,0.05)=3.841 … 有意水準5%で有意
異質性(ズレ)の検定:χLOF2=0.854(p=0.8365)<χ2(3,0.05)=7.815 … 有意水準5%で有意ではない

重回帰分析と同様にロジスティック回帰分析も記述統計学的手法のため、推測統計学的手法である検定とは相性が良くありません。 それでも有意症患者のために、回帰誤差が近似的に正規分布するという性質を利用して検定を行うことができます。

直線性の検定は出現率の変動の中にロジスティック曲線で説明できる変動があるかどうか、つまり寄与率が0かどうかの検定です。 異質性(ズレ)の検定は出現率の変動の中にロジスティック曲線では説明できない変動があるかどうか、つまり実際の出現率とロジスティック曲線のズレが0かどうかの検定です。 ただしこれらの検定で検出差を指定するのは困難なので、たいていの場合は有意性検定になります。 そのため実質的な意味はほとんどないと考えた方が良いと思います。

ちなみに表10.3.1のデータにコクラン・アーミテージの傾向分析を適用し、説明変数xと出現率pを直線で回帰すると次のようになります。 そしてこの回帰直線と最小2乗法を利用して計算したロジスティック回帰曲線をグラフ化すると、図10.3.1の黒色の直線青色の曲線になります。

回帰モデル:p=α + βx + ε
直線回帰式:p=-0.262 + 0.0216x
寄与率:r2=0.945(94.5%)   回帰係数の95%信頼区間:下限 βL=0.016 上限 βU=0.027
直線性(回帰)の検定:χβ2=68.573(p=2.3037×10-15)>χ2(1,0.05)=3.841 … 有意水準5%で有意
異質性(ズレ)の検定:χLOF2=4.014(p=0.2599)<χ2(3,0.05)=7.815 … 有意水準5%で有意ではない

また表10.3.1のデータをプロビット変換してから回帰分析を行うと次のようになります。 そしてこの変数−プロビット直線を変数−出現率曲線に変換すると、図10.3.1の赤色の曲線になります。 表10.3.1の反応有の閾値分布が近似的に正規分布するなら、ロジスティック回帰曲線よりもこのプロビット回帰曲線の方が正確です。

回帰モデル:y(pのプロビット)=α + βx + ε
変数−プロビット直線回帰式:y=2.445 + 0.071x
寄与率:r2=0.972(97.2%)
直線性(回帰)の検定:χβ2=58.857(p=1.7347×10-14)>χ2(1,0.05)=3.841 … 有意水準5%で有意
異質性(ズレ)の検定:χLOF2=1.680(p=0.6414)<χ2(3,0.05)=7.815 … 有意水準5%で有意ではない
図10.3.1 回帰直線と回帰曲線

(2) 最尤法を利用する方法

最尤法を利用する方法は回帰誤差が近似的に正規分布しなくても適用できます。 その代わり厳密な解を代数的に求めることができないため、適当な初期値を用いた繰り返し計算によって近似解を求めます。 そのため初期値と繰り返し計算の収束条件を変えると計算結果が変わってしまいます。 したがって使用する統計ソフトや計算条件によって結果が変わることがあり、最小2乗法を利用する方法と違って誰が計算しても同じ結果になるというわけではありません。

判別分析も最尤法によって判別関数を求めます。 しかし判別分析の場合は2群を研究者が任意に指定し、それらの群の検査項目の値が近似的に正規分布すると仮定して最尤法を適用します。 つまり目的変数である群の出現率には誤差がなく、説明変数である検査項目に誤差があると仮定して最尤法を適用するわけです。 それに対してこの場合の最尤法は説明変数には誤差がなく、目的変数であるロジットまたは出現率に誤差があると仮定して最尤法を適用します。

表10.3.1のデータについて、最尤法を利用する方法でロジスティック回帰分析を行うと次のようになります。 ここでは最小2乗法を利用して計算した解を初期値にしたため、結果はほとんど変わりません。 したがってロジスティック回帰曲線もほとんど変わらず、図10.3.1の青色の曲線と重なってしまいます。 (注2)

ロジスティック回帰式:l=-4.446 + 0.124x   ロジスティック曲線:
寄与率:r2=0.982(98.2%)   回帰係数の95%信頼区間:下限 β1L=0.089 上限 β1U=0.158
回帰係数の検定:χβ12=48.759(p=2.8938×10-12)>χ2(1,0.05)=3.841 … 有意水準5%で有意
直線性(回帰)の検定:χβ2=79.369(p=2.5952×10-15)>χ2(1,0.05)=3.841 … 有意水準5%で有意
異質性(ズレ)の検定:χLOF2=0.783(p=0.8535)<χ2(3,0.05)=7.815 … 有意水準5%で有意ではない

この計算法では回帰誤差の分布を確定しないため、原則として検定と推定はできません。 しかし何が何でも検定しなければ気がすまない有意症患者のために、最尤法による解が漸近的に正規分布するという性質を利用して、上記のように回帰係数の検定と推定、そして直線性(回帰)の検定と異質性(ズレ)の検定を行うことがあります。

最尤法を利用する方法は回帰誤差が近似的に正規分布しなくても適用できるところに特徴があります。 ところが上記の検定と推定は最尤法による解が近似的に正規分布する性質を利用しています。 最尤法による解が厳密に正規分布するのは回帰誤差が正規分布する時であり、その時は最小2乗法による解と最尤法による解が一致し、各種の検定結果も一致します。 したがってどうせ近似的な正規性を利用するなら、回帰誤差が近似的に正規分布すると仮定して計算する方が便利です。

また最尤法を利用する方法は適当な初期値を用いた繰り返し計算によって近似解を求め、近似解が収束したところで計算を終了する近似計算法です。 そのためデータが極端に偏っていたり誤差が多かったりすると繰り返し計算が発散して異常な解になったり、解が求められなかったりすることがあります。 例えば表9.1.1のデータを前向き研究で得られたデータと見なしてロジスティック回帰分析を適用すると、次のように判別分析から得られるロジスティック曲線とはかなり異なった結果になります。

l:ロジット  p:動脈硬化の出現率  x1:TC  x2:TG
ロジスティック回帰分析の回帰式:l=-448.932 + 2.270x1 - 0.341x2
ロジスティック曲線:
判別分析の判別関数:z=-61.164 + 0.335x1 - 0.075x2
ロジスティック曲線:

この原因は表9.1.1のデータがきれいすぎる、つまり極端に偏っているからです。 図9.6.1を見るとわかるように、表9.1.1のデータに判別分析を適用すると25例のデータが全て正しく判別されてしまいます。 そのためロジスティック回帰分析では繰り返し計算が発散してしまい、異常な解になってしまったのです。 これらのことからわかるように、最尤法を利用した計算法よりも最小2乗法を利用した計算法の方が色々な意味で便利です。

しかし説明変数が多数になると、全ての説明変数の値が同じというデータが複数あるのは稀です。 そうするとデータを表10.3.1のような形式にまとめた場合、出現率が0または1というものが大多数になります。 出現率が0または1の時、ロジット変換はできません。 そのためこのような場合、最小2乗法を利用した計算法では反応有の例数と合計例数にそれぞれ0.5を加えてからロジット変換した修正値を用いて計算を行います。 そのような計算法では、当然、結果の精度が低くなってしまいます。

それに対して最尤法を利用した方法ではロジット変換を直接利用せずに特殊な方法で計算を行うため、出現率が0または1のデータも修正せずに計算に入れることができます。 そのためコンピュータが発達したことと相まって、現在は最尤法を利用する方法が主流になっています。


(注1) 表10.3.1を一般化し、線形ロジスティックモデルを当てはめると次のようになります。

表10.3.2 説明変数がp個の一般的データ
説明変数反応有反応無出現率ロジット重み
x11x1jx1pr1n1-r1n1p1l1n1w1
:::::::::
xi1xijxiprini-rinipiliniwi
:::::::::
xm1xmjxmprmnm-rmnmpmlmnmwm
        
  

出現率が0または1の時にも通用するように、次のようなロジットを用いる時もあります。

または
または
線形ロジスティックモデル:
  
  

重みniwiは重み付け最小2乗法に用いるもので、ロジットの分散の逆数です。 この重みは誤差の少ないロジットほど重要視するためのものであり、コクラン・アーミテージの傾向検定の重みと同じように、出現率が二項分布することから導き出されます。 (→5.3 計数値の相関と回帰 (注4)7.1 重回帰モデル (注1))

線形ロジスティックモデルの両辺に重みの平方根√(niwi)をかけ、重み付け最小2乗法によってβの最良不偏推定値(BLUE解)を求めると次のようになります。

1/2=1/2β+1/2ε      εi〜N(0,12)
  =(1/2)2
正規方程式とその解:[']='   =[']-1'
V()=σ2[']-1≒[']-1   σ2≒1

回帰誤差εiが近似的に正規分布することを利用して、偏回帰係数の検定と推定を行うことができます。

検定:≧χ2(1,α)の時、有意水準100α%で有意
推定:100(1-α)%信頼区間=
→ 下限:  上限:
[']jj-1:[']-1の第j対角要素

説明変数が1つの時は次のようになります。

  



  

この時は偏回帰係数の検定がそのまま直線性(回帰)の検定になり、異質性(ズレ)の検定も行うことができます。

Sβ=b12Sxxβ2≧χ2(1,α)の時、有意水準100α%で有意
SLOF=Sll-SβLOF2≧χ2(m-2,α)の時、有意水準100α%で有意
寄与率:

表10.3.1の例題について実際に計算してみましょう。

           
           
           


  

  b0=-0.334-33.181×0.122≒-4.378   
bjの95%信頼区間= → 下限:β1L=0.088 上限:β1U=0.155
Sββ2=0.1222×3419.316≒50.789(p=1.0280×10-12)>χ2(1,0.05)=3.841 … 有意水準5%で有意
SLOFLOF2=51.643-50.789=0.854(p=0.8365)<χ2(3,0.05)=7.815 … 有意水準5%で有意ではない

m=2でx1=0、x2=1とすると、次のように直線性の検定がオッズ比の検定に一致します。 これは直線性の検定はロジットの回帰係数つまり対数オッズの差の検定に相当し、回帰係数を指数変換するとオッズ比になるためです。 (→3.4 2標本の計数値 (2)名義尺度(分類データ) (注5))





(注2) 線形ロジスティックモデルに最尤法を適用する時はデータを次のような形式で整理します。

表10.3.3 説明変数がp個の一般的データ(最尤法用)
群分類y
1:反応有・0:反応無
説明変数
x1xjxp
1x11x1jx1p
::::
1xr1xrjxrp
0x(r+1)1x(r+1)jx(r+1)p
::::
0xn1xnjxnp
線形ロジスティックモデル:
  
        

回帰誤差εiが正規分布すると仮定せず、最尤法によってβの最尤推定値を求めます。

反応有の個体を得る確率:
反応無の個体を得る確率:
両者合わせて表すと:
全体の尤度:
対数尤度:

対数尤度を最大にするニュートン・ラプソン(Newton-Raphson)法によって近似的に求めます。 ニュートン・ラプソン法はニュートンの2階勾配法とも呼ばれ、関数を2次式で近似し、関数の最大値または最小値を反復計算によって近似的に求める方法です。

今、最大値(または最小値)を持つ関数をf(x)とします。 一般にf(x)は最大値の近傍では2次関数つまり放物線で近似できます。 そこでf(x)を最大値に近い値xkテーラー展開(Tayler expansion)し、2次の項まで取ります。


f'(xk):f(x)の1次微係数  f''(xk):f(x)の2次微係数

この近似関数を最大にする時のxは、関数をxで微分して0と置いた方程式の解になります。

f(x)の1次微係数と2次微係数が直接求められない時は、次のような差分法によって近似的に求めます。

ニュートンの前進差分:   
スターリングの中心差分:   
ニュートンの後退差分:   
※|Δ|≪0 (|Δ|=xk×10-13〜xk×10-6等を用いる)

この解xk+1のところでまたテーラー展開を行い、次の解を求めれば次第に最大値に近づくはずです。 このような反復計算を繰り返し、xkとxk+1の差が十分に小さくなったところでxkを最終的な近似解にします。 関数の変数が複数の時は行列とベクトルを利用して次のように計算します。


:f()の1次偏微係数ベクトル(傾斜ベクトル)
:f()の2次偏微係数行列(Hessの行列、f()の曲率を表す計量行列)
k+1=k-k-1k   :ナブラ(ハミルトン演算子)

このニュートン・ラプソン法を対数尤度関数に適用すると次のようになります。

  
k=k   
k+1=k-k-1k
図10.3.2 ニュートン・ラプソン法の模式図

1次偏微分ベクトルkはL(k)の傾斜を表す傾斜ベクトルであり、最大値付近ではmax(ゼロベクトル)になります。 そして2次偏微分係数行列k(ヘスの行列)はL(k)の曲率を表す計量行列です。 最大値付近におけるこの曲率が大きいほどmaxの変動幅が小さくなるので、これはmaxの確実性を表す指標になります。

ただしL(k)は上に凸の曲線なので、この曲率は負の値になります。 そこでこの曲率の符号を反転させた値の期待値を求めると、それは推定量の確実性を表す指標つまりデータが持つ情報量の多さを表す指標になります。 それをフィッシャーの情報量(Fisher's information)と呼び、maxの符号を反転させた行列の期待値E(-max)をfと書いてフィッシャーの情報行列(Fisher's information matix)と呼んでいます。

最尤推定値には漸近的正規性があり、フィッシャーの情報量の逆数が最尤推定量の分散になります。 そこでこれらの性質を利用して偏回帰係数の検定と推定を行うことができます。 この検定をワルド(Wald)の検定といいます。 (→1.4 推定 (注4))

E(-max)=f:フィッシャーの情報行列
V(bj)=[f]jj-1f-1≒[-max]-1の第j対角要素=bjの分散
検定:≧χ2(1,α)の時、有意水準100α%で有意
推定:100(1-α)%信頼区間
→ 下限:  上限:

説明変数がなく、定数項だけの線形ロジスティックモデルの尤度と対数尤度は次のようになります。



r:反応有の例数  n:例数

尤度は「もっともらしさ」ですから、この最も単純なモデルの尤度と説明変数がp個の時の尤度の比が大きいほど説明変数の寄与が少なく、尤度比が小さいほど寄与が大きいことになります。 (→9.3 1変量の場合 (1) 尤度と最尤法)

尤度比の逆数を平方して対数変換した値、つまり対数尤度差に(-2)を掛けた値が近似的にχ2分布することを利用した検定法を尤度比検定(likelihood ratio test)といいます。 この尤度比検定を定数項だけのモデルと説明変数がp個のモデルの対数尤度の差に適用すれば、説明変数全体の回帰の検定つまり直線性の検定を行うことができます。 またモデルの当てはまりの良し悪しを表す指標として、赤池の情報量基準AIC(Akaike's Information Criterion)を計算することもできます。 (→(2) 名義尺度(分類データ) 1) データに対応がない場合9.3 1変量の場合 (注1))

対数尤度差:
尤度比検定:-2{L0 - L(β)} = χβ2≧χ2(p,α)の時、有意水準100α%で有意
AIC = 2{(p+1) - L(β)}

一方、説明変数と定数項を合わせた数(p+1)が例数nと同じ時、線形ロジスティックモデルとデータが完全にフィットします。 これを完全モデル(full model)といいます。 ただし実際のデータでは説明変数が全く同じものが有り得ます。 例えば表10.3.1には178例のデータがありますが、説明変数の値が全く同じグループが5個あり、これを一般化した表10.3.2には説明変数の値が全く同じグループがm個あります。

このような時は説明変数が(m-1)個あれば線形ロジスティックモデルとデータが完全にフィットします。 このモデルも完全モデルになりますが、本来の完全モデルと区別するために飽和モデルと呼ばれることがあります。 この飽和モデルの尤度と対数尤度は次のようになります。



※pi、ri、niは表10.3.2の表記法に従う

説明変数がp個のモデルと、この飽和モデルの対数尤度の差に尤度比検定を適用することによって異質性(ズレ)の検定を行うことができます。 この検定で用いる検定統計量Dをデビアンス(deviance)といい、モデルとデータのズレの大きさを表す指標になります。 Dは重回帰分析における残差平方和を拡張した統計量に相当します。

対数尤度差:
尤度比検定:D = -2{L(β) - Lf} = χLOF2≧χ2(m-1-p,α)の時、有意水準100α%で有意
D:デビアンス、期待値=自由度=m - 1 - p

異質性の検定としては、これ以外にピアソン残差(Pearson residuals)を利用する手法もあります。 ピアソン残差はモデルから求められる理論的反応例数と実際の反応例数の差を、その標準誤差で標準化した値です。 この値を平方して合計した値は重み付け最小2乗法によるロジスティック回帰分析におけるSLOFに相当するため、これを利用して異質性の検定を行うことができます。

ピアソン残差:   :説明変数がiの時の理論的出現率
≧χ2(m-1-p,α)の時、有意水準100α%で有意

ただしこの手法はniが小さいと不正確になります。 そして最尤法を用いるのはniが小さくて最小2乗法による解の精度が悪い時です。 そのためこの手法はあまり実用的ではありません。 そこで理論的出現率pi(i)の値に基づいてデータを複数のグループに分け(通常は10分位で10個に分割)、各グループの理論的反応例数と実際の反応例数の差を用いて異質性の検定を行う方法が考案されています。 その手法をホスマー・レメショウ(Hosmer-Lemeshow)検定といいます。

≧χ2(g-2,α)の時、有意水準100α%で有意
g:グループ数  Ok:グループkの実際の反応例数   Nk:グループkの例数  πk:グループkの理論的出現率

ただしこの手法はグループ数に恣意性があり、しかもグループ数を変えると結果が変わるため、あまり良い手法とはいえません。 したがって異質性の検定としてはデビアンスを用いる手法が最も実用的ということになります。

ニュートン・ラプソン法は計算を開始する初期値と、反復計算を終了する収束条件によって結果が変わります。 特に初期値は非常に重要です。 初期値が本当の最大値に近いほど反復計算が速く収束し、精度の高い近似解が求められます。 しかし図10.3.2の右側のように極値に近い初期値から計算を始めると、近似解が最大値ではなく極値に収束してしまいます。 また初期値が最大値と大きく離れていると、反復計算を繰り返すたびに近似解がどんどんと大きくなって発散し、近似解が求められないことがあります。

そこで初期値と収束条件を色々と変えて何度も計算し、最大値に最も近いと思われる値を最終的な近似解にするのが理想です。 しかし計算を何度繰り返しても、求めた近似解が本当の解にどの程度近似しているかわからないのが普通です。 そのため実際には反復計算が発散した時以外は何度も計算することはありません。 これはニュートン・ラプソン法の限界であり、最尤法を利用してロジスティック回帰分析を行う方法の欠点です。

ちなみにロジスティック回帰分析の初期値は最小2乗法で求めた値か、それとも判別分析で求めた値を用いると良い近似解が得られる可能性が高くなります。 ただし判別分析で求めた値を初期値にする時は、第1節で説明したように群1と群2の例数から事前確率を計算し、それを用いて定数項をロジスティック回帰分析の定数項のように補正する必要があります。

表10.3.1の例題について実際に計算してみましょう。 最初に計算の準備として、表10.3.1を表10.3.4のような形式に直します。

表10.3.4 ロジスティック回帰分析用テストデータ(最尤法用)
群分類y
1:反応有・0:反応無
説明変数x
110
::
150
010
::
050

初期値として(注1)で最小2乗法によって求めた値を使用します。

初期値:






  

この更新された1を用いて同様の計算を繰り返します。


  


  

23の違いが10-5以下になったので、ここで計算を終了して2を最尤推定値の近似値にします。 この時、-2-1の対角要素を利用して回帰係数の検定と推定を行うことができます。 そして回帰係数を用いてこのモデルの尤度を求めることができるので、尤度比検定を利用して回帰の検定とズレの検定を行うことができます。 また寄与率はロジスティック回帰式によって計算したロジット推定値を説明変数にし、実際のロジットを目的変数にした回帰分析を行い、その時の寄与率から求めます。

ロジスティック回帰式:l=-4.446 + 0.124x
回帰係数の検定:>χ2(1,0.05)=3.841 … 有意水準5%で有意
回帰係数の95%信頼区間=0.124±0.035 → 下限:β1L=0.089 上限:β1U=0.158
オッズ比:OR1=exp(0.124)=1.132
オッズ比の95%信頼区間  下限:OR1L=exp(0.089)=1.093 上限:OR1U=exp(0.158)=1.093
このモデルの対数尤度:L(β1)=-80.4286   AIC=2×(2 + 80.4286)=164.8572
定数項だけのモデルの対数尤度:
飽和モデルの対数尤度:
χβ2=-2×(-120.1130 + 80.4286)=79.369(p=2.5952×10-15)>χ2(1,0.05)=3.841 … 有意水準5%で有意
D=χLOF2=-2×(-80.4286 + 80.0371)=0.783(p=0.8535)<χ2(3,0.05)=7.815 … 有意水準5%で有意ではない
           
(30例)   (35例)
(47例)   (21例)
(45例)
r2=0.982

表10.3.3のように説明変数が全て異なる時、実際のロジットを計算できません。 このような時は対数尤度を利用して擬似寄与率を求める方法があります。 ただしこれはあくまでも疑似寄与率ですから、本来の寄与率とは異なり正確な値ではありません。

擬似寄与率:
このデータの擬似寄与率:

また生存時間解析で用いられるC統計量(C-index、Concordance index)という指標をロジスティック回帰分析に適用することがあります。 この指標は予後予測の精度を表すノンパラメトリックな指標であり、ロジスティック回帰分析に適用するとROC分析におけるROC曲線のAUC(曲線下面積)に相当する値になります。 これは予測値と実際の値の順序が一致している程度を表すノンパラメトリックな指標のため、たとえ1になっても予測値と実際の値がぴったり一致しているとは限りません。 そのため予測精度を表す指標としては寄与率ほど良い指標ではありませんが、一応、紹介しておきます。 (→11.4 比例ハザードモデル (注3))

線形ロジスティックモデルにおいて変数が1つだけで、しかもそれが0または1という値を取るダミー変数とすると、次のように重み付け最小2乗法を用いた時と同様に直線性の検定がオッズ比の検定に一致します。 これは、この場合の最尤法が重み付け最小2乗法に相当するからです。 (→3.4 2標本の計数値 (2)名義尺度(分類データ) (注5))

反応の確率: (x1:群1の時は0、群2の時は1となるダミー変数)
対数尤度:



これをg0の式に代入すると
  
  

  b1=ln(OR)