ここんところ仕事に家事に疲弊していて、残り時間はひたすらぼーっと過ごしていた。時間蠅は矢を好むということわざの通りである。(← 疲れているとろくなことを書かない)

Merkle, E.C., Fitzsimmons, E., Uanhoro, J. Goodrich, B. (2021) Efficient Bayesian Structural Equation Modeling in Stan. J. Statistical Software, 100, 1-22.

　Rのblavaanパッケージ、というのはつまりはSEMのための定番パッケージのひとつlavaanのベイズ版なんだけど、その紹介。実戦投入する前の儀式として読んだ。本論文の前にMerkle & Rosseel (2018, 同誌)というのがあるんだけど、そっちはめんどくさいのでパス。
　私はMplusの忠実な信者なので、lavaan/blavaanなど無視してMplusのみと戯れるシンプルライフを送りたいのだが、なかなかそうも云ってられない事情がある。

1. イントロダクション
　ベイジアン潜在変数モデルの人気は高まっているが、JAGSやStanは難しいしサンプリングも時間がかかる。そこで我々はblavaanをつくった。lavaanのシンタクスが使えます。ほかにbrmsパッケージ(混合・多変量モデル)、rstanarmパッケージ(回帰モデル)、ctsemパッケージ(時系列モデル)、edstanパッケージ(IRTモデル)、pcFactorStanモデル(一対比較因子モデル)、とかがある。
　さて、blavaanはJAGSコードを実行時に生成し、runjagsパッケージ経由でJAGSで推定するパッケージであった。その後拡張しまして、rstan経由でStanを呼ぶこともできるようにしたんだけど、JAGSより速いわけでも効率的なわけでもなかった。
　このたびStanの呼び方を改善しまして、速度も効率も良くなりましたのでご紹介します。

2. モデルの定義
　lavaanのSEMモデル表現はLISRELに近い。
　SEMモデルを以下とする。$$ \mathbf{y}_i = \mathbf{\nu} + \mathbf{\Lambda} \mathbf{\eta}_i + \mathbf{\varepsilon}_i, \ \ \mathbf{\varepsilon}_i \sim N_p (\mathbf{0}, \mathbf{\Theta})$$ $$ \mathbf{\eta}_i = \mathbf{\alpha} + \mathbf{B} \mathbf{\eta}_i + \mathbf{\zeta}_i, \ \ \mathbf{\zeta}_i \sim N_m (\mathbf{0}, \mathbf{\Psi}) $$ \(\mathbf{y}_i\)は観察 \(i\) の\(p\)個の連続観察変数。\(\mathbf{\eta}_i \) は\(m \times 1\)の潜在変数ベクトル。残差の共分散行列は対角とすることが多い。そのとき \( \mathbf{y}\)の周辺分布はMVNとなり、その平均と共分散は $$ \mathbf{\mu} = \mathbf{\nu} + \mathbf{\Lambda} \mathbf{\alpha} $$ $$ \mathbf{\Sigma} = \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\Psi} (\mathbf{I} – \mathbf{B}^\top)^{-1} \mathbf{\Lambda}^\top + \mathbf{\Theta} $$ となる(\((\mathbf{I} – \mathbf{B})\)は逆行列を持つとする)。
　伝統的なLISRELフレームワークには外生観察変数の行列もあるが、lavaanでは「分散をすべて説明できる」潜在変数として扱う。
　SEMのベイズ推定では、\(\eta_i\)もほかのパラメータと一緒にサンプリングすることが多い。観察変数はふつう\(\eta_i\)の下で条件付き独立なので単変量正規分布になるというメリットがある。いっぽうMCMCがすごく遅くなり効率も悪くなる。

3. MCMCアプローチ
3つのアプローチを説明しよう。

