統計学入門－第10章

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

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

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

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

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

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

表10.1 ロジスティック回帰分析用テストデータ
説明変数x	反応有	反応無	計	出現率p
10	2	28	30	0.067
20	4	31	35	0.114
30	14	33	47	0.298
40	13	8	21	0.619
50	39	6	45	0.867
合計	72	106	178	0.404

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

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

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

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

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

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

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

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

表10.2 説明変数がp個の一般的データ
説明変数					反応有	反応無	計	出現率	ロジット	重み
x₁₁	…	x_1j	…	x_1p	r₁	n₁-r₁	n₁	p₁	l₁	n₁w₁
:		:		:	:	:	:	:	:	:
x_i1	…	x_ij	…	x_ip	r_i	n_i-r_i	n_i	p_i	l_i	n_iw_i
:		:		:	:	:	:	:	:	:
x_m1	…	x_mj	…	x_mp	r_m	n_m-r_m	n_m	p_m	l_m	n_mw_m

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

・線形ロジスティックモデル

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

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

・正規方程式とその解

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

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

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

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

…有意水準5％で有意ではない

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

表10.3 説明変数がp個の一般的データ(最尤法用)
群分類y 1:反応有・0:反応無	説明変数
群分類y 1:反応有・0:反応無	x₁	…	x_j	…	x_p
1	x₁₁	…	x_1j	…	x_1p
:	:		:		:
1	x_r1	…	x_rj	…	x_rp
0	x_(r+1)1	…	x_(r+1)j	…	x_(r+1)p
:	:		:		:
0	x_n1	…	x_nj	…	x_np

・線形ロジスティックモデル

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

反応有の個体を得る確率：

反応無の個体を得る確率：

両者合わせて表すと

全体の尤度：

対数尤度：

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

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

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

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

この解x_k+1のところでまたテーラー展開を行い、次の解を求めれば、次第に最大値に近づくことになります。このような反復計算を繰り返し、x_kとx_k+1の差が十分に小さくなったところで、x_kを最大値の最終的な近似値にします。

関数の変数が複数の場合は、行列とベクトルを利用して次のように計算します。

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

最尤推定値の漸近的正規性を利用して、偏回帰係数の検定を行うことができます。これを「ワルド(Wald)の検定」といいます。

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

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

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

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

対数尤度差：

尤度比検定：

の時有意水準100α％で有意
AIC=2{(p+1)-L(β)}

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

このようなモデルも完全モデルですが、本来の完全モデルと区別するために飽和モデルと呼ばれることがあります。この飽和モデルの尤度と対数尤度は次のようになります。

※p_i、r_i、n_iは表10.2の表記法に従う

説明変数がp個のモデルとこの飽和モデルの対数尤度の差に尤度比検定を適用することによって、異質性(ズレ)の検定を行うことができます。

対数尤度差：

尤度比検定：

の時有意水準100α％で有意
D：デビアンス(deviance)…モデルとデータのズレの大きさを表す指標、期待値＝自由度(m-1-p)

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

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

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

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