12. 負の二項回帰：カウントデータの過分散対応

負の二項分布の確率モデル：ポアソン-ガンマ混合

カウントデータに対する回帰分析の基本モデルであるポアソン回帰は、条件付き分散が条件付き平均に等しいという等分散性を前提とします。しかし医療・疫学・交通安全などの実データでは、観測された分散が平均を上回る過分散が頻繁に生じます。過分散を無視してポアソン回帰を適用すると、標準誤差が過小になり、仮説検定や信頼区間が信頼できなくなります。負の二項回帰は、この過分散を確率モデルの構造として明示的に組み込んだ一般化線形モデルの一種です。

負の二項分布（NB2型）は、ポアソン分布とガンマ分布の混合として導出されます。各観測単位 $i$ の真の発生率 $\lambda_i$ が、平均 $\mu_i$・形状パラメータ $\theta>0$ のガンマ分布に従うと仮定します。

$$\lambda_i \sim \text{Gamma}\!\left(\theta,\,\frac{\theta}{\mu_i}\right)$$

$\lambda_i$ を所与としたとき $Y_i$ はポアソン分布 $\text{Poisson}(\lambda_i)$ に従います。$\lambda_i$ を周辺化（積分消去）すると、$Y_i$ の周辺分布はガンマ関数 $\Gamma(\cdot)$ を含む次の確率質量関数をもちます。

$$P(Y=y \mid \mu,\theta) = \frac{\Gamma(y+\theta)}{\Gamma(\theta)\,y!}\left(\frac{\theta}{\theta+\mu}\right)^{\theta}\left(\frac{\mu}{\theta+\mu}\right)^{y}, \quad y=0,1,2,\ldots$$

平均と分散はガンマ混合の性質から次のように導かれます。

$$E[Y]=\mu, \qquad \text{Var}[Y]=\mu+\frac{\mu^2}{\theta}$$

分散の第2項 $\mu^2/\theta$ が過分散の構造的起源です。形状パラメータ $\theta$ が大きくなるほど過分散は小さくなり、$\theta \to \infty$ の極限では $\text{Var}[Y] \to \mu$ となってポアソン分布に収束します。この収束性はモデル選択においてポアソン分布と NB2 の連続的な関係を与えます。

潜在的なランダム効果 $\lambda_i$ は観測個体間の測定されていない不均一性（未観測異質性）を表す確率的機構です。個体ごとにリスクが異なる疾患研究や、地域ごとに曝露水準が異なる事故データなどに対して、この機構は分散増大の自然な説明を与えます。負の二項分布は離散型の指数型分布族に属し、対数リンク関数を介した一般化線形モデルとして定式化されます。NB2 モデルは潜在的なランダム効果がガンマ分布に従うことを仮定しており、この仮定が成立しない場合には推定量の安定性と解釈に影響が生じます。

（Fig1. ポアソン分布とNB2分布の確率質量関数の比較（μ=5、θ=2, 5, 20））

NB1とNB2：二つの分散パラメタリゼーション

負の二項分布には分散関数の異なる複数のパラメタリゼーションが存在します。実務および文献で最も頻繁に用いられるのは、分散が平均の二次関数となる NB2 と、線形となる NB1 です。両者の選択はデータの過分散構造の理解に直結します。

NB2 の分散関数は次のとおりです。

$$\text{Var}[Y]=\mu+\frac{\mu^2}{\theta}$$

NB1 の分散関数は、過分散パラメータ $\phi>0$ を用いて次のように表されます。

$$\text{Var}[Y]=\mu(1+\phi)$$

これら両者を包含する統一的な表現として、次の一般化分散関数があります。

$$\text{Var}[Y]=\mu+\alpha\mu^{p}$$

$p=2,\,\alpha=1/\theta$ のとき NB2、$p=1,\,\alpha=\phi$ のとき NB1 に対応します。過分散パラメータ $\alpha$（または $1/\theta$）は非負でなければなりません。$\alpha<0$ は確率的に無意味であり、この制約が成立しない場合はモデルの適用範囲外となります。

NB2 が一般化線形モデルの標準的な負の二項モデルとして広く採用される理由は主に二点あります。第一に、ポアソン-ガンマ混合という確率論的な基礎付けをもつことです。第二に、分散が平均に対して二次的に増大する構造が、生物・医学・社会科学データの実態と整合することが多い点です。一方 NB1 は、分散と平均の比が一定（$\text{Var}/\mu=1+\phi$）であるという線形過分散の状況に適しており、製造業の不良品数や特定の社会科学データでの適用例があります。