3.1 JAGSにおけるパラメータ拡張
　関心あるモデル(推論モデル)を、サンプリングしやすい等価なモデル(作業モデル)に変換する。モデルの共分散パラメータに注目し、それぞれの共分散を「ファントム」潜在変数に変換するのである。こうすることで観察変数は互いに条件付き独立になり、JAGSでのサンプリングが速くなる。詳しくはMarkle & Rosseel (2018 J.Stat.Soft.) をみよ。
　blavaanでは、target = “jags” と指定するとこの方法となる。以前はこれがデフォルトだった。

3.2 Stanにおける尤度単純化
　Stanの場合、ファントム潜在変数を作ってもメリットがない。いっぽう、潜在変数の共分散行列 \( \mathbf{\Psi} \) のcapitalizationというのをやる。
　伝統的なSEMでは、潜在変数の分布は$$ \mathbf{\eta} \sim N( (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\alpha}, (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\Psi} (\mathbf{I} – \mathbf{B}^\top)^{-1} ) $$ となる。逆行列と行列式が出てくるので大変である。
　しかし、多くのモデルではもっと単純化できる。たとえば、\(\mathbf{B}\)は対角がゼロの三角行列、\(\mathbf{\Psi}\)は対角行列であることが多い。このとき… [中略] … というように単純化できる。
　blavaanでは、target = “stanclassic” と指定するとこの方法となる。

3.3 Stanにおける新しいアプローチ
　上記のどちらの手法も、モデルの潜在変数を他のモデルパラメータと一緒にサンプリングする。そうではなくて、潜在変数を周辺化したモデル尤度だけを扱えばいいじゃん、というのがHecht, et al.(2020 SEM)の提案である。潜在変数の事後分布はStanのgenerated quantitiesでサンプリングするのである。
　新しいアプローチでは、blavaanパッケージをインストールする時点で、あるStanプログラムをコンパイルしてしまう。そのプログラムは、ユーザが指定するMVNのSEMはたいていこのプログラムで推定できちゃう、というプログラムである。ちなみに、もしやる気があれば、このStanモデルにご自分で情報をパスしてキックしてもよい。
　このやり方の欠点はふたつある。第一に、多くのパラメータの事前分布は変えられない。\(B\)の事前分布はデフォルトでは\(N(0, 10)\)で、平均やSDを変えることしかできない。第二に、たとえば負荷どうしの等値制約とか切片同士の等値制約は掛けられるが、クラスの異なるパラメータ同士の制約はかけられない。またあるパラメータを他のパラメータの関数にすることもできない。そういうことがやりたいときは、我々の作ったStanファイルをexportして修正するか、上記の古い手法に戻ってください。
　blavaanでは、target = “stan” と指定するとこの方法となる。現在はこれがデフォルト。

4. 課題
　上記の手法は、直交因子のCFAというようなバニラ・モデルに適用するのは簡単だが、次のような難題がある。

• \(\mathbf{\Theta}, \mathbf{\Psi}\)の一部を固定したいときどうするか。blavaanでStanを呼ぶ場合、共分散行列はたとえば\( \mathbf{\Theta} = \mathbf{DRD} \) という風に分解される[原文では\(\mathbf{D}, \mathbf{R}\)に下添字\(\mathbf{\Theta}\)がついているが省略する]。\(\mathbf{D}\)はSDの対角行列で\(\mathbf{R}\)は相関行列である。で、それぞれに独立した事前分布を与える。そうするとときどき非正則な共分散行列がサンプリングされてしまうが、Stanはサンプリングを続けることができる(JAGSだと止まるけど)。他にも事前分布の案はいろいろある。今後の課題である。
• データにおける欠損をどうするか。blavaanは欠損データをMARと考える。JAGSでは勝手にサンプリングしてくれる。Stanでは完全情報尤度を用いる。
• SEMではふつう、識別のために個々の潜在変数のスケールを設定しなければならない。分散を1にするとか、どこかの負荷を1にするとか。ベイジアンでは前者のほうが難しい。負荷の符号が反転してしまうからだ。blavaanでは、JAGSの場合はある負荷について0でtruncateした事前分布を使う。Stanの場合はgenerated quantitiesブロックで、ある負荷が必ず正になるようにサンプリングのたびに反転をかける。

