(注1)　表15.1.1を一般化し、ポアソン回帰モデルを当てはめると次のようになります。

○ポアソン回帰モデル
　　
　　
　　 　　 　　

このモデルでは回帰誤差ε_iが正規分布しないので、最尤法によってβの最尤推定値ｂを求めます。 そのための準備として、まず表15.2.1のデータを説明変数の値が全て同じケースごとにグループ化します。 そのグループがm個あり、i番目のグループの例数をn_i、グループ内のカウントデータの合計数をr_iとすると、それらは表15.2.2のように整理することができます。

本来のカウントデータは非常に稀にしか起こらない事象が起こった時に、それに関する説明変数の値を観測したものなので、yの合計が0つまりカウントが0のものはありません。 しかしポアソン分布はカウントが0の時も含んでいるのでカウントが0のグループが存在しても適用可能です。 r_iがポアソン分布すると仮定すると次のような式が成り立ちます。

カウントデータがr_iになる確率：
※カウントデータが0の時の確率：
グループiの尤度：
グループiの対数尤度：L(λ_i) = ln{ℒ(λ_i)} = r_i ln(λ_i) - ln(r_i!) - λ_i
ここで ln(λ_i)=ｘ_i'β、λ_i = exp(ｘ_i'β) より
グループiの対数尤度：L(β|ｘ_i) = r_i(ｘ_i'β) - ln(r_i!) - exp(ｘ_i'β)
全体の対数尤度：

この全体の対数尤度関数にニュートン・ラプソン法を適用し、最尤解を求めると次のようになります。 (→10.3 ロジスティック回帰分析の計算方法　(2) 最尤法を利用する方法 (注2))

　　
　　
w_ik = exp(ｘ_i'b_k) 　　b_k+1 = b_k - Ｈ_k^-1ｇ_k

ロジスティック回帰分析と同様に、ワルドのχ²値によって偏回帰係数が0かどうかの検定と推定を行うことができます。

V(ｂ) = -Ｈ_k^-1 ≒ Ｉ_f^-1 　　E(-Ｈ_k)=Ｉ_f：情報行列 　　[Ｉ_f]_jj^-1：Ｉ_f^-1の第j対角要素
検定： ＞ χ²(1,α)の時、有意水準100α％で有意
推定：100(1-α)％信頼区間
→ 下限： 　上限：

また偏回帰係数が全て0の時の尤度つまり説明変数が無くて定数項だけのモデルの尤度と、定数項と説明変数がp個のモデルの尤度の比を利用した尤度比検定によって、説明変数全体の回帰の検定を行うことができます。 さらに飽和モデルの尤度を利用してモデルとデータのズレつまり異質性の検定を行うことができます。 定数項だけのモデルの最尤推定値λ₀と尤度、そして飽和モデルの尤度は次のようにして求めることができます。

定数項だけのモデルの対数尤度： 　　尤度の自由度 = 1

∴ 　　
対数尤度差： 　　尤度の自由度 = (p + 1) - 1 = p
回帰の尤度比検定：-2{L(0) - L(β)} = χ_β² ＞ χ²(p,α)の時、有意水準100α％で有意
飽和モデルの対数尤度： 　　尤度の自由度 = m
　　∴λ_i = r_i
対数尤度差： 　　尤度の自由度=m - (p + 1) = m - p - 1
ズレの尤度比検定：D = -2{L(β)) - L_f}=χ_LOF²＞χ²(m-p-1,α)の時、有意水準100α％で有意
D：デビアンス(deviance) … モデルとデータのズレの大きさを表す指標

偏回帰係数の初期値はカウントデータyが小さい時は指数関数を直線によって近似できることを利用して求めます。 つまり表15.2.2のデータに重回帰分析を適用して偏回帰係数を求め、それを初期値として用いるわけです。 切片については、全ての説明変数に平均値を代入した時のλが定数項だけのモデルの最尤推定値λ₀と一致するように調整します。

表15.1.1のデータについて実際に計算してみましょう。 この表の16番と33番は説明変数の値が同じですから、この2つを1つのグループにして発生件数は5として計算します。

○yを目的変数にした重回帰分析の結果
y = 1.65901 - 0.34666x₁ + 0.315884x₂ - 0.00111275x₃
… カウント数と例数を用いた最尤推定値の近似値
ln(2.04545) = b₀₀ - 0.34666×0.545455 + 0.315884×1.95455 - 0.00111275×37.6364
∴b₀₀ = 0.715618 - 0.386444 = 0.329174
初期値：
　　
　　

更新されたｂ₁を用いて同様の計算を繰り返すと、4回目で値が収束します。

　　
　　
ｂ₅ = ｂ₄ - Ｈ₄^-1ｇ₄ ≒ ｂ₄
○ポアソン回帰式：l = ln(λ) = 0.619135 - 0.227517x₁ + 0.157503x₂ - 0.00195712x₃
○偏回帰係数の検定と推定

95％信頼区間　下限：β_1L = -0.651551　上限：β_1U = 0.196516

95％信頼区間　下限：β_2L = -0.100592　上限：β_2U = 0.415598

95％信頼区間　下限：β_3L = -0.0215049　上限：β_3U = 0.0175906
○尤度比検定
このモデルの対数尤度：L(β) = -62.2636
定数項だけのモデルの対数尤度：L(０) = -63.389 　　飽和モデルの対数尤度：L_f = -55.4038
回帰：χ_β² = -2×(-63.389 + 62.2636) = 2.251(p = 0.5220) ＜ χ²(3,0.05) = 7.815
ズレ：χ_LOF² = -2×(-62.2636 - 55.4038) = 13.720(p = 0.9999) ＜ χ²(39,0.05) = 54.572
擬似寄与率：

ポアソン回帰分析でもワルドの検定に用いるχ²値を変数選択用統計量にして、ロジスティック回帰分析と同様の手順で変数選択を行うことが原理的には可能です。 しかし説明変数の組み合わせが変わると説明変数の値が全て同じケースの組み合わせが変わってしまい、最尤法の計算が非常に煩雑になります。 そのため普通はロジスティック回帰分析のような変数選択は行いません。

SASやRといった既存の統計ソフトでは、一般化線形モデル用関数(GLM等)を利用してポアソン回帰分析を行うことができます。 その場合、他の一般化線形モデルと同様に1つのデータを1つのカウントデータとして最尤法を適用します。 そのため表15.1.1のように説明変数の値が同じカウントデータが複数あっても、これをグループ化せずに計算するので結果が不正確になってしまいます。 そこで表15.1.1のデータを正しく解析するには、次のような表にしてから解析する必要があります。