15. ゼロ過剰モデル：ZIP・ZINBの理論と応用

ゼロ過剰とは何か

カウントデータを扱う際、観測ゼロの割合が通常のポアソン分布から期待される割合を大幅に上回ることがあります。この現象をゼロ過剰と呼び、疫学・生態学・経済学のカウントデータ分析において頻繁に遭遇します。

ゼロ過剰には2つのメカニズムが存在します。構造ゼロとは、対象が当該イベントを生起し得ない状態にあることで生じるゼロです。非飲酒者の飲酒量は必ずゼロであり、これは確率的な偶然ではなく対象の性質に起因します。標本ゼロとは、イベントが発生し得る状況にあるにもかかわらず、観測期間内に偶然ゼロが記録された場合です。この2種類のゼロを区別することが適切なモデル選択の出発点となります。

ポアソン分布$\text{Pois}(\lambda)$においてゼロが観測される確率は次式で与えられます。

$$P(Y=0) = e^{-\lambda}$$

標本平均$\bar{y}$を$\lambda$の推定値として用いたとき、観測ゼロ数と期待ゼロ数の比によりゼロ過剰の程度を評価できます。

$$\text{ZI率} = \frac{\text{観測ゼロ数}}{n \cdot e^{-\bar{y}}}$$

この比が1を大幅に超える場合、ゼロ過剰の存在が示唆されます。疫学データにおける医師受診回数や感染件数では、調査期間中に全く受診・感染しなかった個体が多数存在するため、ゼロが多発します。これらのゼロには構造ゼロと標本ゼロの両方が混在する場合があります。

通常のポアソン回帰でゼロ過剰データを当てはめると、ゼロに対応する残差が系統的に正側に偏り、過大なPearsonカイ二乗統計量が観察されます。このような残差の系統的なパターンはモデルの誤設定を示す診断的指標となります。

（Fig1. 医師受診回数データにおける観測度数分布とポアソン期待度数の比較（ゼロ過剰の視覚化））

ZIPモデルの定式化

ゼロ過剰ポアソンモデル（ZIP）は、ゼロを生成する2つのプロセスを明示的に混合することでゼロ過剰を統計的に扱います。観測値$Y$の確率質量関数は次のように定義されます。

$$P(Y=0) = \pi + (1-\pi)e^{-\lambda}$$
$$P(Y=k) = (1-\pi)\frac{e^{-\lambda}\lambda^k}{k!}, \quad k \geq 1$$

ここで$\pi \in [0,1)$は構造的に非リスク状態にある確率（not-at-risk probability）、$\lambda > 0$はリスク群における平均カウントです。ZIPは$\pi$の確率でBernoulliプロセスがゼロを生成し、$(1-\pi)$の確率でポアソンプロセスがカウントを生成する2成分混合モデルとして解釈されます。ポアソン部分もゼロを生成し得るため、ゼロ観測は2つの生成源から構成されます。

各成分のパラメータは共変量を通じて以下のリンク関数でモデル化されます。

$$\text{logit}(\pi_i) = \mathbf{z}_i^\top \boldsymbol{\gamma}$$
$$\log(\lambda_i) = \mathbf{x}_i^\top \boldsymbol{\beta}$$

ゼロ過剰成分の共変量ベクトル$\mathbf{z}$とカウント成分の共変量ベクトル$\mathbf{x}$には、それぞれ異なる説明変数を投入できます。この設計上の柔軟性により、構造ゼロを決定する要因とカウントの大きさを決定する要因を独立してモデル化できます。

ZIPの対数尤度は、ゼロ観測（$y_i=0$）と正値観測（$y_i \geq 1$）を分離した次の形で表現されます。

$$\ell(\boldsymbol{\beta}, \boldsymbol{\gamma}) = \sum_{i: y_i=0} \log\left[\pi_i + (1-\pi_i)e^{-\lambda_i}\right] + \sum_{i: y_i \geq 1} \left[\log(1-\pi_i) + y_i\log\lambda_i – \lambda_i – \log(y_i!)\right]$$

ZIPモデルには次の仮定が伴います。カウント成分はポアソン分布に従い、条件付き期待値と分散が等しい等分散条件（$E[Y|\mathbf{x}]=\text{Var}[Y|\mathbf{x}]$）が成立します。またゼロ過剰メカニズムとカウントプロセスは条件付き独立であることが前提となります。

ZIPモデルの限界として、ポアソン部分の等分散仮定は維持されるため、過分散がカウント成分に残存する場合にはZIPは不適切であり、ZINBへの拡張が必要となります。また2つの回帰方程式構造により、共変量の限界効果の解釈が通常の回帰分析より複雑になります。

ZINBモデルへの拡張

ゼロ過剰と過分散が同時に存在する場合、カウント成分に負の二項分布を採用したゼロ過剰負の二項モデル（ZINB）が適切です。ZINBの確率質量関数は次のように表されます。

