6. ポアソン回帰：カウントデータの統計モデリング

カウントデータとポアソン分布

カウントデータは、ある単位時間や単位面積内で発生する事象の件数を表す非負整数値のデータです。交差点での1時間あたりの交通事故件数、ウェブページの日次アクセス数、製品1ロットあたりの欠陥個数、ある地域での月間感染症罹患数などが代表的な例です。カウントデータは値が0以上の整数に限定されており、期待値が小さい条件下では分布が右に歪む傾向を持ちます。こうした分布の特性は正規分布の仮定に適合しないため、通常の最小二乗法による線形回帰モデルの直接適用は適切ではありません。負の予測値が生じ得ることや、正規誤差の仮定が成立しないことが根本的な問題となります。

ポアソン分布はポアソン過程から導出される離散確率分布です。ポアソン過程とは、事象が互いに独立に、かつ単位区間内で一定の平均発生率のもとで生じるランダム過程を指します。カウント $Y$ がパラメータ $\lambda > 0$ のポアソン分布に従うとき、その確率質量関数は

$$P(Y = k;\, \lambda) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0, 1, 2, \ldots$$

で定義されます。パラメータ $\lambda$ はこの分布の期待値かつ分散を同時に表します。すなわち $E[Y] = \text{Var}[Y] = \lambda$ という等分散性が成立します。$\lambda$ の値が大きくなるにつれて、分布の形状は右歪みから対称に近づきます。

（Fig1. ポアソン分布の確率質量関数（λ = 1, 3, 5, 10））

ポアソン分布が適切なモデルとなるための仮定として、各観測単位内での事象の発生が互いに独立であること、単位区間内での発生率 $\lambda$ が一定であること、同一時点での複数事象の同時発生確率が無視できることが要求されます。ポアソン分布は連続値や負値をとる応答変数には適用できません。また実データにおいては、分散が平均を超える過分散や、ゼロカウントが過剰に発生するゼロ過剰が観察される場合があり、その際はポアソン分布の等分散仮定が成立しなくなります。

ポアソン回帰モデルの定式化

ポアソン回帰は、一般化線形モデルの枠組みで定式化される統計モデルです。一般化線形モデルは確率分布、線形予測子、リンク関数の三要素から構成されます。ポアソン回帰では、応答変数 $Y_i$ の条件付き分布をポアソン分布と仮定し、条件付き期待値 $\mu_i = E[Y_i \mid \mathbf{x}_i]$ と線形予測子 $\eta_i = \mathbf{x}_i^\top \boldsymbol{\beta}$ を対数リンク関数によって結びつけます。対数線形モデルとして

$$\log(\mu_i) = \mathbf{x}_i^\top \boldsymbol{\beta}$$

が成立し、期待値の逆変換は

$$\mu_i = \exp(\mathbf{x}_i^\top \boldsymbol{\beta})$$

で与えられます。対数リンク関数を用いることで、$\mu_i$ が必ず正値をとることが自然に保証されます。線形予測子 $\mathbf{x}_i^\top \boldsymbol{\beta}$ は任意の実数値をとり得ますが、指数変換により正値の期待カウントが得られる構造です。

対数リンクはポアソン分布の正準リンクでもあります。ポアソン分布の自然パラメータは $\eta = \log(\lambda)$ であり、この対応から対数リンクは指数型分布族の理論的枠組みと整合します。対数スケールで線形加法的に表現された説明変数の効果は、元のスケールでは乗法的な効果として期待値に反映されます。すなわち、ある説明変数が $\Delta x$ だけ変化したとき、期待カウントは $\exp(\beta_j \Delta x)$ 倍になります。

（Fig2. ポアソン回帰の当てはめ曲線と95%信頼帯（logリンクによる指数的増加））

ポアソン回帰の仮定として、応答変数 $Y_i$ が条件付きポアソン分布に従うこと、各観測が他の観測に対して条件付き独立であること、リンク関数が正しく特定されていることが要求されます。等分散仮定 $E[Y_i \mid \mathbf{x}_i] = \text{Var}[Y_i \mid \mathbf{x}_i] = \mu_i$ は強い制約であり、実データでの違反が多いことが主要な限界です。また説明変数の効果が対数スケールで線形加法的であることを仮定しており、交互作用や非線形効果は明示的にモデルに組み込む必要があります。

