ゼロ値を多く含むデータの統計モデル：離散データの場合

はじめに

今回は少し趣向を変えて、現在僕が関心を持っている（そこまでメジャーではないと思われる方面の）統計モデリングの話題を書いていこうと思います。
ほとんど覚え書きのような内容なので、まあ、こんなんもあるのか、という程度にみてもらえればと思います。詳しくは参考文献や、Wikipediaなんかで・・・

ゼロを多く含むデータ

世の中には性質上多くのゼロが含まれてしまうデータが色々とあります。
医学分野では、例えば特定の期間におけるがん検診の受診回数や外来受診回数などが挙げられます。「全く受診しない」人が一定数いることがしばしばあるため、これらはゼロを多く含むカウントデータという形式になる可能性があります。

カウントデータのモデリングにはポアソン分布や負の二項分布がよく用いられますが、関心ある事象の頻度が偶然ゼロになるだけでなく、ある属性の人はほぼ確実にゼロになる（検診や受診の例だと、健康な人や若い人は受診しない可能性が極めて高いというケースが考えられます）ような場合、これら一般的な確率分布が想定するよりもゼロの頻度がとても多くなることが考えられます*1。

上記のようなデータは「zero modified data*2」と称されることがあります。

連続型（と通常みなすことのできる）データにおいても、似たような状況が生じることがあります。
つまり、ゼロの一点の確率密度と、正の値をとる連続型分布とを組み合わせた分布型でよく表現できるようなデータがこれにあたります。
医学での例としては、ある期間における医療費の額が挙げられます。この場合も、「医療費がかかっていない」人が相当数いるという場合が十分考えられます。
連続型データの場合には、「semi continuous data」と呼ばれます。

これらzero modified dataやsemi continuous dataのモデリング方法としては、「two-part model」と呼ばれるモデルが有用であると言われています。
two-part modelとは、「ゼロ値のパート」と「ゼロでない値（正の値）のパート」の2つの要素のモデルを組み合わせたものです。
離散データに関するモデルとしては「hurdle model」と「zero-inflated model」、連続データについては「two-part semi continuous model」が例として挙げられます。
今回は離散データに用いられる2つのモデルについて整理していきます。

hurdle model

hurdle modelは、「値がゼロになるか・正の値（整数）になるか」を決める確率と、正の値の確率分布との2つの要素でデータを説明します。正の値の確率分布としては、一般的な離散型分布を1以上の値をとるように修正した「切断分布」を利用します。

具体的には、正の値に使用したい切断前の確率分布を $p(y|\theta)$ とすると

$\displaystyle P(Y = 0) = 1 - \pi, \\ P(Y = y) = \pi \frac{p(y|\theta)}{1 - p(0|\theta)} \hspace{10pt} (y = 1, 2, \cdots)$

の形でゼロおよび正の値の確率を定義します。
ここで $\pi = P(Y>0)$ 、つまり $\pi$ は、がん検診の例で言えば、検診を一度でも受ける確率を意味します。
まとめると以下のようにも書けます：

$\displaystyle P(Y = y) = (1 - \pi) I(y = 0) + \pi \frac{p(y|\theta)}{1 - p(0|\theta)} I(y > 0) \hspace{20pt} (1)$

ただし $I(y = 0), I(y > 0)$ は指示関数（カッコの中の論理式が真の場合1、偽の場合0を返す関数のこと）です。

f:id:mstour:20210818201007j:plain

$\pi$ と $1 - p(0|\theta)$ の大小によって、元の確率分布 $p(y|\theta)$ と比較してゼロが過剰か、あるいはゼロが過小かが決まります。
$\pi < 1 - p(0|\theta)$ であればゼロ過剰、 $\pi > 1 - p(0|\theta)$ であればゼロ過小であることを表現することができます。等号が成り立つ場合、元の分布に等しくなります。

