欠測データ解析としての混合モデル

はじめに

反復測定データの解析方法として、これまでCovariance pattern modelと混合モデルを紹介してきた。いずれも、目的変数の確率分布として多変量正規分布を考え、モデルのパラメータを最尤法(Maximum likelihood method)によって推定する。

最尤法を用いるこれらの方法は特に経時測定データ(Longitudinal data)で利用されることが多いのだが、広く使われている大きな理由の一つとして「目的変数に欠測(Missing)がある場合でも、一定の条件のもとであれば適切な推測ができる」という点が考えられる。というのも、経時測定データを得るためには長期間の追跡が必要であり、特に人間を対象にしたデータ（臨床検査とか、質問票とか、分野を問わず多くのものが該当する）ではデータ収集期間の途中で対象者とコンタクトがとれなくなることは日常的に起こりうるからである*1。
f:id:mstour:20210303210050j:plain

今回は、混合モデルの「欠測データの解析」という側面を掘り下げてみたい。なお、Covariance pattern modelは厳密には混合モデルとは異なるが、ここで述べる内容は全く同じように当てはまるので、今回は混合モデルだけを念頭に置いて話を進めることにする。

覚えておきたいポイント

まず最初に、実用上はこれだけ知っていればとりあえずは問題ないと思われる事柄をまとめておく。
なお、欠測メカニズムがMCARであるとは、関心あるデータが欠測するかどうかは他のいかなる変数（情報）にも影響されない、つまり欠測は完全にランダムに起こることをいう。また欠測メカニズムがMARであるとは、欠測するかどうかは「観測されている変数だけ」に影響される状態を指す。
統計学的な欠測メカニズムに関する説明は以下の記事も読んでいただけると幸いである。
mstour.hatenablog.com

欠測メカニズムがMCARまたはMARであれば、最尤法によってモデルパラメータを正しく推定することができる*2。ただし、MARの場合には欠測するかどうかを説明できる変数を全てモデルに含める必要がある。
一方、混合モデルは、治療の効果やリスク因子の影響を表すモデルパラメータを最尤法によって推定する。
したがって、経時測定データに欠測が生じている場合でも、欠測メカニズムがMCARまたはMARと想定できるのであれば（必要な変数を含んだ）混合モデルによる解析は適切である。

1.については次のセクションで詳しく見ていくが、問題点を再度確認しておこう。
データの欠測メカニズムがMCARであれば、そもそもデータに欠測のある個人を解析から除外しても、推定結果にバイアスが入らないことが知られている。例えば平均値を計算する場合、完全にランダムにいくつかのデータを計算から外したとしても大きく値が変わることは（平均的には）ない。
一方、欠測メカニズムがMARである場合、欠測に関与する情報を無視して解析すると時に大きなバイアスが入りうる。例えば日本人全体の平均身長を見積もりたいとすると、男性のデータばかりが欠測している時に単純な平均をとると、本来期待するものとはかけ離れた結果になりそうである。もし仮に欠測は「男性であること」だけが原因で起きているとすると（例えば、一部の男性だけ測定をし忘れた）、直感的には性別の情報を考慮することによって適切な推測が行えそうに思える。
性別などの変数を考慮して解析を行うとなると、使用されることが多いのは線形回帰モデル（一般化線形モデルや、混合モデルなども含めるものとする）であるが、問題は欠測データがあるような場合にパラメータの推定結果が本当に正しいのか、ということである。
結論としては、最尤法で推定する場合にはパラメータ推定結果は正しい*3ことが示されている。ただし、線形モデルの中に欠測するかどうかに関係のある変数を全て含めている必要がある*4。

2.については以前の記事でも触れた通りである。

1.と2.の事実によって、3.の通り欠測のある経時測定データの解析には混合モデルが妥当な方法であることが保証される*5。
ただし、妥当であるためには欠測に関与する変数（モデルの説明変数と、途中まで観測されている目的変数との両方を指す）がモデルに含まれていることが必要となる。

欠測メカニズムがMARの時に最尤法が適切であることの詳細

ここからは、MARの時に最尤法による推定が正当化されることの理論的背景をざっと確認していく。
ここではN人の個人に対してJ回の測定が計画された（が、一部は欠測が生じている）経時測定データ $\boldsymbol{y} = (y_1, \cdots, y_J)^T$ を解析することを考えよう*6。
これらのデータの欠測状況を示す変数（以下、「欠測指標」という）を $r_j$ で表し、 $y_j$ が観測されていれば $r_j = 1$ 、欠測となっていれば $r_j = 0$ となるよう定義する。J個の経時測定データ全ての欠測指標をまとめて $\boldsymbol{r} = (r_1, \cdots, r_J)^T$ とする。データが欠測するかどうかは事前にはわからないため、欠測指標も確率変数となる。