分散-平均プロットはモデル選択の補助的指標として有用です。観測データの分散と平均の関係を対数スケールで示したとき、両対数プロット上の傾きが1に近ければ NB1 の構造を、2に近ければ NB2 の構造を示唆します。小標本では NB1 と NB2 の識別が困難になります。分散関数の形状は漸近理論に基づいて推定されるため、標本サイズが小さい場合には両者の情報量規準の差が統計的に有意にならず、モデル選択の根拠が弱くなります。

（Fig2. 分散-平均関係の比較：ポアソン・NB1・NB2）

最尤推定とスコア方程式

負の二項回帰の回帰係数 $\beta$ と形状パラメータ $\theta$ は最尤推定によって求められます。$n$ 個の独立な観測 $(y_i,\mathbf{x}_i)$ に対して、NB2 の対数尤度は次のように定式化されます。

$$\ell(\beta,\theta \mid \mathbf{y}) = \sum_{i=1}^{n}\Bigl[\log\Gamma(y_i+\theta)-\log\Gamma(\theta)-\log(y_i!)+y_i\log\mu_i-(y_i+\theta)\log(\theta+\mu_i)+\theta\log\theta\Bigr]$$

ただし $\mu_i=\exp(\mathbf{x}_i^\top\beta)$ です。対数リンク関数を通じた連鎖律により、$\beta$ に関するスコア方程式は次のように得られます。

$$\frac{\partial\ell}{\partial\beta}=\sum_{i=1}^{n}\frac{y_i-\mu_i}{\theta+\mu_i}\,\mathbf{x}_i=\mathbf{0}$$

$\theta$ に関するスコア方程式には、ガンマ関数の対数微分であるディガンマ関数 $\psi(\cdot)$（$\psi(x)=d\log\Gamma(x)/dx$）が現れます。

$$\frac{\partial\ell}{\partial\theta}=\sum_{i=1}^{n}\left[\psi(y_i+\theta)-\psi(\theta)+\log\frac{\theta}{\theta+\mu_i}+1-\frac{\theta+y_i}{\theta+\mu_i}\right]=0$$

Fisher 情報行列は $\beta$ と $\theta$ のブロック構造をもちます。$\beta$ に関するブロックは IRLS（反復重み付き最小二乗法）の重み行列に対応し、$\theta$ に関するブロックにはポリガンマ関数 $\psi^{(1)}(\cdot)$（トリガンマ関数）が含まれます。この行列の正定値性が推定量の漸近正規性の基礎となります。

実際の推定では $\beta$ と $\theta$ を交互に更新するアルゴリズムが用いられます。$\theta$ を固定した上で $\beta$ を IRLS により更新し、次いで $\beta$ を固定した上で $\theta$ を Newton-Raphson 法で更新するサイクルを繰り返します。収束判定にはスコアベクトルのノルムや対数尤度の変化量が用いられます。初期値の設定や収束速度はデータの過分散の程度に依存します。

正則条件として $\theta>0$ および観測の独立性が要求されます。大標本漸近理論に依存するため、小標本では推定量の標準誤差が実際の不確かさを十分に反映しない場合があります。過分散パラメータの正確な特定には十分な標本数が必要であり、$n<100$ 程度の小標本では信頼区間の被覆確率が名目水準を下回ることが報告されています。$\theta$ が非常に大きい（ポアソン分布に近い）状況では対数尤度の $\theta$ 方向の曲率が小さくなり、$\theta$ の推定精度が低下します。

モデルの仮定と診断

負の二項回帰の統計的仮定が成立しているかは、残差と逸脱度を用いた診断によって評価されます。残差の分布やパターンを通じて、モデルの適合度・分散構造の妥当性・ゼロ過剰の有無を確認できます。

逸脱度残差 $d_i$ は、飽和モデルの対数尤度への $i$ 番目の寄与 $\ell_s(y_i)$ と当てはまりモデルの寄与 $\ell_f(\hat{\mu}_i)$ を用いて次のように定義されます。

$$d_i=\text{sign}(y_i-\hat{\mu}_i)\sqrt{2\bigl[\ell_s(y_i)-\ell_f(\hat{\mu}_i)\bigr]}$$

Pearson 残差 $r_i$ と過分散統計量 $\hat{\phi}$ は次のとおりです。$p$ は推定されたパラメータ数です。

$$r_i=\frac{y_i-\hat{\mu}_i}{\sqrt{\hat{\mu}_i+\hat{\mu}_i^2/\hat{\theta}}}, \qquad \hat{\phi}=\frac{\displaystyle\sum_{i=1}^{n}r_i^2}{n-p}$$

