Asparouhov, T., Muthen, B. (2010) Bayesian Analysis Using Mplus: Technical Implementation. Mplus Technical Appendices.
Mplusのベイズ推定についての技術資料。共分散行列の事前分布を指定する方法について考えていたらわけがわからなくなり、遡って読んでみた。いささか古いテクニカル・ペーパーなのでこれまで無視していたんだけど…
膨大な内容のうちいま必要な部分に目を通しただけなので、実のところ読了とは言いがたいが、たぶん死ぬまで読み終えることはないだろうから、心の整理の都合上、読了ということにしておく。読了とはかくも柔軟な概念なのだ。
1. イントロダクション
Mplusにおける潜在変数モデルのベイズ推定の実装について詳細を述べる。アルゴリズムはGibbsサンプラーである。
Gibbsサンプラーというのは、パラメータとか潜在変数とか欠損値なんかの系列を反復的に生成して事後分布を構築する方法である。Gibbsサンプラーは、パラメータ、潜在変数、欠損値をいくつかのブロックにまとめて系列的に更新する。イテレーション\(t\)で生成されるすべてのパラメータ、潜在変数、欠損値からなるベクトルを\(\theta_t\)としよう。これを\( \theta_t = (\theta_{1t}, \ldots, \theta_{dt})\)と\(d\)個のブロックにわけ、$$ [\theta_{1t} | \theta_{2,t-1}, \ldots, \theta_{d,t-1}, data, priors] $$ $$ [\theta_{2t} | \theta_{1t}, \theta_{3,t-1}, \ldots, \theta_{d, t-1}, data priors]$$ $$ … $$ $$ [\theta_{dt} | \theta_{1t}, \ldots, \theta_{d-1, t-1}, data, priors]$$というように順繰りに生成するわけである。無事収束した暁には、Mplusは\(\theta_n, \ldots, \theta_m\)を使って事後分布を構築する。
収束は\(\theta_t\)のブロック化に強く依存する。Mplusは個々のモデルについてもっとも最適な分割を選び、状況に応じて異なる更新アルゴリズムを使う。おおまかにいえば、\(\theta_t\)において高い相関を持つ要素は同じグループにしておかないといけない。
MCMC系列の収束というのは、残念ながら簡単に診断できるものではない。[…中略…]
2. 連続変数のSEMモデル
2.1 モデル
\(p\)個の観察従属変数のベクトルを\(y\)とする。\(m\)個の潜在変数のベクトルを\(\eta\)とする。\(q\)個の観察変数のベクトルを\(x\)とする。以下の構造方程式モデルを考える。$$ y = \nu + \Lambda \eta + K x + \varepsilon $$ $$ \eta = \alpha + B \eta + \Gamma x + \zeta$$ \(\epsilon\)と\(\zeta\)は平均0の正規残差で、分散共分散行列は\(\Theta\), \(\Psi\)とする。いずれもブロック対角とする。なお、ブロック対角でなくても良いが、そのときはALGORITHM = GIBBS(RW); とする。するとランダム・ウォーク・アルゴリズムを使った推定が行われる。このアルゴリズムはどんな構造の分散共分散行列でも推定できるが、Mplusのデフォルトのアルゴリズムに比べるとmixingの質が悪く、ちょっと効率が低い。
2.2 モデル拡張
さっきの\(y\)の式の右辺に架空の潜在変数を入れることでより柔軟にできる。たとえば、\(Y_1\)を\(Y_2\)に回帰したいなら、架空の潜在変数\(\eta_1, \eta_2\)を考え、\(Y_1 = \eta_1, Y_2 = \eta_2\)として\(\eta_1 = \alpha + \beta \eta_2 + \ldots \)とすればよい。Mplusはなにもいわなくてもこの3本の式を自動で実装してくれる。
2.3 事前分布
\(\gamma\)をすべての自由パラメータ\(\nu, \alpha, \Lambda, B, \Gamma, K\)のベクトルとする。事前分布として正規分布 $$ \gamma \sim N(\gamma_0, \Omega_\gamma) $$ を仮定する。さらに、ブロック対角な共分散行列\(\Psi, \Theta\)のすべてのブロックについて、事前分布として逆ウィシャート分布を仮定する。たとえば$$ \Psi_{jj} \sim IW(\Omega_{\Psi_{jj}}, d_{\Psi_{jj}}) $$ というように。ただし、ブロックのサイズが1であるときは逆ガンマ分布 $$ \Psi_{jj} \sim IG(\alpha_{\Psi_{jj}}, \beta_{\Psi_{jj}})$$ となる。
以上のすべての事前分布は共役事前分布である。つまり、Gibbsサンプラーが生成する条件付き分布は事前分布と同じ族の分布である。
2.4 推定
パラメータは3つのブロックに分割される。すなわち、(1)切片と負荷 \(\gamma\), (2)分散共分散 \(\Psi, \Theta\), (3)潜在変数\(\eta\)である。
推定は3つのステップに分かれる。
- ステップ1、\(\eta\)の更新。さっきの式$$ \eta = \alpha + B \eta + \Gamma x + \zeta$$ を次のように書き換える: $$ \eta = B^{-1}_0 (\alpha + \Gamma x) + \zeta_0 $$ ただし\(\zeta_0 = B^{-1}_0 \zeta\)で、その共分散行列は\(\Psi_0 = (B^{-1}_0) \Psi (B^{-1}_0)^\top \)となる。\(\eta\)の条件付き分布は下式で与えられる。 $$ [\eta | *] \sim N(Dd, D) $$ $$ D = (\Lambda^\top \theta^{-1} \Lambda L \Psi^{-1}_0) $$ $$ d = \Lambda^\top \theta^{-1} (y – \nu – Kx) + \Psi^{-1}_0 B^{-1}_0 (\alpha + \Gamma x) $$
- ステップ2、\(\gamma\)の更新。さっきの式 $$ y = \nu + \Lambda \eta + K x + \varepsilon $$ $$ \eta = \alpha + B \eta + \Gamma x + \zeta$$ をこう書き換える。$$ z = F \gamma + \epsilon $$ ここで \(\epsilon = (\varepsilon, \zeta), z = (y, \eta)\)である。\(F\)は\(p+m\)行\((p+m)(1+m+q)\)の行列で$$ F = I \otimes (1, \eta, x)$$ である。[ここ、前にもぽかーんとしたんだよな… まあとにかく、\(y\)と\(\eta\)を縦に積んだベクトルが左辺にある一次式に書き換えられるのね、はい、信じます]
\(\gamma\)の条件付き分布は下式で与えられる。ケースのインデクスを\(i\)として、$$ [\gamma | * ] \sim N(Dd, D)$$ $$ D = (\sum_i^n F_i V^{-1} F_i + \Omega^{-1}_\gamma)^{-1} $$ $$ d = \sum_i^n F_i V^{-1} z_i + \Omega^{-1}_\gamma \gamma_0$$ ただし、\(V\)はブロック対角に\(\Theta, \Psi\)を持つブロック対角行列である。 - ステップ3, 分散共分散行列の生成。このステップについていえば、もはや\(\eta\)と\(Y\)に本質的な違いは無いので、\(\Theta\)の条件付き分布についてのみ述べる。
ブロックサイズが1のとき。事前分布は$$ \theta_{11} \sim IG(\alpha_0, \beta_0) $$ 条件付き分布は$$ [\theta_{11}|*] \sim IG(\alpha_0 + \frac{n}{2}, \beta_0+\beta_1)$$ $$ \beta_1 = \frac{1}{2} \sum_i^n (y_{i1} – \nu_1 – \sum_j^m \lambda_{1j} \eta_{ij} – \sum_j^q K_{1j} x_{ij})^2 $$
ブロックサイズが1より大きいとき。事前分布は $$ \Theta_{11} \sim IW(\Omega, f)$$ 条件付き分布は$$ [\Theta_{11}|* ] \sim IW(E + \Omega, n + f)$$ ただし\(E\)は、以下の行列の最初の対角ブロックである。$$ \sum_i^n (y_i – \nu – \Lambda \eta_i – K x_i)(y_i – \nu – \Lambda \eta_i – K x_i)^\top $$ [うーん、よくわからん。\(\Theta\)は観察方程式の誤差項の共分散行列だから、ふつうは対角行列だよね？ ってことは\(\theta_{11}, \ldots, \theta_{pp}\)のそれぞれについて上記をそれぞれ生成するってことだよね？ ブロックサイズが1より大きいのってどういうときなの？ あっそうか、\(\Psi\)も\(\Theta\)と同様に扱われるわけだから、そっちでは因子間に相関があるのは珍しくないな]
2.5 収束
[ここは関心あるんだけど、いまちょっと時間が無いのでパス]
[以下、章見出しのみメモする]
3. 連続変数とカテゴリカル変数のSEMモデル
4. 欠損データのあるSEMモデルの推定
5. パラメータの固定とパラメータ間の等価性
6. 非共役事前分布
7. 観察された媒介変数
8. 混合モデル
9. マルチレベルモデル
10. Posterior Predictive Checking
11. 多重代入
12. 情報事前分布と無情報事前分布
- まず正規事前分布 \(N(\mu, \upsilon)\)について。無情報事前分布なら\(N(0, 10^{10})\)がよいだろう。これは数値的にいえばもはや\((-\infty, \infty)\)にわたって密度1という非正則事前分布に等しい。
Mplusでは、正規変数の切片・負荷・傾きのデフォルト事前分布は\(N(0, 10^{10})\)である。カテゴリカル変数の切片・負荷・傾きのデフォルト事前分布は\(N(0, 5)\)である。閾値のデフォルト事前分布は\(N(0, 10^{10})\)である。情報事前分布を指定するのも簡単である。 - 逆ガンマ事前分布\(IG(\alpha, \beta)\)について。密度関数を定義しておくと、\(x \gt; 0\)について$$ f(x) \sim x^{-\alpha-1} \exp(-\beta/x) $$ である。最頻値は\(\beta / (\alpha+1)\)である。
無情報事前分布なら、\(IG(-1, 0)\)がよいだろう。これは\((-\infty, \infty)\)にわたって密度1である。\(IG(0, 0)\)とか\(IG(0.001, 0.001)\)というのも良い。小標本でランダム効果の分散が潰れちゃったり爆発しちゃったりしないようにするには、経験的に \(IG(1,2)\)がおすすめである。[← へー！] - 一様分布\(U(a,b)\)について。無情報事前分布なら\(U(-10^{10}, 10^{10})\)とか、分散パラメータなら\(U(0, 10^{10})\)とか、がいいんじゃないでしょうか。
- ガンマ分布\(G(\alpha, \beta)\)について。定義しておくと、密度関数は\(x > 0\)について $$ f(x) \sim x^{\alpha – 1} \exp(-\beta x)$$ である。
- 対数正規分布\(LN(\mu, \upsilon)\)について。密度関数は$$ f(x) \sim x^{-1} \exp(-(\log(x) – \mu)^2 / (2\upsilon)) $$ である。
- 逆ウィシャート分布\(IW(\Omega, d)\)について。\(\Omega\)はサイズ\(p\)の正則行列、\(d\)は整数である。密度関数は $$ f(\Sigma) \sim |\Sigma|^{-\frac{d+p+1}{2}} \exp(-\frac{Tr(\Omega \Sigma^{-1})}{2}) $$ である。最頻値は\(\frac{\Omega}{d + p + 1}\)となる。
無情報事前分布としては\(IW(0, -p-1)\)がよい。つまり、\(\Omega\)の全要素を0とし\(d = -p-1\)とする。これは本質的に、すべての\(\Sigma\)パラメータについて\((-\infty, \infty)\)上の一様事前分布を与えていることになる。
その他の無情報事前分布として、\(IW(0, 0)\)とか\(IW(I, p+1)\)もある。\(IW(I, p+1)\)というのは、相関パラメータの周辺事前分布が\((-1, 1)\)について一様になり、分散パラメータの周辺事前分布が\(IG(1, 0.05)\)となることを意味している。
\(\Sigma\)に対して情報事前分布\(IW(\Omega, d)\)を与える際には、\(i = 1, \ldots, p, j= 1, \ldots, i\)について \(\sigma_{ij} \sim IW(\Omega_{ij}, d)\)と指定すること。ただし\(\Omega_{ij}\)とは\(\Omega\)行列の\((i,j)\)要素である。
分散共分散行列のデフォルト事前分布は、連続変数だったら\(IW(0, -p-1)\), カテゴリカル変数が混じっていたら\(IW(I, p+1)\)である。 - 最後にディリクレ分布について。密度関数は、パラメータを\(\alpha_1, \ldots, \alpha_K\)として$$ f(x_1, x_2, \ldots, x_K) = x^{\alpha_1-1}_1 x^{\alpha_2-1}_2 \ldots x^{\alpha_K-1}_K$$ ただし\(x_1 + \ldots + x_K = 1\)である。Mplusでは混合モデルのクラス割合の事前分布に用いられる。
一般的な無情報事前分布は \(\alpha_i = 1\)なのだが、Mplusのデフォルトは\(\alpha_i = 10\)である(クラスが小さい解を作りたくないから)。もっとも、小標本の場合はあまり良いデフォルトではないかもしれない。
情報事前分布を与えるときには、パラメータ[C#i]にラベル\(p_i\)を与えて、\(i=1, \ldots, K-1\)について\(p_i \sim D(\alpha_1, \alpha_K)\)と指定すること。
—————-
これまでさんざんMplusを使ってきたにも関わらず、全然理解していなかったことばかりだ…