$$P(Y=0) = \pi + (1-\pi)\left(\frac{\theta}{\theta+\mu}\right)^{\theta}$$
$$P(Y=k) = (1-\pi)\frac{\Gamma(k+\theta)}{k!\,\Gamma(\theta)}\left(\frac{\theta}{\theta+\mu}\right)^{\theta}\left(\frac{\mu}{\theta+\mu}\right)^{k}, \quad k \geq 1$$

ここで$\mu$はカウント成分の平均、$\theta > 0$は分散を制御する形状パラメータです。NB2型分散関数（$\text{Var}=\mu+\mu^2/\theta$）を採用するZINBの全体分散は次式で与えられます。

$$\text{Var}(Y) = (1-\pi)\left[\mu + \frac{\mu^2}{\theta}\right] + \pi(1-\pi)\mu^2$$

右辺第1項は負の二項分布由来の分散に$(1-\pi)$を乗じたもの、第2項はゼロ過剰混合に由来する追加分散です。$\theta \to \infty$の極限では$\mu^2/\theta \to 0$となり、カウント成分のポアソン収束によりZINBはZIPに一致します。したがってZIPはZINBのネスト特殊ケースとして位置づけられます。

ZINBモデルには次の仮定が伴います。カウント部分は負の二項分布に従い、形状パラメータ$\theta$はデータから推定されます。ゼロ過剰部分は独立した二値プロセスとして機能します。

ZINBモデルの限界として、パラメータ数の増加（$\boldsymbol{\beta}$・$\boldsymbol{\gamma}$・$\theta$）により小標本での収束不安定リスクが増します。また$\theta$の識別可能性はデータのゼロ以外の分布の広がりに依存するため、ゼロ以外の観測が少ない場合は$\theta$の推定が困難になります。

パラメータ推定：EMアルゴリズムと最尤法

ZIPおよびZINBのパラメータ推定は、EMアルゴリズム（期待値最大化法）または直接数値最適化によって実施されます。EM視点では、各ゼロ観測が構造ゼロか標本ゼロかは観測されない潜在状態であり、ZIP問題を不完全データ問題として捉えます。潜在指示変数$Z_i$をi番目観測が構造ゼロである場合に1、リスク群からのカウントである場合に0をとるBernoulli変数として定義します。

E-stepでは、現在のパラメータ推定値のもとで$Z_i$の事後期待値を計算します。$Y_i=0$が観測された場合の事後期待値は次式となります。

$$E[Z_i | Y_i=0] = \frac{\pi_i}{\pi_i + (1-\pi_i)e^{-\lambda_i}}$$

$Y_i \geq 1$の場合は$Z_i$は必ず0であるため、事後期待値は0となります。M-stepでは、E-stepで得られた期待値を重みとして各成分のパラメータを更新します。$\pi$の更新はゼロ観測に対する事後期待値の平均として得られ、$\lambda$（またはZINBにおける$\mu$と$\theta$）の更新は重み付き最尤推定として実施されます。この反復によりEMアルゴリズムの単調増加性が保証され、対数尤度は収束まで単調増加します。

標準誤差の算出にはLouis公式に基づく観測情報行列が用いられます。完全データ情報行列と潜在変数の欠損に由来する情報損失の差として観測情報行列が得られ、その逆行列から共分散行列が推定されます。直接数値最適化と比較した場合、EMはHessian計算を回避できますが収束速度は一般に遅くなります。

パラメータ推定の仮定として、対数尤度面が実用上一山型でありEMが実用的な解に収束することが前提となります。限界として、EMは対数尤度の単調増加を保証しますが大域的最大値は保証しないため、複数の初期値から推定を実施することが推奨されます。またHessianの直接計算が困難な場合には、標準誤差に数値誤差が混入する可能性があります。

モデル診断と検定

ゼロ過剰モデルの適切性を診断するために、Vuong検定・尤度比検定・AIC（赤池情報量規準）/BIC（ベイズ情報量規準）・Rootgramが利用されます。

Vuong検定はネストしていないモデルの比較に適用されます。ZIPとポアソン回帰の比較においては、各観測の対数尤度比の標本平均を標準化した統計量が用いられます。

$$V = \frac{\frac{1}{\sqrt{n}}\displaystyle\sum_{i=1}^{n} m_i}{\hat{\sigma}_m}, \quad m_i = \log\frac{f_{\text{ZIP}}(y_i)}{f_{\text{Poisson}}(y_i)}$$

ここで$\hat{\sigma}_m$は$m_i$の標本標準偏差です。$V$は漸近的に標準正規分布に従い、$V > 1.96$でZIP優位、$V < -1.96$でPoisson優位、中間は不確定と判定されます。Vuong検定の帰無仮説は両モデルが等価な当てはまりを持つことであり、モデルの構造的優劣ではなく当てはまりの比較に限定されます。

ZIPとZINBのようなネストモデルの比較には尤度比検定が適用されます。

$$LR = 2(\ell_{\text{ZINB}} – \ell_{\text{ZIP}}) \sim \chi^2(1)$$

