星野本「調査観察データの統計科学」はかなり難しい本だと思います。
第一章はいいのですが、第二章からいきなりつまづいてしまいました。
ここでは、星野本の最初のつまづきポイントである、選択モデルに対して理解できるように読み解いていきたいと思います。
最初に、前提知識として下記 PDF を読んでください。
この知識があるのとないのでは、理解度が全然違うと思います。
本記事では、星野本で省略された式変形について詳しく書いていますが、逆に星野本に書いてあることについては省略しています。
星野本の補足資料としてご利用ください。
それでは、星野本 p.28 の選択モデルのところから始めましょう。
関心のある変数を \(y\) とし、\(y = (y_{obs}^t, y_{miss}^t)^t\) を完全データのベクトルとし、\(y_{obs}\) を観測されているデータのベクトル、\(y_{miss}\) を欠測しているデータのベクトルとする。
ここで関心があるのは、完全データベクトルについての同時分布 \(p(y|\theta)\) である。
しかし、欠測データ \(y_{miss}\) が無いため、直接はパラメータ推定できない。
そこで、遠回りに思えるかもしれないが、欠測インディケータ \(m\) を用意して、\(y\) と \(m\) の同時分布を考える。
\[ p(y,m | \theta, \phi) \]
これは欠測が起こる仕組みを組み込んだモデルになっている。ここで、\(\phi\) は欠測メカニズムに関するパラメータである。
このモデルは、確率の乗法定理により次のように表される。
\[ p(y,m | \theta, \phi) = p(y|\theta,\phi)p(m|y,\theta,\phi) \]
このとき、\(y\) は \(\phi\) に無関係であり、\(m\) は \(\theta\) に無関係であるという仮定を置くと、次のようになる。
\[ p(y,m | \theta, \phi) = p(y|\theta)p(m|y,\phi)\]
これを、選択モデルという。
このままではまだ欠測データを含んでいるため、\(\theta\) を推定できない。
しかし、欠測が起こる仕組みをモデルに組み込んだため、次のようにできる。
\[ p(y_{obs},m|\theta,\phi) = \int p(y|\theta)p(m|y,\phi) dy_{miss} \]
欠測の仕組みをもとに欠測値を積分消去している。 これを完全尤度(フル尤度)という。
これは上記 PDF で言及されている FIML に使われる尤度関数のことである(ただし、上記 PDF では MAR の場合だけを論じているが、この式はより一般的である)。
この数式だけではわかりにくいので、具体的な例を提示する。
今、サンプルサイズ \(N\) 回のコイン投げを行う。
コインの表が出たとき \(y=1\)、裏が出たとき \(y=0\) を観測する。
コインの表が出る確率は \(\theta\) である。
ただし、\(1-\phi\) の確率で \(y\) は欠測するとする。
このとき、
\[ p(y|\theta) = \theta^y (1 - \theta)^{1-y}\]
\[ p(m|y,\phi) = \phi^m (1-\phi)^{1-m} \]
であるため、完全尤度は
\[ \begin{eqnarray} p(y_{obs},m|\theta,\phi) &=& \int p(y|\theta)p(m|y,\phi) dy_{miss}\\ &=& \int p(y_{obs}, y_{miss}|\theta)p(m|y_{obs}, y_{miss},\phi) dy_{miss} \\ &=& \int p(y_{obs}|\theta)p(y_{miss}|\theta) p(m=1|y_{obs},\phi)p(m=0|y_{miss},\phi) dy_{miss}\\ &=& p(y_{obs}|\theta) p(m=1|y_{obs},\phi) \int p(y_{miss}|\theta) p(m=0|y_{miss},\phi) dy_{miss}\\ &=& \prod_{i:m_i=1} \theta^{y_i} (1 - \theta)^{1-y_i} \phi^{m_i} \prod_{i:m_i=0} \sum_{y_{miss}\in\{0,1\}} \theta^{y_{miss}} (1 - \theta)^{1-y_{miss}} (1-\phi)^{1-m_i} \\ &=& \prod_{i:m_i=1} \theta^{y_i} (1 - \theta)^{1-y_i} \phi \prod_{i:m_i=0} \{ \theta (1-\phi) + (1 - \theta) (1-\phi)\} \\ &=& \prod_{i:m_i=1} \theta^{y_i} (1 - \theta)^{1-y_i} \phi \prod_{i:m_i=0} (1-\phi) \end{eqnarray} \]
となる。(\(y\) が離散のため積分は総和になる)
ランダムな欠測(MAR : missing at random)である場合を考えよう。
このとき、MAR の定義より、欠測するかどうかは欠測値には依存せずに観測値にのみ依存するので、
\[ p(m|y,\phi) = p(m|y_{obs},\phi)\]
が成り立つ。
したがって、完全尤度は
\[ \begin{eqnarray} p(y_{obs},m|\theta,\phi) &=& \int p(y|\theta)p(m|y,\phi) dy_{miss} \\ &=& \int p(y|\theta)p(m|y_{obs},\phi) dy_{miss} \\ &=& \int p(y|\theta) dy_{miss} \times p(m|y_{obs},\phi) \\ &=& \int p(y_{obs}, y_{miss}|\theta) dy_{miss} \times p(m|y_{obs},\phi) \\ &=& \int p(y_{obs}|\theta) p(y_{miss} | y_{obs}, \theta) dy_{miss} \times p(m|y_{obs},\phi) \\ &=& p(y_{obs}|\theta) \times \int p(y_{miss} | y_{obs}, \theta) dy_{miss} \times p(m|y_{obs},\phi) \\ &=& p(y_{obs}|\theta) \times 1 \times p(m|y_{obs},\phi) \\ &=& p(y_{obs}|\theta) p(m|y_{obs},\phi) \end{eqnarray} \]
となる。
これにより、対数完全尤度は
\[ \log p(y_{obs},m|\theta, \phi) = \log p(y_{obs}|\theta) + \log p(m|y_{obs},\phi) \]
と表される。
ここで、普通は \(\theta\) のみに関心があるため、\(\theta\) の最尤推定量を得るためには、\(\log p(y_{obs}|\theta)\) を最大化する \(\theta\) を求めればよいことがわかる。
このとき、\(p(y_{obs}|\theta)\) を「観測データの尤度」と呼ぶ。
入試得点 \(y_1\) と入学後の成績 \(y_2\) の関係を調べたい。
このとき、入学試験で合格点に達しなかった学生は入学できないためにデータの欠測が生じる。
具体的には \(y_1\) が \(C\) 以上の場合には合格して \(C\) 未満の場合には不合格なため \(y_2\) は欠測する。
この場合、欠測するかどうかは常に観測される \(y_1\) にのみ依存し、欠測している値 \(y_2\) には依存しないため、ランダムな欠測(MAR)である。
入試得点 \(y_1\) によって入学後の成績 \(y_2\) を予測する回帰モデル
\[ y_2 = \theta_1 + \theta_2 y_1 + \varepsilon \] \[ \varepsilon \sim Normal(0, \sigma) \]
を考えた場合、パラメータ \(\theta_1\)、\(\theta_2\) を推定するには、「観測データの尤度」
\[ p(y_{obs}|\theta) \]
を最大化すればよい。 観測データの尤度は、欠測インジケータ \(m\)(\(m=1\) のとき \(y_2\) を観測、\(m=0\) のとき \(y_2\) は欠測)を導入して、
\[ p(y_{obs}|\theta) = \prod_{i:m_i=1} p(y_1,y_2| \theta_1, \theta_2) \prod_{i:m_i=0} \int p(y_1,y_2|\theta_1,\theta_2) dy_2 \]
と表される。
また、
\[ \begin{eqnarray} \int p(y_1,y_2|\theta_1,\theta_2) dy_2 &=& \int Normal(y_2 | \mu=\theta_1 + \theta_2 y_1, \sigma) dy_2\\ &=& 1 \end{eqnarray} \]
であるため、結局
\[ p(y_{obs}|\theta) = \prod_{i:m_i=1} p(y_1,y_2| \theta_1, \theta_2) \]
となり、観測データの尤度を最大化するパラメータは、欠測データを無視して欠測していないデータのみで求めた値と一致する。