前口上 | 目次 | 第1章 | 第2章 | 第3章 | 第4章 | 第5章 | 第6章 | 第7章 | 第8章 | 第9章 | 第10章 |
第11章 | 第12章 | 第13章 | 第14章 | 第15章 | 第16章 | 第17章 | 第18章 | 第19章 | 第20章 | 付録 |
1 | 2 | 3 | 4 | 5 |
ロジスティック回帰分析ではロジットの回帰誤差が特殊な分布になり、普通はその分布を理論的に確定することができません。 そのため回帰誤差が近似的に正規分布すると見なして重回帰分析と同じように最小2乗法を利用して回帰分析を行う方法と、最尤法を利用した繰り返し近似計算によって回帰分析を行う方法の2種類があります。
最初に最小2乗法を利用する方法について説明しましょう。 この場合の線形ロジスティックモデルは次のようになります。 このモデルでは説明変数xjの値を研究者が任意に指定した時、ロジットの回帰誤差εが近似的に正規分布すると仮定します。 ただしロジスティック回帰式を計算するには回帰誤差の正規性は必要ではなく、検定を行う時だけ回帰誤差の正規性が必要になります。 そして説明変数は研究者が任意の値を指定するので誤差がなく、正規性という概念そのものが当てはまりません。 重回帰分析と同様に、これはよく誤解されていることです。
このモデルは目的変数が少し特殊なだけで、重回帰モデルとほとんど同じです。 そのため計算方法も重回帰分析とほとんど同じです。 例として表10.3.1のデータにこのモデルを当てはめ、最小2乗法を利用して解を求めると次のようになります。 (注1)
説明変数x | 反応無 | 反応有 | 計 | 出現率p |
---|---|---|---|---|
10 | 28 | 2 | 30 | 0.067 |
20 | 31 | 4 | 35 | 0.114 |
30 | 33 | 14 | 47 | 0.298 |
40 | 8 | 13 | 21 | 0.619 |
50 | 6 | 39 | 45 | 0.867 |
合計 | 106 | 72 | 178 | 0.404 |
重回帰分析と同様にロジスティック回帰分析も記述統計学的手法なので、推測統計学的手法である検定とは相性が良くありません。 それでも有意症患者のために回帰誤差が近似的に正規分布するという性質を利用して検定を行うことができます。
直線性の検定は出現率の変動の中にロジスティック曲線で説明できる変動があるかどうか、つまり寄与率が0かどうかの検定です。 異質性(ズレ)の検定は出現率の変動の中にロジスティック曲線では説明できない変動があるかどうか、つまり実際の出現率とロジスティック曲線のズレが0かどうかの検定です。 ただしこれらの検定で検出差を指定するのは困難なので、たいていの場合は有意性検定になります。 そのため実質的な意味はほとんどないと考えた方が良いと思います。
ちなみに表10.3.1のデータにコクラン・アーミテージの傾向分析を適用し、説明変数xと出現率pを直線で回帰すると次のようになります。 そしてこの回帰直線と最小2乗法を利用して計算したロジスティック回帰曲線をグラフ化すると、図10.3.1の黒色の直線と青色の曲線になります。
また表10.3.1のデータをプロビット変換してから回帰分析を行うと次のようになります。 そしてこの変数−プロビット直線を変数−出現率曲線に変換すると、図10.3.1の赤色の曲線になります。 表10.3.1の反応有の閾値分布が近似的に正規分布するなら、ロジスティック回帰曲線よりもこのプロビット回帰曲線の方が正確です。
最尤法を利用する方法は回帰誤差が近似的に正規分布しなくても適用できます。 その代わり厳密な解を代数的に求めることができないので、適当な初期値を用いた繰り返し計算によって近似解を求めます。 そのため初期値と繰り返し計算の収束条件を変えると計算結果が変わってしまいます。 したがって使用する統計ソフトや計算条件によって結果が変わることがあり、最小2乗法を利用する方法と違って誰が計算しても同じ結果になるというわけではありません。
判別分析も最尤法によって判別関数を求めます。 しかし判別分析の場合は2群を研究者が任意に指定し、それらの群の検査項目の値が近似的に正規分布すると仮定して最尤法を適用します。 つまり目的変数である群の出現率には誤差がなく、説明変数である検査項目に誤差があると仮定して最尤法を適用するわけです。 それに対してこの場合の最尤法は説明変数には誤差がなく、目的変数であるロジットまたは出現率に誤差があると仮定して最尤法を適用します。
表10.3.1のデータについて最尤法を利用する方法でロジスティック回帰分析を行うと次のようになります。 ここでは最小2乗法を利用して計算した解を初期値にしたので、結果はほとんど変わりません。 したがってロジスティック回帰曲線もほとんど変わらず、図10.3.1の青色の曲線とほとんど重なってしまいます。 (注2)
この計算法では回帰誤差の分布を確定しないので、本来は検定と推定はできません。 しかし何が何でも検定しなければ気がすまない有意症患者のために、最尤法による解が漸近的に正規分布するという性質を利用して上記のように回帰係数の検定と推定、そして直線性(回帰)の検定と異質性(ズレ)の検定を行うことがあります。
最尤法を利用する方法は回帰誤差が近似的に正規分布しなくても適用できるところに特徴があります。 ところが上記の検定と推定は最尤法による解が近似的に正規分布する性質を利用しています。 最尤法による解が厳密に正規分布するのは回帰誤差が正規分布する時であり、その時は最小2乗法による解と最尤法による解が一致し、各種の検定結果も一致します。 したがってどうせ近似的な正規性を利用するなら、回帰誤差が近似的に正規分布すると仮定して計算する方が便利です。
また最尤法を利用する方法は適当な初期値を用いた繰り返し計算によって近似解を求め、近似解が収束したところで計算を終了する近似計算法です。 そのためデータが極端に偏っていたり誤差が多かったりすると繰り返し計算が発散して異常な解になったり、解が求められなかったりすることがあります。 例えば表9.1.1のデータを前向き研究で得られたデータと見なしてロジスティック回帰分析を適用すると、次のように判別分析から得られるロジスティック曲線とはかなり異なった結果になります。
この原因は表9.1.1のデータがきれいすぎる、つまり極端に偏っているからです。 図9.6.1を見るとわかるように、表9.1.1のデータに判別分析を適用すると25例のデータが全て正しく判別されてしまいます。 そのためロジスティック回帰分析では繰り返し計算が発散して異常な解になってしまったのです。 これらのことからわかるように、最尤法を利用した計算法よりも最小2乗法を利用した計算法の方が色々な意味で便利です。
しかし説明変数が多数になると、全ての説明変数の値が同じというデータが複数あるのは稀です。 そうするとデータを表10.3.1のような形式にまとめた場合、出現率が0または1というテータが大多数になります。 出現率が0または1の時、ロジット変換はできません。 そのためそのようなデータがあると、最小2乗法を利用した計算法では期待ロジットを用いた繰り返し補完法によって近似解を求めます。 しかし繰り返し補完法は出現率が0または1のデータが少数の時の計算法であり、そのようなデータが大多数の時は結果の精度が低くなってしまいます。
それに対して最尤法を利用した方法ではロジット変換を直接利用せずに特殊な方法で計算を行うので、出現率が0または1のデータも計算に入れることができます。 そのためコンピュータが発達したことと相まって、現在は最尤法を利用する方法が主流になっています。
説明変数 | 反応有 | 反応無 | 計 | 出現率 | ロジット | 重み | ||||
---|---|---|---|---|---|---|---|---|---|---|
x11 | … | x1j | … | x1p | r1 | n1-r1 | n1 | p1 | l1 | n1w1 |
: | : | : | : | : | : | : | : | : | ||
xi1 | … | xij | … | xip | ri | ni-ri | ni | pi | li | niwi |
: | : | : | : | : | : | : | : | : | ||
xm1 | … | xmj | … | xmp | rm | nm-rm | nm | pm | lm | nmwm |
重みniwiは重み付け最小2乗法に用いるもので、ロジットの分散の逆数です。 この重みは誤差の少ないロジットほど重要視するためのものであり、コクラン・アーミテージの傾向検定の重みと同じように出現率が二項分布することから導き出されます。 (→5.3 計数値の相関分析と回帰分析 (注4)、7.1 重回帰モデル (注1))
線形ロジスティックモデルの両辺に重みの平方根√(niwi)をかけ、重み付け最小2乗法によってβの最良不偏推定値(BLUE解)bを求めると次のようになります。
回帰誤差εiが近似的に正規分布することを利用して、偏回帰係数の検定と推定を行うことができます。
説明変数が1つの時は次のようになります。
出現率が0または1の時はロジットが計算できず、重みが0になってしまいます。 そこで出現率が0または1のデータがある時は、それらを計算から除外せずに次のような実用ロジットを用いて近似解を求める方法があります。
ただしこの近似計算法は例数が少ない時は誤差が大きくなってしまいます。 そこで出現率が0または1の時のロジットを期待ロジットで補完して繰り返し計算する方法もあります。
この方法では、最初は出現率が0または1のデータを除外するか、あるいは上記の実用ロジットを用いて初期解b0を求めます。 次にそのb0とXを用いて出現率が0または1のデータのロジット推定値つまり期待ロジットliを求め、liから期待出現率piを逆算して期待重みniwi = nipi(1 - pi)を求めます。 そしてそのliとniwiを用いてb1を求め、初期解b0を更新します。 この計算を繰り返して、bが収束したらそれを最終的な解にします。 この繰り返し補完法による解は、長い目で見れば上記の実用ロジットを用いて計算した解よりも精度が高くなります。
表10.3.1の例題について実際に計算してみましょう。 この例題には出現率が0または1のデータがないので補完計算は不必要です。
m = 2でx1 = 0、x2 = 1とすると、次のように直線性の検定がオッズ比の検定に一致します。 これは直線性の検定はロジットの回帰係数つまり対数オッズの差の検定に相当し、回帰係数を指数変換するとオッズ比になるためです。 (→3.4 2標本の計数値 (2)名義尺度(分類データ) (注5))
群分類y 1:反応有・0:反応無 | 説明変数 | ||||
---|---|---|---|---|---|
x1 | … | xj | … | xp | |
1 | x11 | … | x1j | … | x1p |
: | : | : | : | ||
1 | xr1 | … | xrj | … | xrp |
0 | x(r+1)1 | … | x(r+1)j | … | x(r+1)p |
: | : | : | : | ||
0 | xn1 | … | xnj | … | xnp |
回帰誤差εiが正規分布すると仮定せず、最尤法によってβの最尤推定値bを求めます。
対数尤度を最大にするbをニュートン・ラプソン(Newton-Raphson)法によって近似的に求めます。 ニュートン・ラプソン法はニュートンの2階勾配法とも呼ばれ、関数を2次式で近似し、関数の最大値または最小値を反復計算によって近似的に求める方法です。
今、最大値(または最小値)を持つ関数をf(x)とします。 一般にf(x)は最大値の近傍では2次関数つまり放物線で近似できます。 そこでf(x)を最大値に近い値xkでテーラー展開(Tayler expansion)し、2次の項まで取ります。
この近似関数を最大にする時のxは、関数をxで微分して0と置いた方程式の解になります。
f(x)の1次微係数と2次微係数が直接求められない時は、次のような差分法によって近似的に求めます。
この解xk+1のところでまたテーラー展開を行い、次の解を求めれば次第に最大値に近づくはずです。 このような反復計算を繰り返し、xkとxk+1の差が十分に小さくなったところでxkを最終的な近似解にします。 関数の変数が複数の時は行列とベクトルを利用して次のように計算します。
このニュートン・ラプソン法を対数尤度関数に適用すると次のようになります。
1次偏微分ベクトルgkはL(bk)の傾斜を表す傾斜ベクトルであり、最大値付近ではgmax ≒ 0(ゼロベクトル)になります。 そして2次偏微分係数行列Hk(ヘスの行列)はL(bk)の曲率を表す計量行列です。 最大値付近におけるこの曲率が大きいほどbmaxの変動幅が小さくなるので、これはbmaxの確実性を表す指標になります。
ただしL(bk)は上に凸の曲線なので、この曲率は負の値になります。 そこでこの曲率の符号を反転させた値の期待値を求めると、それは推定量の確実性を表す指標つまりデータXが持つ情報量の多さを表す指標になります。 それをフィッシャーの情報量(Fisher's information)と呼び、Hmaxの符号を反転させた行列の期待値E(-Hmax)をIfと書いてフィッシャーの情報行列(Fisher's information matix)と呼んでいます。
最尤推定値には漸近的正規性があり、フィッシャーの情報量の逆数が最尤推定量の分散になります。 そこでこれらの性質を利用して偏回帰係数の検定と推定を行うことができます。 この検定をワルド(Wald)の検定といいます。 (→1.4 推定 (注4))
説明変数がなく、定数項だけの線形ロジスティックモデルの尤度と対数尤度は次のようになります。
尤度は「もっともらしさ」ですから、この最も単純なモデルの尤度と説明変数がp個の時の尤度の比が大きいほど説明変数の寄与が少なく、尤度比が小さいほど寄与が大きいことになります。 (→9.3 1変量の場合 (1) 尤度と最尤法)
尤度比の逆数を平方して対数変換した値、つまり対数尤度差に(-2)を掛けた値が近似的にχ2分布することを利用した検定法を尤度比検定(likelihood ratio test)といいます。 この尤度比検定を定数項だけのモデルと説明変数がp個のモデルの対数尤度の差に適用すれば、説明変数全体の回帰の検定つまり直線性の検定を行うことができます。 またモデルの当てはまりの良し悪しを表す指標として、赤池の情報量基準AIC(Akaike's Information Criterion)を計算することもできます。 (→(2) 名義尺度(分類データ) 1) データに対応がない場合、9.3 1変量の場合 (注1))
一方、説明変数と定数項を合わせた数(p + 1)が例数nと同じ時、線形ロジスティックモデルとデータが完全にフィットします。 これを完全モデル(full model)といいます。 ただし実際のデータでは説明変数が全く同じものが有り得ます。 例えば表10.3.1には178例のデータがありますが、説明変数の値が全く同じグループが5個あり、これを一般化した表10.3.2には説明変数の値が全く同じグループがm個あります。
このような時は説明変数が(m - 1)個あれば線形ロジスティックモデルとデータが完全にフィットします。 このモデルも完全モデルになりますが、本来の完全モデルと区別するために飽和モデルと呼ばれることがあります。 この飽和モデルの尤度と対数尤度は次のようになります。
説明変数がp個のモデルと、この飽和モデルの対数尤度の差に尤度比検定を適用することによって異質性(ズレ)の検定を行うことができます。 この検定で用いる検定統計量Dをデビアンス(deviance)といい、モデルとデータのズレの大きさを表す指標になります。 Dは重回帰分析における残差平方和を拡張した統計量に相当します。
異質性の検定としては、これ以外にピアソン残差(Pearson residuals)を利用する手法もあります。 ピアソン残差はモデルから求められる理論的反応例数と実際の反応例数の差を、その標準誤差で標準化した値です。 この値を平方して合計した値は重み付け最小2乗法によるロジスティック回帰分析におけるSLOFに相当するので、これを利用して異質性の検定を行うことができます。
ただしこの手法はniが小さいと不正確になります。 そして最尤法を用いるのはniが小さくて最小2乗法による解の精度が悪い時です。 そのためこの手法はあまり実用的ではありません。 そこで理論的出現率pi(xi)の値に基づいてデータを複数のグループに分け(通常は10分位で10個に分割)、各グループの理論的反応例数と実際の反応例数の差を用いて異質性の検定を行う方法が考案されています。 その手法をホスマー・レメショウ(Hosmer-Lemeshow)検定といいます。
ただしこの手法はグループ数に恣意性があり、しかもグループ数を変えると結果が変わるのであまり良い手法とはいえません。 したがって異質性の検定としてはデビアンスを用いる手法が最も実用的ということになります。
ニュートン・ラプソン法は計算を開始する初期値と、反復計算を終了する収束条件によって結果が変わります。 特に初期値は非常に重要です。 初期値が本当の最大値に近いほど反復計算が速く収束し、精度の高い近似解が求められます。 しかし図10.3.2の右側のように極値に近い初期値から計算を始めると、近似解が最大値ではなく極値に収束してしまいます。 また初期値が最大値と大きく離れていると、反復計算を繰り返すたびに近似解がどんどんと大きくなって発散し、近似解が求められないことがあります。
そこで初期値と収束条件を色々と変えて何度も計算し、最大値に最も近いと思われる値を最終的な近似解にするのが理想です。 しかし計算を何度繰り返しても、求めた近似解が本当の解にどの程度近似しているかわからないのが普通です。 そのため実際には反復計算が発散した時以外は何度も計算することはありません。 これはニュートン・ラプソン法の限界であり、最尤法を利用してロジスティック回帰分析を行う方法の欠点です。
ちなみにロジスティック回帰分析の初期値は最小2乗法で求めた値か、それとも判別分析で求めた値を用いると良い近似解が得られる可能性が高くなります。 ただし判別分析で求めた値を初期値にする時は、第1節で説明したように群1と群2の例数から事前確率を計算し、それを用いて切片をロジスティック回帰分析の切片のように補正する必要があります。
表10.3.1の例題について実際に計算してみましょう。 最初に計算の準備として、表10.3.1を表10.3.4のような形式に直します。
群分類y 1:反応有・0:反応無 | 説明変数x |
---|---|
1 | 10 |
: | : |
1 | 50 |
0 | 10 |
: | : |
0 | 50 |
初期値として(注1)で最小2乗法によって求めた値を使用します。
この更新されたb1を用いて同様の計算を繰り返します。
b2とb3の違いが10-5以下になったので、ここで計算を終了してb2を最尤推定値の近似値にします。 この時、-H2-1の対角要素を利用して回帰係数の検定と推定を行うことができます。 そして回帰係数を用いてこのモデルの尤度を求めることができるので、尤度比検定を利用して回帰の検定とズレの検定を行うことができます。 また寄与率はロジスティック回帰式によって計算したロジット推定値を説明変数にし、実際のロジットを目的変数にした回帰分析を行い、その時の寄与率から求めます。
表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))