正の値のモデリングに使用する確率分布 $p(y|\theta)$ としては、ポアソン分布や負の二項分布のほか、より複雑ないくつかの分布が提案されています。
通常、分散が平均よりも大きい「over dispersion」の場合にはポアソン分布よりも負の二項分布が好まれることが多いですが、ポアソン分布をhurdle model化することによって、ゼロ過剰の場合にはover dispersion、ゼロ過小の場合にはunder dispersionの性質も持つようになります。

zero-inflated model

zero-inflated modelもhurdle modelと同様、ゼロの頻度が非常に多い、つまりゼロの値をとる確率が一般的な確率分布よりも大きいような確率モデルです。
ベースとなるのはPoisson分布のような基本的な離散型分布ですが、hurdle modelとは異なり、切断分布ではなく基本的な離散型分布そのものを利用します。

一般的な定義は以下の通りです。

$\displaystyle P(Y = 0) = (1 - \phi) + \phi p(0|\theta), \\ P(Y = y) = \phi p(y|\theta) \hspace{10pt} (y = 1, 2, \cdots)$

例えば、比較的耳にすることの多いzero-inflated Poisson modelは、次のように定義されます。

$\displaystyle P(Y = 0) = (1 - \phi) + \phi e^{-\mu}, \\ P(Y = y) = \phi \frac{\mu^y e^{-\mu}}{y!} \hspace{10pt} (y = 1, 2, \cdots)$

ここで確率 $\phi$ で1、確率 $1-\phi$ で0となる変数 $Z$ を用意すると、上記の確率関数は次のようにまとめられます。

$\displaystyle P(Y = y) = (1 - \phi) I(Z = 0) + \phi \frac{\mu^y e^{-\mu}}{y!} I(Z = 1) \hspace{20pt} (2)$

このモデルの応用例を考えてみると、理解がしやすいかもしれません。

いま、zero-inflated Poisson modelによってある週の飲酒回数を表現したいものとします。
週の飲酒回数が0だった人たちの集団には、たまたまその週は一度も飲酒をしなかった人と、一切お酒を飲まない人とが混在していると考えられます。
ここで $\phi$ を「飲酒する人の割合（このデータを得るための調査において、飲酒する人が選ばれる確率、ともいえます）」とすると、 $Y=0$ となる確率は「全く飲酒しない人の割合*3（ $1-\phi$ ）」と、「普段飲酒する人が、たまたま週の飲酒回数が0である確率（ $\phi e^{-\mu}$ ）」を合計したものと解釈できます。

f:id:mstour:20210818201029j:plain

このように、zero-inflated modelでは「値が偶然ゼロとなる場合」と、「何らかの理由で必然的に値がゼロとなる場合」とを明示的に分けたモデリングになっているといえるでしょう。

回帰モデル

最後に、zero modified dataで表現された変数*4と、それらの値と相関関係のある変数（説明変数）との関係式を構成する回帰モデルについて説明します。
外来の受診の例で言うと、年齢や病歴といった変数が「受診する確率」にも「受診する回数」にも強く関連していることが想像できます。そこで、年齢や病歴を回帰モデルの説明変数とし、それらが「どの程度」受診状況に影響しているかを推測すれば有益な情報が得られるでしょう。

式(1)のhurdle modelに対しては、「正の値になる確率 $\pi_i$ 」を例えばロジスティックモデルで、「モデルに使用した離散型分布（例えばポアソン分布）の平均 $\mu_i$ 」を対数線形モデルで表します（ $i$ は個人を表すインデックスです）：

$\displaystyle \log \left( \frac{\pi_i}{1-\pi_i} \right) = \beta_{10} + \beta_{11} x_{i11} + \cdots + \beta_{1p} x_{i1p}, \\ \log (\mu_i) = \beta_{20} + \beta_{21} x_{i21} + \cdots + \beta_{2q} x_{i2q}$