$\hat{\phi}\approx 1$ であればモデルの分散特定が適切な可能性を示し、$\hat{\phi}\gg 1$ は残存する過分散を示唆します。この統計量は NB2 特定後の残留過分散のチェックに用いられます。

ゼロ過剰の検査として、観測されたゼロの個数 $n_0$ と NB2 モデルによる期待ゼロ数の比較が有効です。

$$\hat{n}_0^{\text{NB2}}=\sum_{i=1}^{n}P(Y_i=0\mid\hat{\mu}_i,\hat{\theta})=\sum_{i=1}^{n}\left(\frac{\hat{\theta}}{\hat{\theta}+\hat{\mu}_i}\right)^{\hat{\theta}}$$

$n_0 \gg \hat{n}_0^{\text{NB2}}$ となる場合は ZINB（ゼロ過剰負の二項モデル）への拡張が検討されます。

モデルの主要な仮定は以下の三点です。第一に、観測の条件付き独立性：共変量 $\mathbf{x}_i$ を条件付けたとき $Y_1,\ldots,Y_n$ は互いに独立です。この仮定はクラスター構造や空間相関をもつデータでは成立しない場合があります。第二に、対数リンクによる線形予測子と条件付き平均の関係：$\log(\mu_i)=\mathbf{x}_i^\top\beta$ が成立します。第三に、分散関数が NB2 の形式 $\text{Var}[Y_i\mid\mathbf{x}_i]=\mu_i+\mu_i^2/\theta$ で正しく特定されていることです。

形状パラメータ $\theta$ は個体間不均一性の代理的指標として解釈されます。$\theta$ が小さいほど過分散が大きく、観測間の不均一性が高いことを示します。この解釈はあくまで分布の構造的特徴を記述するものであり、因果的な解釈を与えるには別途の設計が必要です。残差プロットでは $d_i$ を予測値 $\hat{\mu}_i$ に対してプロットし、系統的なパターンが見られる場合は対数線形性の仮定または分散構造の誤設定を示唆します。ゼロ過剰が別途存在する場合、NB2 のみではモデルが不適合となります。また集計単位（週単位と月単位など）が異なるデータを混在させると $\theta$ の推定値が変化し、モデル間の比較が困難になります。

疫学への応用：発生率比の推定

疾患発生数・交通事故件数・入院件数等に対するリスク因子の発生率比（IRR）推定は、負の二項回帰の主要な応用領域です。観察期間や曝露量が個体間で異なる場合、オフセット項を含むモデルが用いられます。

観察人年 $t_i$ を考慮したオフセット項付きモデルは次のように定式化されます。

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

$\log(t_i)$ はオフセット項であり、その係数は1に固定されます。これにより $\mu_i/t_i$ が単位時間あたりの発生率として解釈されます。回帰係数 $\beta_j$ の解釈は IRR として与えられます。

$$\text{IRR}_j=\exp(\beta_j)$$

これは他の共変量を一定に保った上で $x_j$ が1単位増加したとき、発生率が $\exp(\beta_j)$ 倍になることを示します。$\text{IRR}_j>1$ は発生率の増加、$\text{IRR}_j<1$ は減少を意味します。信頼区間の構成にはデルタ法（Wald 区間）またはプロファイル尤度に基づく区間が用いられます。標本サイズが十分に大きい場合は両者の差は小さいですが、標本サイズが小さい場合やパラメータが境界近傍にある場合はプロファイル尤度法がより信頼性の高い区間を与えます。

（Fig3. 疫学データへのNB2回帰：発生率比（IRR）と95%信頼区間（フォレストプロット））

曝露期間 $t_i$ が全個体で独立かつ既知であることが仮定されます。この仮定が成立しない場合（追跡不能や打ち切りが多い場合など）、オフセット補正の精度が低下します。観察研究では未測定の交絡変数の存在により IRR の推定値に偏りが生じる可能性があります。未測定交絡因子が存在すると曝露の効果が過大または過小に評価されます。このため観察研究から得られた IRR は因果効果の推定値ではなく、交絡調整後の統計的関連性の指標として解釈する必要があります。また過小報告や記録漏れが存在する場合、観測ゼロの増加がゼロ過剰と混同されるリスクがあります。

準ポアソンとの比較とモデル選択

過分散への対処として、負の二項回帰と並んで用いられる手法が準ポアソンモデルです。両者の理論的位置付けと実践的選択基準を整理します。