最尤推定と係数のレート比解釈

ポアソン回帰のパラメータ $\boldsymbol{\beta}$ は最尤法によって推定されます。$n$ 個の独立観測 $(y_i, \mathbf{x}_i)$ に対して、対数尤度関数は

$$\ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left[ y_i \log(\mu_i) – \mu_i – \log(y_i!) \right]$$

と表されます。ここで $\mu_i = \exp(\mathbf{x}_i^\top \boldsymbol{\beta})$ です。最尤推定量 $\hat{\boldsymbol{\beta}}$ はスコア方程式

$$\frac{\partial \ell(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = \sum_{i=1}^{n} (y_i – \mu_i) \mathbf{x}_i = \mathbf{0}$$

を解くことで得られます。この方程式には閉形式解が存在しないため、IRLS（反復重み付き最小二乗法）と呼ばれる反復アルゴリズムにより数値的に解かれます。IRLSの詳細な手続きは後続記事で扱います。

推定された係数 $\hat{\beta}_j$ はレート比へ変換して解釈されます。レート比は

$$RR_j = \exp(\hat{\beta}_j)$$

で定義されます。連続変数 $x_j$ の場合、$x_j$ が1単位増加したときの期待カウントの倍率を表します。カテゴリ変数の場合は、参照カテゴリに対する期待カウントの比として解釈されます。統計的推論にはWald検定統計量 $z_j = \hat{\beta}_j / \text{SE}(\hat{\beta}_j)$ および $\exp\!\bigl(\hat{\beta}_j \pm 1.96 \cdot \text{SE}(\hat{\beta}_j)\bigr)$ による95%信頼区間が用いられます。点推定のみに依拠することは推定不確実性を無視することになるため、信頼区間の提示が統計的推論において必要です。

大標本条件下では、最尤推定量 $\hat{\boldsymbol{\beta}}$ は真値 $\boldsymbol{\beta}_0$ に対して一致性と漸近正規性を持ちます。小標本環境では漸近近似に基づく信頼区間や検定の精度が低下する可能性があります。カウントがほぼ0しか存在しない共変量パターンでは対数尤度が平坦化し、推定が不安定になる場合があります。

オフセット項：曝露の違いを補正する

異なる観察期間や母集団規模を持つデータでは、観測されたカウント数を直接比較することは適切ではありません。観察期間が2年の集団で10件の事象が観察された場合と、観察期間が1年の集団で10件が観察された場合とでは、年間発生率として解釈すると大きく異なります。こうした曝露量の違いを補正するためにオフセット項が導入されます。曝露変数の例として、人年、観察期間、母集団規模があります。

オフセット項 $\log(t_i)$ をモデルに加えると

$$\log(\mu_i) = \log(t_i) + \mathbf{x}_i^\top \boldsymbol{\beta}$$

となります。ここで $t_i$ は観察期間・人年などの曝露変数です。この定式化により単位曝露あたりの発生率は

$$\text{rate}_i = \frac{\mu_i}{t_i} = \exp(\mathbf{x}_i^\top \boldsymbol{\beta})$$

と表され、説明変数の効果が発生率に対する乗法的効果として解釈されます。オフセット項は係数が1に固定された共変量とみなすことができ、推定の対象とはなりません。人口規模や観察期間が異なる集団間での発生率の公平な比較において、オフセット項の適切な設定が不可欠です。

オフセット変数 $t_i$ は既知であり、モデルパラメータとして推定されないことが前提となります。観察期間中の真の発生率が一定であることも仮定されます。オフセット変数の誤設定（対数変換の失念など）は系統的な推定偏りをもたらします。観察期間内に発生率が変化する場合には、時間分割法や生存時間分析への拡張が必要となります。

適合度診断と等分散性の確認

ポアソン回帰の適合度を評価し等分散仮定の妥当性を診断するために、いくつかの統計量が用いられます。Pearsonカイ二乗統計量は

$$X^2 = \sum_{i=1}^{n} \frac{(y_i – \hat{\mu}_i)^2}{\hat{\mu}_i}$$

で定義されます。適合が良好な場合、$X^2$ は自由度 $n – p$ の $\chi^2$ 分布に近似的に従います。逸脱度は

$$D = 2\sum_{i=1}^{n}\left[y_i \log\!\left(\frac{y_i}{\hat{\mu}_i}\right) – (y_i – \hat{\mu}_i)\right]$$

で定義され、飽和モデルと当てはめモデルの対数尤度差を2倍したものです。

過分散の診断指標として、分散パラメータの推定量

$$\hat{\phi} = \frac{X^2}{n – p}$$

が広く用いられます。$\hat{\phi} \approx 1$ であれば等分散仮定が概ね成立していると判断されますが、$\hat{\phi} \gg 1$ は過分散の兆候を示します。真の条件付き分散がポアソン分布の分散関数 $V(\mu) = \mu$ に従うことがポアソン回帰の中心的仮定であり、この確認が診断の出発点となります。

Pearson残差プロットでは、予測値 $\hat{\mu}_i$ に対してPearson残差 $(y_i – \hat{\mu}_i)/\sqrt{\hat{\mu}_i}$ を描画します。等分散仮定が成立する場合、残差は予測値に対して無相関であり一定の散布を示すことが期待されます。過分散が存在する場合には標準誤差が過小評価され、第一種誤りの確率が膨らみます。ゼロ過剰への対処はゼロ過剰ポアソンモデルや負の二項回帰への拡張が必要であり、標準的なポアソン回帰では対応できません。

（Fig3. Pearson残差プロット：等分散モデルと過分散モデルの比較）

疫学への応用：疾患発生率の分析

疫学研究では、異なる集団間での疾患発生率の比較が中心的な分析課題となります。ポアソン回帰は、人年を分母とした疾患発生率のモデリングに広く用いられます。地域別・職業別の疾患発生数を、年齢・性別・曝露群を共変量とし人年をオフセットとしてモデリングすることで、交絡因子を制御した発生率の比較が可能となります。

集団 $i$ における疾患件数 $Y_i$ を人年 $T_i$ のもとで

$$\log(\mu_i) = \log(T_i) + \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}$$

とモデリングします。標準化発生比は観察件数と期待件数の比であり、年齢・性別構成が統制された状況でのポアソン回帰の推定値はこの標準化発生比と対応する構造を持ちます。交絡因子を共変量として制御したもとでの調整レート比 $\exp(\hat{\beta}_j)$ は曝露効果の推定量として用いられます。ただし、観察研究で得られるレート比は統計的な関連であり、因果関係を直接示すものではありません。因果推論には適切な研究デザインと交絡の制御に関する別途の議論が必要です。

各年齢層・性別・曝露群の内部で疾患発生率が一定であるという均質性仮定、および各個人の観察が独立であることが前提となります。集団レベルで集計されたデータから個人レベルの関連を推論する際には、生態学的誤謬のリスクが伴います。疫学データはクラスタリングや空間相関・不均一性により過分散を示す場合が多く、分析前に $\hat{\phi}$ を確認し、必要に応じて準ポアソンや負の二項回帰への切り替えを検討することが実務上不可欠です。

関連モデルとの比較と選択指針

ポアソン回帰の適用範囲と限界は、隣接するモデルとの対比により明確になります。通常の最小二乗法による線形回帰モデルとの違いは、応答変数の分布仮定、リンク関数、係数の解釈にあります。線形回帰は正規分布を仮定し恒等リンクを用いるため、予測値が負値をとり得ます。ポアソン回帰は対数リンクにより正値制約を自然に満たし、係数の解釈も加法的効果ではなく乗法的効果（レート比）となります。

過分散が検出された場合の対処として、準ポアソン回帰と負の二項回帰が選択肢となります。準ポアソン回帰はポアソン回帰の分散をスケールパラメータ $\phi$ で拡大することで過分散に対処しますが、完全な確率モデルではないため、AIC（赤池情報量規準）に基づくモデル比較には制約があります。負の二項回帰は個体差に由来する過分散を確率モデルとして明示的に扱い、尤度に基づく推論が可能です。ゼロカウントが過剰に多い場合には、ゼロ過剰ポアソンモデルがカウント生成過程を2段階でモデル化する枠組みを提供します。ポアソン回帰が適切に用いられる条件として、等分散性の成立、ゼロ過剰がないこと、観察の独立性が揃っていることが再確認されます。