欠測データ解析に必要な尤度を準備する

経時測定データの確率分布*7を $f(\boldsymbol{y}|\boldsymbol{\theta})$ とする。やや抽象的な表現になっているが、正規分布を仮定するモデルを考える場合には、 $\boldsymbol{\theta}$ は平均ベクトル $\boldsymbol{\mu}$ と共分散行列 $\Sigma$ の各成分を並べたベクトルになる。
通常であれば、N人分の確率分布 $f(\boldsymbol{y}|\boldsymbol{\theta})$ の積から導出される尤度関数を最大化することによってモデルの推測を行うことができるが、欠測データを扱うためには、さらに欠測指標も合わせた確率分布を考える必要がある。

そこで、経時測定データと欠測指標との同時確率分布を $f(\boldsymbol{y}, \boldsymbol{r}|\boldsymbol{\theta}, \boldsymbol{\psi})$ とする。ただし、パラメータは経時測定データに関するパラメータ $\boldsymbol{\theta}$ と、それ以外のパラメータ $\boldsymbol{\psi}$ からなる*8。
最尤法は、確率分布（正確には確率密度関数）に実際に得られたデータの値を代入してパラメータの関数とみなした尤度関数を最大化するのであった。ところが欠測データがある場合には、一部の値が存在せず、したがって観測されたデータだけを用いて最尤法を行わなければならない。

いま、経時測定データのうち、実際に観測されたものからなるベクトルを $\boldsymbol{y}_{obs}$ 、欠測したものからなるベクトルを $\boldsymbol{y}_{mis}$ とする。また、欠測指標 $\boldsymbol{r}$ は $\boldsymbol{y}_{obs}$ と $\boldsymbol{y}_{mis}$ に対応した値をとるものとする。例えば1〜3回目の値が観測、4〜5回目の値が欠測だった場合、 $\boldsymbol{y}_{obs} = (y_1, y_2, y_3)^T$ 、 $\boldsymbol{y}_{mis} = (y_4, y_5)^T$ 、 $\boldsymbol{r} = (1, 1, 1, 0, 0)^T$ である。
欠測に対応した最尤法を行うには、経時測定データと欠測指標との同時確率分布 $f(\boldsymbol{y}, \boldsymbol{r}|\boldsymbol{\theta}, \boldsymbol{\psi})$ を「欠測したデータ $\boldsymbol{y}_{mis}$ 」に関して積分したもの：

$\displaystyle f(\boldsymbol{y}_{obs}, \boldsymbol{r}|\boldsymbol{\theta}, \boldsymbol{\psi}) = \int f(\boldsymbol{y}, \boldsymbol{r}|\boldsymbol{\theta}, \boldsymbol{\psi}) d \boldsymbol{y}_{mis} \hspace{30pt} (1)$

を最大化する。左辺に観測されたデータと欠測指標の値を代入すれば未知のパラメータの関数となるので、通常の最尤法と同様の手続きをとることができる。このようにして構成した尤度を「完全尤度(Full likelihood)」あるいは「観測データの尤度(Observed data likelihood)」という。

式(1)の完全尤度を用いることで、観測されたデータだけから推定を行うことができそうな気がするが、実際には同時確率分布 $f(\boldsymbol{y}, \boldsymbol{r}|\boldsymbol{\theta}, \boldsymbol{\psi})$ の具体的な形はわからないので（少なくとも混合モデルで想定するのは目的変数の確率分布だけである）、完全尤度をさらに分解していく必要がある。
その前に、後の展開で必要となるMARの数学的定義を確認しておこう。

MARの数学的定義

欠測指標 $\boldsymbol{r}$ の確率分布について以下が成り立つとき、データ $\boldsymbol{y}$ の欠測メカニズムはMARである。

$\displaystyle f(\boldsymbol{r}|\boldsymbol{y}, \boldsymbol{\theta}, \boldsymbol{\psi}) = f(\boldsymbol{r}|\boldsymbol{y}_{obs}, \boldsymbol{\theta}, \boldsymbol{\psi}) \hspace{30pt} (2)$

つまり、 $\boldsymbol{y}$ の各要素が欠測するかどうかの確率は、観測されたデータ $\boldsymbol{y}_{obs}$ だけに影響されることを意味する。
なお、目的変数だけではなく、さまざまな変数も欠測確率に影響を与えるようなことがありえる（むしろ、ほぼそうであると言ってよいだろう）。そのような変数の集まりを $\boldsymbol{x}$ とすると、もし $\boldsymbol{y}_{obs}$ だけでMARが成立しなくても

