玄関雑学の部屋雑学コーナー統計学入門

15.2 ポアソン回帰分析結果の解釈

(1) ポアソン回帰分析の計算結果

ポアソン回帰モデルでは理論的発生件数λがポアソン分布すると仮定するため、ln(λ)の回帰誤差εはポアソン分布を対数変換した分布をし、正規分布しません。 そこでロジスティック回帰分析と同様に、最小2乗法ではなく最尤法を利用して回帰分析を行います。 そして寄与率の代わりに尤度を利用した擬似寄与率を求め、各種の検定を行うことができます。 表15.1.1のデータにポアソン回帰モデルを当てはめ、最尤法を利用して解を求めると次のようになります。 (注1)

○ポアソン回帰モデル
ln(λ)=β0 + β1x1 + β2x2 + β3x3 + ε
λ:理論的発生件数  x1:診療科(0:内科系 1:外科系)   x2:処方薬剤数  x3:診療科職員数  ε:回帰誤差
○ポアソン回帰式
l=ln(λ)=0.619135 - 0.227517x1 + 0.157503x2 - 0.00195712x3
λ=exp(l)=exp(0.619135 - 0.227517x1 + 0.157503x2 - 0.00195712x3)
擬似寄与率:R2=0.141(14.1%)
○偏回帰係数の検定と推定
χβ12=1.106(p=0.2930)<χ2(1,0.05)=3.841 … 有意水準5%で有意ではない
x 95%信頼区間 下限:β1L=-0.651551 上限:β1U=0.196516
χβ22=1.431(p=0.2317)<χ2(1,0.05)=3.841 … 有意水準5%で有意ではない
95%信頼区間 下限:β2L=-0.100592 上限:β2U=0.415598
χβ32=0.039(p=0.8444)<χ2(1,0.05)=3.841 … 有意水準5%で有意ではない
95%信頼区間 下限:β3L=-0.0215049 上限:β3U=0.0175906
○尤度比検定
回帰の検定:χβ2=2.251(p=0.5220)<χ2(3,0.05)=7.815 … 有意水準5%で有意ではない
ズレの検定:χLOF2=13.720(p=0.9999)<χ2(39,0.05)=54.572 … 有意水準5%で有意ではない
図15.2.1 ポアソン回帰による指数曲線

図15.2.1は横軸をポアソン回帰式から求めたl=ln(λ)の値にし、縦軸を実際の発生件数にして表15.1.1の各データをプロットし、さらに指数曲線λ=exp(l)を描いたものです。 このグラフを見ると、実際の発生件数のデータと指数曲線があまりフィットしていないことがわかります。 このことは解析結果の擬似寄与率が14.1%と非常に小さく、回帰の検定結果が有意ではないことからもわかります。

実は表15.1.1のデータは、順序ロジスティック回帰分析の例題として用いた表10.5.1のデータを項目名を変えてそのまま流用したものです。 そのため発生件数のデータがポアソン分布せず、ポアソン回帰モデルがあまり当てはまりません。 したがってポアソン回帰分析の例題としてはあまり良いデータではありませんが、重回帰分析と比較するために敢えてこのデータを用いました。

このデータに重回帰分析を適用すると、第10章で求めたように次のような結果になります。 (→10.5 順序ロジスティック回帰分析)

○重回帰式
y=1.659 - 0.347x1 + 0.316x2 - 0.00111x3
重寄与率:R2=0.140(14%)  重相関係数:R=0.374
○偏回帰係数の検定と推定
Fβ1=2.288(p=0.1382)<F(1,41,0.05)=4.079 … 有意水準5%で有意ではない
95%信頼区間 下限:β1L=-0.810 上限:β1U=0.116
Fβ2=4.926(p=0.0322)>F(1,41,0.05)=4.079 … 有意水準5%で有意
95%信頼区間 下限:β2L=0.028 上限:β2U=0.603
Fβ3=0.0115(p=0.9153)<F(1,41,0.05)=4.079 … 有意水準5%で有意ではない
95%信頼区間 下限:β3L=-0.022 上限:β3U=0.020

