15.3 出現率のポアソン回帰分析

(1) ダミーデータを用いたポアソン回帰分析

疫学分野では発生件数つまりカウントデータではなく、被験者ごとに疾患の有無を調べたデータにポアソン回帰分析を適用することがあります。 その場合、カウントデータを「0：無　1：有」というダミーデータにし、データの合計例数nを用いてカウントを出現率p = λ/nに変換し、これを目的変数にしてポアソン回帰分析を適用します。 例として第10章の表10.3.1にポアソン回帰分析を適用すると次のようになります。 _(注1)

また第10章で求めたように、同じ表にロジスティックモデル、プロビットモデル、直線モデルを当てはめると次のようになります。 (→10.3 ロジスティック回帰分析の計算方法)

○ロジスティック回帰式(最小2乗法による解)：図15.3.1の青色の曲線
l = -4.378 + 0.122x　　
寄与率：r² = 0.983(98.3%)
直線性(回帰)の検定：χ_β² = 50.789(p = 1.0280×10^-12) ＞ χ²(1,0.05) = 3.841 … 有意水準5％で有意
異質性(ズレ)の検定：χ_LOF² = 0.854(p = 0.8365) ＜ χ²(3,0.05) = 7.815 … 有意水準5％で有意ではない
○変数−プロビット直線回帰式(プロビット曲線)：図15.3.1の赤色の曲線
y(pのプロビット) = 2.445 + 0.071x
寄与率：r² = 0.972(97.2%)
直線性(回帰)の検定：χ_β² = 58.857(p = 1.7347×10^-14) ＞ χ²(1,0.05) = 3.841 … 有意水準5％で有意
異質性(ズレ)の検定：χ_LOF² = 1.680(p = 0.6414) ＜ χ²(3,0.05) = 7.815 … 有意水準5％で有意ではない
○直線回帰式：図15.3.1の緑色の直線
p = -0.262 + 0.0216x
寄与率：r² = 0.945(94.5%)
直線性(回帰)の検定：χ_β² = 68.573(p = 2.3037×10^-15) ＞ χ²(1,0.05) = 3.841 … 有意水準5％で有意
異質性(ズレ)の検定：χ_LOF² = 4.014(p = 0.2599) ＜ χ²(3,0.05) = 7.815 … 有意水準5％で有意ではない

(2) ポアソン回帰分析とロジスティック回帰分析の比較

ポアソン回帰分析の結果とロジスティック回帰分析の結果を比較すると、回帰とズレの検定結果はよく似ているものの、回帰式の切片と回帰係数が少し違っていることがわかります。 そして図15.3.1を見ると出現率が0.2(20％)くらいまでは両者の曲線はあまり乖離しておらず、出現率が0.5より高くなると指数曲線が急激に上昇し、説明変数の値が50くらいで1(100％)を超えてしまうことがわかります。 ポアソン回帰分析は稀に起きる現象が近似的にポアソン分布する性質を利用した分析手法ですから、これは当然のことです。

またロジスティック回帰分析の回帰係数を指数変換するとオッズ比になります。 そしてオッズ比は出現率が低い時(だいたい10％未満)は相対リスクと近似することから、オッズ比を相対リスクと解釈することがあります。 それに対してポアソン回帰分析の回帰係数を指数変換するとリスク比になり、これをそのまま相対リスクと解釈することがあります。 そこで表10.3.1のデータについて説明変数の値が10の時を基準にして、20〜50の相対リスクを実測値、ロジスティック回帰分析の結果、ポアソン回帰分析の結果に基いて計算すると次のようになります。

表15.3.1を見ると、ロジスティック回帰分析の結果に基づいた出現率と相対リスクは実測値に基づいたそれらの値とある程度近似していて、オッズ比は定数倍(3.383倍)になっていることがわかります。 ロジスティック回帰分析は「対数オッズと説明変数の間に線形関係がある」つまり「オッズ比が説明変数と正比例する」というモデルに基づいた手法ですから、これは当然です。

それに対してポアソン回帰分析の結果に基づいた出現率と相対リスクも実測値に基づいたそれらの値とある程度近似しているものの、相対リスクが定数倍(1.8390倍)になっているので説明変数の値が60以上になると出現率が1を超えてしまうことがわかります。 ポアソン回帰分析は「対数出現率と説明変数の間に線形関係がある」つまり「リスク比が説明変数と正比例する」というモデルに基づいた手法ですから、これも当然です。

そして第10章で説明したように、個体が反応する時の説明変数の閾値が正規分布すると説明変数−出現率関数は必然的に累積正規分布(シグモイド曲線)になり、オッズ比と説明変数が近似的に正比例します。 この時、説明変数が一定間隔で増加しても相対リスクは一定の倍率で増加しないので、説明変数の値を指定しないと相対リスクは正確には求められません。 もし相対リスクが説明変数と比例すると表15.3.1のポアソン回帰分析の出現率のように出現率がどこかで必ず1を超えてしまうので、これは当然のことです。 (→10.2 各種のシグモイド曲線)