5. 適用例
　[実行例1。指標11個、因子3つのSEMで、因子をまたいだ残差相関がある。blavaan::bsem()で推定できる。メモ省略]

5.1 性能比較
　実行例1の場合、実行速度はJAGS, 新Stanが速く、旧Stanががくっと遅い。しかし秒あたりESSでみると新Stanがダントツである。
　実行例2 [3因子のCFA。blavaan::bcfa()で推定できる。メモ省略]。ここでも秒あたりESSで新Stanが圧勝。
　実行例3 [ちょっとめんどくさそうな成長モデル。blavaan::blavaan()で推定]。秒あたりESSで新Stanが圧勝。

5.2 検証とシミュレーション・ベースのカリブレーション
　lavaanを正とすると、上記3つのどのMCMC法でも、事後平均はML推定とほぼ同じ、事後SDはML推定よりちょっと大きめである。
　そこで新Stanについて、事後分布そのものについて検討してみよう。事前分布からドローして、たくさんのデータセットを生成する[つまり、パラメータの事前分布から500回ドローしたら500個のデータセットができるということだろう]。それぞれのデータセットにモデル当てはめ、MCMCの事後サンプルの、事前分布からのサンプルに対する順位を調べる。もしMCMCアルゴリズムがカリブレートされていたら、順位は一様分布になるはずである。
　[ええええ？ 頭がわるくてわからないよ… 事前分布から500回パラメータベクトルをドローして、それぞれについてデータセットをひとつ生成して、各データセットに対してモデルを当てはめ、サイズ1000の事後MCMCサンプルを得たとするじゃないですか。500×1000個の事後サンプルについて、事前分布でのパーセンタイル点を求めるということ？ それはなんだか変だよね… 理屈がよくわからない。Talts, et al.(2018, arXiv)というのがreferされている]

　結果。
　事前分布としてblavaanのデフォルトの弱情報事前分布を用いた場合、残念ながら、順位は0と1に集中する。これは、共分散行列が非正則になっちゃうことが多いからである。事後分布は、事前分布のうち共分散行列が正則になる領域にカリブレートされているというべきである。
　事前分布としてもう少し情報的な事前分布を与えると、順位は一様になる。
　教訓。研究者は情報的な事前分布を与え、MCMCサンプリング中に共分散行列が正則であり続けるようにしないといけない。いっぽう、無情報事前分布に含まれる情報の程度は、モデルの尤度として周辺尤度を使うのか条件付き尤度を使うのかによって変わってくる。
　この結果は、MCMCアルゴリズムについてのより詳細な研究におけるblavaanの価値を示している。さらに、blavaanのデフォルトの事前分布に改善の余地があることを示している。[華麗に逃げやがった…]

6. 結論
　ベイジアンSEMの推定の際に潜在変数を積分消去することでサンプリングの効率を改善できることを示した。この方法だと、標本サイズがパラメータ空間の次元に影響しなくなる。
　この周辺尤度アプローチは有用な反面、頻度主義者がSEMで直面していた問題が再浮上する。順序変数のような非正規の観察変数について周辺尤度が閉形式にならないとか、潜在変数の相互作用が扱えないとか。前者についてはdata augmentationで解決できるだろう。それに、DICやWAICを求めるためにはどのみち周辺尤度が必要なわけで、非正規なモデリングであっても、周辺尤度アプローチで進めたほうが有望だろう。
　[…後略…]
———————–
　よくわからんところもあったが、勉強になりました。
　Stanで因子分析とかやるとき、潜在変数をいちいちサンプリングしないといけないものなの？ Mplusだと因子得点はパラメータじゃなくて、あとでモデルパラメータの事後分布に基づいて生成する量なんだけど？ と思ってたんだけど、Stanであってもやっぱりそういう発想があるのね。さらに、そういう手法にはそれはそれで難しい点があるわけね。ふーん。