$\displaystyle f(\boldsymbol{r}|\boldsymbol{y}, \boldsymbol{x}, \boldsymbol{\theta}, \boldsymbol{\psi}) = f(\boldsymbol{r}|\boldsymbol{y}_{obs}, \boldsymbol{x}, \boldsymbol{\theta}, \boldsymbol{\psi})$

であれば変数の集合 $\boldsymbol{y}, \boldsymbol{x}$ に関してMARが成立する。このような場合には、例えば線形モデルの説明変数に $\boldsymbol{x}$ を含めればよい。

完全尤度の分解

では、再度、式(1)の完全尤度を見てみよう。なお以降では説明変数 $\boldsymbol{x}$ の記載は省略するが、 $\boldsymbol{x}$ で条件付けしなければMARとならないような場合にはモデルに含める必要がある。その時には、以下の展開ではそれぞれの確率分布の条件部に $\boldsymbol{x}$ が加わることになる。

条件付き確率の定義を用いると、式(1)の右辺は次のように変形できる。

$\displaystyle \int f(\boldsymbol{y}, \boldsymbol{r}|\boldsymbol{\theta}, \boldsymbol{\psi}) d \boldsymbol{y}_{mis} = \int f(\boldsymbol{y}|\boldsymbol{\theta}) f(\boldsymbol{r}|\boldsymbol{y}, \boldsymbol{\theta}, \boldsymbol{\psi}) d \boldsymbol{y}_{mis}$

さらに式(2)のMARの定義を用いると

$\displaystyle \int f(\boldsymbol{y}|\boldsymbol{\theta}) f(\boldsymbol{r}|\boldsymbol{y}, \boldsymbol{\theta}, \boldsymbol{\psi}) d \boldsymbol{y}_{mis} = \int f(\boldsymbol{y}|\boldsymbol{\theta}) f(\boldsymbol{r}|\boldsymbol{y}_{obs}, \boldsymbol{\theta}, \boldsymbol{\psi}) d \boldsymbol{y}_{mis} = f(\boldsymbol{y}_{obs}|\boldsymbol{\theta}) f(\boldsymbol{r}|\boldsymbol{y}_{obs}, \boldsymbol{\theta}, \boldsymbol{\psi})$

と変形できる。まとめると

$\displaystyle f(\boldsymbol{y}_{obs}, \boldsymbol{r}|\boldsymbol{\theta}, \boldsymbol{\psi}) = f(\boldsymbol{y}_{obs}|\boldsymbol{\theta}) f(\boldsymbol{r}|\boldsymbol{y}_{obs}, \boldsymbol{\theta}, \boldsymbol{\psi})$

である。実際の解析の場面ではN人分の尤度を用いることと（通常は異なる個人のデータは独立と仮定するので、単純にN人の尤度を掛け算する）、さらに対数をとったもの（対数尤度）を最大化するので、N人分の尤度について両辺の対数をとると

$\displaystyle \sum_{i=1}^{N} \log f(\boldsymbol{y}_{i, obs}, \boldsymbol{r}_i|\boldsymbol{\theta}, \boldsymbol{\psi}) = \sum_{i=1}^{N} \log f(\boldsymbol{y}_{i, obs}|\boldsymbol{\theta}) + \sum_{i=1}^{N} \log f(\boldsymbol{r}_i|\boldsymbol{y}_{i, obs}, \boldsymbol{\theta}, \boldsymbol{\psi})$

となる。

このように、欠測メカニズムがMARであれば完全尤度を「経時測定データに関する部分（第1項）」と「欠測指標に関する部分（第2項)」に分解することができる。経時測定データに関する部分であれば、確率分布が多変量正規分布であることを想定しているのでパラメータの推定は比較的簡単に行えそうである。

完全尤度の一部を用いれば適切な推測が可能である

長くなってきたので結論だけを述べると、完全尤度のうちの「経時測定データに関する部分」だけを用いて最尤法を行うことによって、パラメータ $\boldsymbol{\theta}$ 、つまり混合モデルであれば固定効果パラメータ・誤差の分散パラメータ・ランダム効果の分散パラメータは適切に推定できる。より正確には、

推定結果は $\boldsymbol{\theta}$ の一致推定量であり、
推定結果の従う確率分布はサンプルサイズが大きくなれば正規分布に近づく

ということが知られている。

駆け足で紹介してきたが、このような数学的背景が、欠測の可能性がつきまとう経時測定データ解析で混合モデルが広く用いられている一つの理由であるかと思う。
このように、欠測メカニズムがMARである限り、混合モデルなどの最尤法による推定を行う解析方法で関心のあるパラメータを適切に推定できることが保証されるが、そもそもMARと想定することが妥当であるかには注意が必要である。
式(2)のMARの定義について*9、左辺は