モデルは異なりますが、偏回帰係数の符号と値がだいたい似ていることと、検定結果と重寄与率もだいたい似ていることがわかります。 このことから、このデータには計算が簡単で結果の解釈も容易な重回帰分析を適用した方が実用的であることがわかります。 そもそもこのデータは順序ロジスティック回帰分析の説明用ですから、これは致し方ありません。

(2) 各種パラメーターの意味

ポアソン回帰式の偏回帰係数はx他の説明変数が一定で注目している説明変数だけが「1」増加した時にln(λ)がいくつ変化するかを表す値、つまり理論的発生件数λを対数変換した値の変化量を表す値です。 そのため偏回帰係数を指数変換するとλの比、つまりλが相対的に何倍になるかを表す値になります。 例えば(1)の解析結果で各変数の偏回帰係数を指数変換すると次のようになります。

x1(診療科):exp(-0.227517)=0.796509 … 内科系に比べて外科系の医療事故件数は約0.8倍ある
x2(処方薬剤数):exp(0.157503)=1.17058 … 処方薬剤数が1個増加すると医療事故件数は約1.2倍になる
x3(診療科職員数):exp(-0.00195712)=0.998045 … 職員数が1名増加すると医療事故件数は約0.998倍になる

これらの値は医療事故の発生件数の比であり、分母が共通ですから発生率の比つまりリスク比と解釈することが可能です。 そしてリスク比は近似的に相対リスク(相対危険度)と解釈することができるので、これらの値を調整相対リスクと解釈することができます。

また重回帰分析と同様にポアソン回帰分析も記述統計学的手法であり、推測統計学的手法である検定とは相性が良くありません。 そのため偏回帰係数の検定は、たいてい単なる有意性検定になります。 したがって検定結果よりもポアソン回帰式や疑似寄与率を科学的に検討する方が有意義です。 しかし検定に用いるχ2値や有意確率p値を各変数の相対的な重要度の指標として利用することはできます。 (1)の解析結果中のそれらの値を見ると、あまり大きな違いはないものの、一応、x2(処方薬剤数)が一番大きく影響していると解釈することができそうです。


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

表15.2.1 一般的カウントデータ
目的変数
カウントデータ
説明変数
x1xjxp
y1x11x1jx1p
::::
yixi1xijxip
::::
ynxn1xnjxnp
○ポアソン回帰モデル
  
  
        

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

表15.2.2 グループ別
カウントデータ
グループ説明変数
yの合計x1xjxp
r1x11x1jx1p
::::
rixi1xijxip
::::
rmxm1xmjxmp

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

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

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

  
  
wik=exp(i'bk)   bk+1=bk-k-1k

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

V()=-k-1f-1   E(-k)=f:情報行列   [f]jj-1f-1の第j対角要素
検定:>χ2(1,α)の時、有意水準100α%で有意
推定:100(1-α)%信頼区間
→ 下限:  上限:

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

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

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

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

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

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

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

  
  
5=4 - 4-144
○ポアソン回帰式:l=ln(λ)=0.619135 - 0.227517x1 + 0.157503x2 - 0.00195712x3
○偏回帰係数の検定と推定

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   飽和モデルの対数尤度:Lf=-55.4038
回帰:χβ2=-2×(-63.389 + 62.2636)=2.251(p=0.5220)<χ2(3,0.05)=7.815
ズレ:χLOF2=-2×(-62.2636 - 55.4038)=13.720(p=0.9999)<χ2(39,0.05)=54.572
擬似寄与率:

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

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

表15.2.3 医療事故の発生件数(説明変数の値でグループ化)
医療機関ID発生件数診療科 (0:内科系 1:外科系)処方薬剤数診療科職員数
110121
210130
310137
410246
511124
611156
711158
811224
911238
1011258
1111326
1211341
1320123
1420143
1520147
1720241
1820245
1920253
2020340
2121122
2221139
2321152
2421223
2521228
2621232
2721243
2821324
2921327
3021342
3130120
3230144
3350235
3430237
3530341
3630355
3730351
3830336
3931134
4031242
4131251
4231321
4331335
4431336