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

はじめに

以前、ゼロの多い離散データの統計モデルの話を書きましたが、今回はその続きとして、連続型データの場合のモデルを紹介していきたいと思います。
mstour.hatenablog.com

離散データの場合には飲酒回数を例として考えました。ゼロの多い連続データの典型例としても、ある期間における飲酒の量が挙げられます。
ただし、こちらは飲んだ「頻度」ではなく「量（mlなどの単位で表される）」が該当します。
全くお酒を飲まない（飲めない）人の飲酒量は常にゼロのため、ヒストグラムを描くとゼロに高い山ができます。
一方、飲酒の習慣がある人についても、飲む量には個人間のばらつきが大きいと考えられ、非常に右に裾の長い分布形になることが考えられます。このような分布形のデータは「Semi-continuous data」と呼ばれます。

その他には、

エクササイズにかける時間（運動を全くしない人の時間はゼロ）
喫煙本数（喫煙しない人の本数はゼロ）

などもゼロの多い連続データと言えるでしょう。

このようなデータをうまく表現する方法の1つとして、「ゼロのパート」と「正の値のパート」の2つの要素からなる「Two-part model」があります。
ゼロのパートに関しては「ゼロでない確率 $\pi$ *1をパラメータとする二項分布」、正の値のパートに関しては「正の値のみをとる連続型確率分布」が具体的な構成要素になります。

他には離散型の場合と同様にZero-inflated modelも選択肢になりますが、手元に詳しい資料がないので今回は割愛します。こちらも離散型データの時と考え方は同じですが、連続型分布として「0以上の値をとる」分布を用意する必要があります。0を含む分布は基本的な確率分布の中ではあまり思い当たらないのですが、切断分布(Truncated distribution)を使えば問題ありません。
切断正規分布を用いた例として、例えば以下のようなものがあります。
www.ncbi.nlm.nih.gov

ということで、今回はTwo-part modelに関して話を進めていきたいと思います。

Two-part model

Semi-continuus data $Y$ に対するTwo-part modelは、次のような確率密度関数で表現できます。

$\displaystyle f(y) = \left( 1-\pi \right) ^{(1-s)} \times \left[ \pi g(y) \right]^s \hspace{10pt} (y \geq 0)$

なお $\pi$ は $Y \neq 0$ になる確率、 $s$ は $Y = 0$ のときに0、 $Y \neq 0$ のときに1をとる変数とします。
また $g(y)$ は正の値をとる何らかの連続型確率分布の確率密度関数です。選択肢としては対数正規分布やガンマ分布、さらに複雑な分布も考えることができます。

Two-part modelは、離散型データのモデルとして紹介したHurdle modelと同じ形式になっています。
上の定義の通り、確率 $1-\pi$ で値が0になるのですが、この0の値は「常に0になる」場合と「たまたま0になった」場合との両方を含んでいます。
したがって、常に値が0になる人（例えば、お酒を全く飲まない人）と、たまたま値が0だった人（お酒を飲むことはあるが、調査期間中にたまたま全く飲まなかった人）とを区別することができません。
一方、Zero-inflated modelでは「お酒を全く飲まない人」と「たまたま飲酒量が0だった人」の確率分布はモデル内の別々の部分で表されています。ここがTwo-part modelとZero-inflated modelの大きな相違点かと思います。

回帰

離散型データと同様にして回帰モデル(Two-part regression model)を考えることができます。
以下では、モデルの数式を具体的に特定するために、使用する連続型分布として対数正規分布を選択することにします。
対数正規分布は、 $\mu$ と $\sigma$ の2つのパラメータで特徴付けられます。
対数正規分布に従う変数を $Y$ とすると、 $\mu$ は対数変換した時の平均 $E(\log(Y))$ 、 $\sigma^2$ は同じく対数変換後の分散 $V(\log(Y))$ にあたります。

回帰モデルでは例えば以下のように、平均パラメータ $\mu$ と、Two-part modelにおいて値が正になる確率 $\pi$ それぞれをいくつかの説明変数で表現します：

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

離散型データの場合と同様、説明変数は $\mu$ に関するモデルと $\pi$ に関するモデルとで同じでも異なっていても構いません。