$\displaystyle f(\boldsymbol{r}|\boldsymbol{y}) = f(\boldsymbol{r}|\boldsymbol{y}_{obs}, \boldsymbol{y}_{mis}) = \frac{f(\boldsymbol{y}_{mis}|\boldsymbol{r}, \boldsymbol{y}_{obs}) f(\boldsymbol{r}|\boldsymbol{y}_{obs})}{f(\boldsymbol{y}_{mis}|\boldsymbol{y}_{obs})}$

のように変形できるので、両辺を $f(\boldsymbol{r}|\boldsymbol{y}_{obs})$ で割ることにより

$\displaystyle f(\boldsymbol{y}_{mis}|\boldsymbol{r}, \boldsymbol{y}_{obs}) = f(\boldsymbol{y}_{mis}|\boldsymbol{y}_{obs})$

となる*10。
つまり、 $\boldsymbol{y}_{obs}$ で条件付けすると、欠測が生じている変数の確率分布は欠測パターンに依存しない。
くだけた言い方をすると、欠測メカニズムがMARであると想定することは、「データに欠測がある個人の欠測した値は、データに欠測のない他のよく似た誰かとだいたい（確率的には）同じだろう」と想定していることになる。この想定が正しいかどうかは、そもそも欠測している値自体を見ることができないので、検証のしようがない。
f:id:mstour:20210303210111j:plain
この点が混合モデルをはじめとするMARを想定した解析の難点の一つだと言えるだろう*11。

まとめ

人に関するデータを経時的に測定する場合、さまざまな理由でデータの測定が途中でできなくなる可能性があるため、データの欠測という問題が避けられない。このような場合の解析方法として、混合モデルが用いられることが多くある。
今回は、混合モデルが欠測のある経時測定データに対して用いられている理論的な背景を簡単に追ってみた。一つの記事では全てを書ききれないので、関心をお持ちの方はしっかりとした書籍や論文を参考にしていただければと思うが、関連するトピックはまた随時書いていきたい。

本記事の作成に際しては、以下の資料を参考にした。
[1] 高井啓二・星野崇宏・野間久史(2016), "欠測データの統計科学", 岩波書店.
[2] 日本製薬工業協会(2016), "欠測のある連続量経時データに対する統計手法について Ver2.0".
http://www.jpma.or.jp/medicine/shinyaku/tiken/allotment/statistics.html

*1:単なる転居によって連絡がつかなくなることもあるし、研究や調査に参加するのが嫌になって拒否されてしまう場合もある。

*2:後で述べるように、厳密にはサンプルサイズが大きくなれば真の値に近づくということが保証されているということだが、統計学的にはこれで十分と考えることが多いので「正しく推定できる」と表現した。ただしサンプルサイズが小さい場合には注意が必要である。

*3:繰り返しになるが、最尤法によるパラメータ推定結果はサンプルサイズが大きくなれば真の値に収束する、というのがより正確である。このような性質を持つ推定結果は「一致推定量」と呼ばれる。

*4:MARは変数の集合に対して定義される概念であるため、例えば5つの変数でMARが成立していたとしても、そこから1つの変数を除くとMARではなくなる可能性がある。

*5:今後紹介する予定だが、正規分布を仮定できない経時測定データに対する解析方法の1つに一般化推定方程式(Generalized estimating equation; GEE)というものがある。しかし、これは最尤法を用いた推定を行わないため、混合モデルとは異なりMARの場合の妥当性は保証されない。

*6:表記を見やすくするために、個人を示すインデックス $i$ はしばらく用いないものとする。

*7:厳密には「確率密度関数」だが、直感的なわかりやすさのためにこのように表現することにする。

*8:欠測指標に関連するパラメータは $\boldsymbol{\psi}$ だけではなく $\boldsymbol{\theta}$ の一部も含む可能性があるので、「欠測指標に関するパラメータ」ではなく「それ以外のパラメータ」という表現とした。なお欠測メカニズムがMARであり、かつ欠測指標に関するパラメータが $\boldsymbol{\psi}$ のみである（つまりパラメータが「経時測定データだけに関するもの」と「欠測指標だけに関するもの」に完全に分けられる）場合を「無視可能(Ignorable)な欠測メカニズム」という。

*9:煩雑なのでパラメータの記載は省略する。

*10:説明変数 $\boldsymbol{x}$ をモデルに含める場合も同様である（条件部に入れればよい）

*11:「感度分析(Sensitivity analysis)」と称して、欠測メカニズムがMARでない（MNAR）というシナリオのもとでも結論が大きく変わらないかを検討するような事例が増えている。