ここで $x_{i11}, \cdots, x_{i1p}$ は正の値になる確率（外来の例だと「受診する確率」）に関係する説明変数、 $x_{i21}, \cdots, x_{i2q}$ は平均値*5（外来の例だと「受診する回数の平均」）に関係する説明変数です。
これらはそれぞれ異なっていても、全く同じでもかまいません。その分野の知識であったり、統計的なモデル選択基準によって決めたりすることになるかと思います。

式(2)のzero-inflated modelについても、全く同様にして回帰モデルを考えることができます。

ただ、いずれの場合も説明変数によって説明される平均パラメータ $\mu_i$ は母集団平均ではないことに注意する必要があります。 $\mu_i$ はzero inflated modelの場合は「正の値を取りうる集団（例えば、何らかの理由で絶対に受診しない人を除いた集団）の平均」、hurdle modelの場合は「その平均に比例する値」に相当するため、分析の目的によっては結果の解釈が難しい可能性があります。

そこで、母集団全体の平均を直接モデル化する方法も提案されています。例えばzero-inflated Poisson modelについては以下の"marginalized zero-inflated Poisson(MZIP) model"が知られています：

$\displaystyle \log \left( \frac{\phi_i}{1-\phi_i} \right) = \beta_{10} + \beta_{11} x_{i11} + \cdots + \beta_{1p} x_{i1p}, \\ \log (\nu_i) = \beta_{20} + \beta_{21} x_{i21} + \cdots + \beta_{2q} x_{i2q}$

ここで $\nu_i = E(Y_i) = \phi_i \mu_i$ は母集団全体の平均であり、 $\mu_i$ は通常のzero-inflated Poisson modelにおけるPoisson分布のパラメータ、つまり「正の値を取りうる集団での平均」です。
こうすることで、説明変数の母集団平均への影響を直接定量化することができます。

おわりに

今回はゼロの多い（あるいは少ない）離散型データに対する統計モデルとして用いられることのある、hurdle modelとzero-inflated modelについて紹介しました。
StanのUser's Guide(5.6 Zero-Inflated and Hurdle Models | Stan User’s Guide
)でも実行例が載っていますので、また何かのデータを使ってRとStanを用いた分析例も書いていければと思います。

参考文献

[1] Neelon, B., O'Malley, A. J., & Smith, V. A. (2016). Modeling zero‐modified count and semicontinuous data in health services research part 1: background and overview. Statistics in Medicine, 35(27), 5070-5093.
[2] Neelon, B., O'Malley, A. J., & Smith, V. A. (2016). Modeling zero‐modified count and semicontinuous data in health services research part 2: case studies. Statistics in medicine, 35(27), 5094-5112.
[3] 岩崎学, & 大道寺香澄. (2009). ゼロ過剰な確率モデルとそのテスト得点の解析への応用. 行動計量学, 36(1), 25-34.
https://www.jstage.jst.go.jp/article/jbhmk/36/1/36_1_25/_article/-char/ja/
[4] Stan User's Guide : Zero-Inflated and Hurdle Models
https://mc-stan.org/docs/2_20/stan-users-guide/zero-inflated-section.html

*1:逆に、何らかの理由でゼロの頻度が極めて少ない、というケースも考えることができます。

*2:ゼロが多い場合はzero inflatedと呼ばれることが多いですが、ゼロが少ない（zero deflated）場合も含めて「zero modified」と表現されます。

*3:全く飲酒しない人が母集団から選ばれる確率、とした方が適切かもしれません。

*4:回帰モデルでは目的変数、結果変数、応答変数などと呼ばれます。

*5:正確に言うと上記のモデルは「切断分布にする前のオリジナルな離散分布」の平均パラメータ（例えばポアソン分布のパラメータ）に関する関係式を立てており、母集団平均の直接のモデリングにはなっていません。ただ説明変数の値が大きくなると受診回数も増える、といった関係性は知ることができます。後述するように、母集団平均を直接モデリングする方法もあります。