なお平均パラメータ $\mu$ は正の値に関する分布（ここでは対数正規分布）に関するパラメータであることに注意する必要があります。
上の式でも出てきているように、このパラメータは「目的変数の値がゼロではない（正の値）という条件のもとでの（言い換えると、 $i=1,\cdots,N$ のうち $Y_i>0$ だった人に限った）」（対数変換後の）平均を意味しています。
したがってj番目の回帰係数 $\beta_{2j}$ は、「目的変数の値がゼロでない人たちに対する」説明変数 $x_{i2j}$ の効果を示しており、値がゼロだった人は対象に含まれていません。

推測の対象を値がゼロだった人も含む母集団全体へ広げるためには、これも離散型の場合と同様に、母集団平均 $E(Y_i)$ を直接モデリングする「Marginalized model」を考えます。

Marginalized modelでは、値が正になる確率 $\pi$ に関する部分は通常のモデルと同じですが、以下のように母集団平均と説明変数とを関連づけます。

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

目的変数の値は0以上であるため、右辺を指数関数で変換しています。
このモデルを用いると、他の説明変数の値が同じである場合、

$\displaystyle \frac{E(Y_i|x_{i2j}=x+1)}{E(Y_i|x_{i2j}=x)} = \frac{\exp(\alpha_{2j} x) \exp(\alpha_{2j})}{\exp(\alpha_{2j} x)} = \exp(\alpha_{2j})$

より、「j番目の説明変数 $x_{i2j}$ が1単位増加したとき、母集団平均は $\exp(\alpha_{2j})$ 倍になる」と解釈することができます*2。

f:id:mstour:20210904144730j:plain

Marginalized modelの回帰係数を推定するためには、回帰係数と対数正規分布のパラメータ $\mu_i$ とを結びつける式を導出して、尤度へ代入する必要があります。
次の式の(1)はMarginalized modelの定義式、(2)はTwo-part modelの確率密度関数を用いて $Y_i$ の期待値を求めています。 $g(y)$ はここでは対数正規分布の確率密度関数です。

$\displaystyle (1) E(Y_i) = \exp \left( \alpha_{20} + \alpha_{21} x_{i21} + \cdots + \alpha_{2q} x_{i2q} \right) \\ (2) E(Y_i) = \int_0^{\infty} y_i f(y_i) dy_i = \pi_i \int_0^{\infty} y_i g(y_i) dy_i = \pi_i \exp \left( \mu_i + \frac{\sigma^2}{2} \right)$

(1)(2)それぞれの右辺が等しいので、 $\mu_i$ について整理すると

$\displaystyle \mu_i = \alpha_{20} + \alpha_{21} x_{i21} + \cdots + \alpha_{2q} x_{i2q} - \log(\pi_i) - \frac{\sigma^2}{2}$

となり、Two-part modelの尤度を回帰モデルのパラメータで表すことができます（ $\pi_i$ の部分にも前に定義したモデルを代入します）。

おわりに

0の多い連続型データに対しても、離散型データの場合の「Hurdle model」や「Zero-inflated model」と同様の統計モデルを用いることができます。
今回はそのうちの一つ「Two-part model」を紹介しました。
正規分布などの標準的な分布を無理に当てはめたり、データの変数変換やカテゴリ化を行うことなく、現象の振る舞いをより自然にモデルで表現できることが大きな魅力かと思います。
なお正の値の確率分布としては対数正規分布のほか「Skew lognormal distribution（対数歪正規分布？）」や「一般化ガンマ分布」などが提案されています。

参考文献

[1] Nguyen, CD, Moreno-Betancur, M, Rodwell, L, Romaniuk, H, Carlin, JB, Lee, KJ. Multiple imputation of semi-continuous exposure variables that are categorized for analysis. Statistics in Medicine. 2021; 1– 14.
[2] Smith, V. A., Preisser, J. S., Neelon, B., and Maciejewski, M. L. (2014), A marginalized two-part model for semicontinuous data, Statist. Med., 33, 4891– 4903.
[3] 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.
[4] 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.
[5] Zhang, X., Guo, B., & Yi, N. (2020). Zero-Inflated gaussian mixed models for analyzing longitudinal microbiome data. Plos one, 15(11), e0242073.

*1:「ゼロのパート」と言っていますが、解釈のしやすさと参考文献の表記を踏まえて「ゼロでない確率」を以降では考えます。

*2:Marginalizedでない回帰モデルの場合にはこのような解釈はできず、説明変数の効果は他の説明変数の水準に依存するという特徴があります。