準ポアソンは準尤度アプローチに基づき、分散関数を線形過分散 $\text{Var}[Y]=\phi\mu$（$\phi>1$）と仮定します。完全な確率モデルを指定しないため、AIC（赤池情報量規準）による比較は原理的に実施できません。一方 NB2 は完全尤度アプローチに基づくため AIC が定義され、モデル比較に活用できます。係数推定値の一致性については、平均構造 $\log\mu_i=\mathbf{x}_i^\top\beta$ が正しく特定されていれば分散関数の誤設定が生じても $\hat{\beta}$ は一致推定量となる準最尤推定の頑健性をもちます。ただし標準誤差の精度は分散設定の正確さに依存します。

ポアソン回帰 vs NB2 の選択には、帰無仮説 $H_0:\theta\to\infty$（すなわちポアソン分布）に対するスコア検定（Lagrange Multiplier 検定）が有効です。その検定統計量の概形は次のとおりです。

$$T_{\text{LM}}=\frac{\left[\displaystyle\sum_{i=1}^{n}\frac{(y_i-\hat{\mu}_i)^2-y_i}{2\hat{\mu}_i^2}\right]^2}{\displaystyle\sum_{i=1}^{n}\frac{1}{2\hat{\mu}_i^2}} \;\xrightarrow{d}\;\chi^2(1) \quad (H_0 \text{ 下})$$

当てはまりの粗い指標として、Pearson $\chi^2$ を自由度で割った過分散指標も参照されます。

$$\hat{\phi}_P=\frac{\chi^2}{n-p}=\frac{\displaystyle\sum_{i=1}^{n}r_i^2}{n-p}$$

モデル選択の実践的指針として、過分散パターンが線形に近い場合は準ポアソンまたは NB1 が、二次的増大を示す場合は NB2 が適しています。AIC による比較が必要な場合は準ポアソンは選択肢から外れます。準ポアソン：分散関数を線形（$\text{Var}=\phi\mu$）と仮定すること、NB2：ガンマ混合による構造的過分散（$\text{Var}=\mu+\mu^2/\theta$）を仮定することが両者の本質的な差異です。推論の目的が係数の点推定のみであれば準ポアソンで十分な場合があります。一方、モデル比較・予測確率の計算・ゼロ確率の推定を必要とする場合は、完全尤度をもつ NB2 が適しています。小標本では NB2 の $\theta$ 推定が不安定になるため、準ポアソンの方が安定した結果を与える場合があります。

ゼロ過剰・拡張モデルへの接続

NB2 モデルでも対応が困難な残留問題として、観測ゼロが期待値を大幅に超えるゼロ過剰があります。ゼロ過剰は疾患の非罹患群・事故の無発生期間・生態学的な不在データなど、データ生成機構が二段階になっている状況で生じます。この場合、ZINB（ゼロ過剰負の二項モデル）が候補として検討されます。

ZINB の確率質量関数の概形は次のとおりです。

$$P(Y=0)=\pi+(1-\pi)\,P_{\text{NB}}(0\mid\mu,\theta)$$

$$P(Y=y)=(1-\pi)\,P_{\text{NB}}(y\mid\mu,\theta), \quad y=1,2,\ldots$$

$\pi\in[0,1)$ はゼロ成分の混合比率であり、$P_{\text{NB}}(0\mid\mu,\theta)=(\theta/(\theta+\mu))^{\theta}$ は NB2 分布のゼロ確率です。構造ゼロ（真の事象発生が原理的にゼロである状況）と標本ゼロ（偶然に観測されなかったゼロ）を区別することが ZINB の動機となります。この区別が曖昧な場合、モデルの推定と解釈が困難になります。

NB 分布と Tweedie 分布族の関係も位置付けとして重要です。Tweedie 分布族のパラメータ $p\in(1,2)$ は複合ポアソン-ガンマ分布を表し、$p$ が2に近づくにつれて NB2 に近似的に対応します。この観点からは負の二項回帰は指数型分布族の広い枠組みの中に位置付けられます。連続的なカウントデータや半連続データへの拡張を検討する場合、Tweedie 分布族は NB2 の自然な拡張先の一つです。

ZINB ではパラメータ数の増加により識別可能性が低下します。混合比率 $\pi$ と平均 $\mu$ の識別には十分な標本数とゼロ過剰の明確なメカニズムが必要です。「ゼロ過剰か NB2 で十分か」の判断はゼロ確率の観測数と期待数の比較や尤度比検定によって評価されますが、小標本での検出力は限られます。解釈上は、構造ゼロのメカニズムが明確に仮定できる場合に ZINB の適用が正当化され、そうでない場合は過適合のリスクを考慮する必要があります。ゼロ過剰モデルの詳細は後続の記事で取り上げます。