この検定では$\theta \to \infty$が境界値であるため、カイ二乗近似が保守的になる場合があります。AICおよびBICはモデル選択の補助的指標として有用ですが、情報量規準はモデルの構造的適切性そのものを検証するものではなく、検定による評価と併用する必要があります。

ルートグラムはゼロを含むすべてのカウント値における適合度を視覚的に診断する手法です。各カウント値$k$に対する観測度数$y_k$の平方根$\sqrt{y_k}$と、モデルから予測された期待度数の平方根を並べることで、ゼロやその近傍のカウントへの当てはまりの差異を評価します。

（Fig2. Poisson・ZIP・ZINBモデルによる予測分布と観測度数の比較（rootgramスタイル））

（Fig3. Vuong検定統計量の漸近標準正規分布と棄却域（ZIPとPoisson回帰の比較例））

Vuong検定の仮定として、各観測の対数尤度差$m_i$が漸近正規性を持つことが前提となります。限界として、両モデルが共に誤設定の場合に偽陽性が発生し得ます。また小標本では漸近理論に基づく棄却域が不正確になるため、十分な標本サイズが確保されていない状況での適用には注意が必要です。

疫学データへの応用

疫学研究における医師受診回数データのモデリングはZIPの典型的な適用場面です。健康な非受診者（構造ゼロの候補）と観察期間内に偶然未受診だった個体（標本ゼロ）が混在するため、通常のポアソン回帰では適切な推定が困難です。

ZIPモデルの係数は2成分それぞれで独立して解釈されます。ゼロ過剰成分の係数$\gamma_j$については、対応する共変量が1単位増加したときの構造ゼロとなるオッズの変化率として解釈されます。

$$\text{構造ゼロとなるオッズ比} = e^{\gamma_j}$$

カウント成分の係数$\beta_j$については、リスク群（受診者）内における平均受診回数の比として解釈されます。

$$\text{レート比} = e^{\beta_j}$$

例えば慢性疾患の有無は非受診者確率（構造ゼロ確率）に影響する要因として$\boldsymbol{\gamma}$で捉えられ、疾患重症度は受診者の受診頻度に影響する要因として$\boldsymbol{\beta}$で捉えられます。この2成分係数の独立した解釈が、単一回帰式では分離できない医療行動のパターンを明示化する実践的価値をもたらします。

周辺効果の計算では、$Y$の周辺期待値$E[Y]=(1-\pi)\lambda$に対する偏微分として導出されます。

$$\frac{\partial E[Y]}{\partial x_j} = (1-\pi)\lambda\beta_j + (1-\pi)\lambda\left(-\frac{\partial \pi}{\partial x_j}\cdot\frac{1}{1-\pi}\right)$$

ゼロ過剰成分とカウント成分の両方が周辺期待値に寄与するため、限界効果は通常の回帰係数とは異なる解釈が必要です。予測値の報告にあたっては、構造ゼロ確率$\hat{\pi}$と受診頻度$\hat{\lambda}$を分離して可視化し、それぞれの信頼区間を提示することが推奨されます。

この適用における限界として次の点が挙げられます。横断研究設計では構造ゼロと標本ゼロの判別が不可能であり、$\pi$の因果的解釈は調査設計の制約に依存します。未測定交絡変数がゼロ過剰成分とカウント成分の両方にバイアスをもたらす可能性があります。また健康な非受診者と保険未加入等の意図的非受診者を観測データから区別することはできません。

モデル比較と拡張の展望

カウントデータに対するモデルの選択は、過分散の有無とゼロ過剰のメカニズムの解釈によって決まります。以下の比較表に各モデルの特性をまとめます。

HurdleモデルはZIPと類似しますが構造的に異なります。Hurdleモデルでは、ゼロはすべて二値プロセスから生成され、正カウントは切断ポアソンまたは切断負の二項分布から生成されます。ZIPではポアソン部分もゼロを生成し得るのに対し、Hurdleでは正カウントへの移行に一度の障壁（hurdle）が設定されます。どちらのモデルが適切かの判断基準は、ゼロに対して実質的な構造的意味付けが可能か否かに基づきます。

クラスター相関がある場合（例：同一病院内の患者間相関）には、ランダム効果を取り込んだゼロ過剰混合効果モデルへの拡張が必要です。ベイズZIPモデルでは$\pi$と$\lambda$の事後分布から不確実性を定量化でき、小標本や事前情報が豊富な状況での推定の安定性に利点があります。連続変数への拡張としてゼロ過剰ガンマモデルやゼロ過剰ベータモデルが、医療費（ゼロが多いが正値は連続）のような半連続データに適用されます。

より複雑なモデル（ZINB・混合効果ZIP等）は解釈可能性とのトレードオフが大きくなります。データが少ない場合は複雑なモデルで推定が不安定になり、単純な負の二項回帰モデルが実用的に優位な場合があります。モデルの複雑さはデータの特性と標本の大きさのバランスによって決定する必要があります。