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

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

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

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

Ｗ^1/2ｌ=Ｗ^1/2Ｘβ+Ｗ^1/2ε 　　 　　ε_i〜N(0,1²)
　　Ｗ=(Ｗ^1/2)²
正規方程式とその解：[Ｘ'ＷＸ]ｂ=Ｘ'Ｗｌ 　　ｂ=[Ｘ'ＷＸ]^-1Ｘ'Ｗｌ
V(ｂ)=σ²[Ｘ'ＷＸ]^-1≒[Ｘ'ＷＸ]^-1 　　σ²≒1

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

検定：＞χ²(1,α)の時、有意水準100α％で有意
推定：100(1-α)％信頼区間=
→ 下限： 　上限：
[Ｘ'ＷＸ]_jj^-1：[Ｘ'ＷＸ]^-1の第j対角要素

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

　　



　　

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

S_β=b₁²S_xx=χ_β²＞χ²(1,α)の時、有意水準100α％で有意
S_LOF=S_ll-S_β=χ_LOF²＞χ²(m-2,α)の時、有意水準100α％で有意
寄与率：

出現率が0または1の時はロジットが計算できず、重みが0になってしまいます。 そこで出現率が0または1のデータがある時は、それらを計算から除外せずに次のようなロジットを用いて近似解を求める方法があります。

または
または

ただしこの近似計算法は例数が少ない時は誤差が大きくなってしまいます。 そこで出現率が0または1の時のロジットを期待ロジットで補完して繰り返し計算する方法があります。

この方法では、最初は出現率が0または1のデータを除外するか、あるいは上記のロジットを用いて初期解ｂ₀を求めます。 次にそのｂ₀とＸを用いて出現率が0または1のデータのロジット推定値つまり期待ロジットl_iを求め、l_iから期待出現率p_iを逆算して期待重みn_iw_i=n_ip_i(1-p_i)を求めます。 そしてそのl_iとn_iw_iを用いてｂ₁を求め、初期解ｂ₀を更新します。 この計算を繰り返し、ｂが収束したらそれを最終的な解にします。 この繰り返し補完法による解は、長い目で見れば上記のロジットを用いて計算した解よりも精度が高くなります。

表10.3.1の例題について実際に計算してみましょう。 この例題には出現率が0または1のデータがないので補完計算は不必要です。

　　 　　 　　 　　
　　 　　 　　 　　
　　 　　 　　 　　


　　

　　b₀=-0.334-33.181×0.122≒-4.378 　　
b_jの95％信頼区間= → 下限：β_1L=0.088　上限：β_1U=0.155
S_β=χ_β²=0.122²×3419.316≒50.789(p=1.0280×10^-12)＞χ²(1,0.05)=3.841 … 有意水準5％で有意
S_LOF=χ_LOF²=51.643-50.789=0.854(p=0.8365)＜χ²(3,0.05)=7.815 … 有意水準5％で有意ではない

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

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

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

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

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

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

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

f'(x_k)：f(x)の1次微係数　　f''(x_k)：f(x)の2次微係数

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

→

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

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

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

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

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

　　
Ｈ_k=∇ｇ_k 　　
ｂ_k+1=ｂ_k-Ｈ_k^-1ｇ_k

１次偏微分ベクトルｇ_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(b_j)=[Ｉ_f]_jj^-1：Ｉ_f^-1≒[-Ｈ_max]^-1の第j対角要素＝b_jの分散
検定：＞χ²(1,α)の時、有意水準100α％で有意
推定：100(1-α)％信頼区間
→ 下限： 　上限：

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

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

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

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

対数尤度差：
尤度比検定：-2{L₀ - L(β)} = χ_β²＞χ²(p,α)の時、有意水準100α％で有意
AIC = 2{(p+1) - L(β)}

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

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

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

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

対数尤度差：
尤度比検定：D = -2{L(β) - L_f} = χ_LOF²＞χ²(m-1-p,α)の時、有意水準100α％で有意
D：デビアンス、期待値＝自由度＝m - 1 - p

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

ピアソン残差： 　　：説明変数がｘ_iの時の理論的出現率
＞χ²(m-1-p,α)の時、有意水準100α％で有意

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

＞χ²(g-2,α)の時、有意水準100α％で有意
g：グループ数　　O_k：グループkの実際の反応例数 　　N_k：グループkの例数　　π_k：グループkの理論的出現率

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

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

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

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

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

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

初期値：






　　

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

　　


　　

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

ロジスティック回帰式：l=-4.446 + 0.124x
回帰係数の検定：＞χ²(1,0.05)=3.841 … 有意水準5％で有意
回帰係数の95％信頼区間=0.124±0.035 → 下限：β_1L=0.089　上限：β_1U=0.158
オッズ比：OR₁=exp(0.124)=1.132
オッズ比の95％信頼区間 　下限：OR_1L=exp(0.089)=1.093　上限：OR_1U=exp(0.158)=1.171
このモデルの対数尤度：L(β₁)=-80.4286 　　AIC=2×(2 + 80.4286)=164.8572
定数項だけのモデルの対数尤度：
飽和モデルの対数尤度：
χ_β²=-2×(-120.1130 + 80.4286)=79.369(p=2.5952×10^-15)＞χ²(1,0.05)=3.841 … 有意水準5％で有意
D=χ_LOF²=-2×(-80.4286 + 80.0371)=0.783(p=0.8535)＜χ²(3,0.05)=7.815 … 有意水準5％で有意ではない
　　 　　 　　 　　
(30例) 　　(35例)
(47例) 　　(21例)
(45例)
r²=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))

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


∴
これをg₀の式に代入すると
　　
∴ 　　

∴　　b₁=ln(OR)