疫学分野などでは、次のような理由で出現率のデータに対してポアソン回帰分析を適用することがあります。

「ロジスティック回帰分析によって求めたオッズ比は出現率が低い時しか相対リスクに近似しないのに対して、ポアソン回帰分析では相対リスクを正確に求めることができる」

しかし表15.3.1からわかるように、出現率のデータではそもそも説明変数と相対リスクが比例しないので、ポアソン回帰分析よりもロジスティック回帰分析かプロビット回帰分析の方が適しています。 説明変数と相対リスクが近似的に比例するのは、出現率が低くて――だいたい10％未満――ロジスティック曲線と指数曲線が直線で近似できる部分だけです。 そしてその部分ではオッズ比と相対リスクが近似するので、どの手法を用いても結果はあまり変わらないことになります。

ポアソン回帰分析が適しているのは、やはり発生例数はカウントできるものの、発生率を求める時の分母になる全体の例数が不特定多数または非常に膨大で、事実上無限大に近い時のカウントデータでしょう。

(注1)　第2節の(注1)で説明したように、ポアソン回帰分析ではカウントが0のデータも計算に入れることができます。 そのため全例のカウントデータが0/1のダミー変数で表されている時、説明変数の値でグループ化せず、全例を別々のデータとしてポアソン回帰分析を行うとλは理論的出現率を表します。 ただしその場合、尤度比検定による回帰の検定とズレの検定は不正確になります。

そこでそのような場合は説明変数の値が同じグループにおいて、カウントデータが1のケースはカウント数rを求めるために利用し、カウントデータが0のケースはグループの例数nを求めるために利用し、次のようなポアソン回帰モデルを用いて出現率に関するポアソン回帰分析を行うことができます。 このモデルは表10.3.1のように説明変数の値別にカウント数と例数が観測されたデータに対して、出現率に関するポアソン回帰分析を適用する時にも利用できます。

ln(p) = ln(λ/n) = ln(λ) - ln(n) = β₀ + β₁x₁ + … + β_jx_j + … + β_px_p + ε
ln(λ) = ln(n) + β₀ + β₁x₁ + … + β_jx_j + … + β_px_p + ε

p：出現率　　λ：理論的出現例数　　n：例数

このモデルを適用する時はデータを説明変数の値が全て同じケースごとにグループ化し、グループ内のカウントデータの合計数と、グループの合計例数を次のように整理します。

表15.3.2 グループ別カウントデータ

グループ		説明変数
カウント合計	合計例数	x1	…	xj	…	xp
r1	n1	x11	…	x1j	…	x1p
:	:	:		:		:
ri	ni	xi1	…	xij	…	xip
:	:	:		:		:
rm	nm	xm1	…	xmj	…	xmp

グループiの尤度： ← 尤度はカウントデータモデルと同じ
グループiの対数尤度：L(p_i) = L(λ_i) = ln{ℒ(λ_i)} = r_i ln(λ_i) - ln(r_i!) - λ_i
ここで ln(λ_i) = ln(n_i) + ln(p_i) = ln(n_i) + ｘ_i'β、 λ_i = n_i exp(ｘ_i'β) より
グループiの対数尤度：L(β|ｘ_i) = r_i{ln(n_i) + (ｘ_i'β)} - ln(r_i!) - n_i exp(ｘ_i'β)
全体の対数尤度：
　　
　　
w_ik = exp(ｘ_i'b_k) 　　b_k+1 = b_k - Ｈ_k^-1ｇ_k

これらの値を利用して最尤解と尤度を求め、出現率に関するポアソン回帰分析を行うことができます。

このモデルにおいて変数が1つだけで、しかもそれが0または1という値を取るダミー変数とすると、次のようにxの回帰係数は対数リスク差になり、その検定はリスク比の検定に相当します。 ただしこの時の対数リスク差の分散は第3章で説明した分散とは少し異なっています。 これは第3章で説明した対数リスク差の分散がデルタ法によって求めた近似値であるのに対して、こちらは最尤法によって求めた近似値だからです。 (→3.4 2標本の計数値　(2)名義尺度(分類データ) (注5))

(x₁：群1の時は0、群2の時は1となるダミー変数)
　　

∴

∴ 　→　

∴　　b₁=ln(RR)
h₀₀ = -∑{n_i exp(b₀ + b₁)} = -{n₁ exp(b₀) + n₂ exp(b₀ + b₁)} = -(r₁ + r₂)
h₀₁ = h₁₀=-∑{n_i exp(b₀ + b₁)x_i} = -{n₂ exp(b₀ + b₁)} = -r₂
h₁₁ = -∑{n_i exp(b₀ + b₁)x_i²} = h₀₁=-r₂
　　
h₀₀h₁₁ - h₀₁² = (r₁ + r₂)r₂ - (-r₂)² = r₁r₂
　　
※デルタ法によるリスク比の近似分散：

最終更新日：2017年9月4日 　第2節へ 第16章へ

　webmaster@snap-tck.com 　Copyleft (C) 2000 SNAP(Sugimoto Norio Art Production)