Lu, Z.H., Chow, S.M., Loken, E. (2016) Bayesian Factor Analysis as a Variable-Selection Problem: Alternative Priors and Consequences. Multivariate Behavioral Research, 51(4), 519–539.

　仕事の都合で因子分析モデルをベイズ推定するとき、いつも悩むのは因子負荷の事前分布の設定である。導師Muthenが提案するベイジアンSEMは確かにすごく有用な手法だと思うし、日本語での紹介が少ないことに義憤を感じて書籍で紹介させて頂いたりもしたんだけど(嗚呼、自己満足)、個別具体的な事例においては、どうしても困っちゃうわけです。いったい負荷の事前分布とはなんなのか… 我々はデータについて何を知っているのか… この世界のなりたちとは… 人生の意味とは… 眠い… 今日は寝よう… っていう風になります。
　
　というわけで、たまたまみつけた論文を読んでみた(目先の仕事からの現実逃避であるともいえる)。ベイジアン因子分析において、因子負荷にspike-and-slab事前分布を与えるのがよろしいのではないでしょうか、という論文。

[イントロダクション]
　次の因子分析モデルを考える。標本サイズを\(N\)として、\(i=1,\ldots,N\)について$$ \mathbf{y}_i = \mathbf{\mu} + \mathbf{\Lambda} \mathbf{\omega}_i + \mathbf{\epsilon}_i$$ $$\mathbf{\omega}_i \sim N_q(\mathbf{0}, \mathbf{\Phi})$$ $$\mathbf{\epsilon}_i \sim N_p(\mathbf{0}, \mathbf{\Psi})$$ ただし\(\mathbf{y}_i, \mathbf{\mu}, \mathbf{\epsilon}_i\)はいずれも長さ\(p\)のベクトル、\(\mathbf{\omega}_i\)は長さ\(q\)の潜在変数ベクトル。\(\mathbf{\Phi}\)は正定値行列, \(\mathbf{\Psi}\)は対角行列とする。
　EFAではなんらか推定してから回転して\(\mathbf{\Lambda}\)を倹約的にする。いっぽうCFAでは、はなから\(\mathbf{\Lambda}\)の一部を0に固定する。その分パラメータのSEは下がる。でも、十分な事前知識がないまま\(\mathbf{\Lambda}\)に0を埋めちゃうと\(\mathbf{\Phi}\)にバイアスが生じちゃいますよね。EFAの柔軟さとCFAの倹約性のいいとこどりをできないものか？
　ひとつのアイデアは修正指標 (MI) を使うこと。しかし、回帰におけるステップワイズ変数選択と同じで、MIでもってモデルを繰り返し修正していったとき、そこから生まれる分散を、標本抽出による分散からどうやって分離するのかという問題が生じる。
　そこでこう提案したい。因子負荷推定をベイジアン変数選択問題として捉えよう。導師MuthenらによるベイジアンSEMアプローチ(BSEM)はその特殊ケースと捉えられるのだ！
　[導師っていうのは私が勝手にそう呼んでるだけで、原文にはない、もちろん]

変数選択
　次の回帰モデルを考える。$$ y_i = \mathbf{x}_i \mathbf{\beta} + \epsilon_i$$ \(\mathbf{x}_i, \mathbf{\beta}\)は長さ\(p\)のベクトル。ふつう\(\epsilon_i \sim N(0, \sigma^2)\)と仮定する。簡略のため切片項はなしってことにする。
　\(\epsilon_i\)がiidならML推定量(兼OLS推定量)は $$ \hat{\mathbf{\beta}}_{MLE} = \mathop{\rm arg~min}\limits_\beta \left[ \sum_i^N (y_i – \mathbf{x}_i \mathbf{\beta})^2 \right] $$ ですわね。
　古典的な変数選択手続きはいろいろあるけど欠点もある。サブセット選択は\(p\)がちょっと大きくなるだけで非現実的になる。系列探索(前向き選択とか後ろ向き削除とかステップワイズとか)はデータセットごとに全然ちがう結果になっちゃうし、モデル選択の不確実性を考慮してないせいでSEの過小評価につながるし…[後略]

　そこで出てくるのが正則化。\(\mathbf{\beta}\)の非ゼロ要素数を決めるんじゃなくて、損失関数に\(\mathbf{\beta}\)のノルムをいれる。$$ \hat{\mathbf{\beta}}_R = \mathop{\rm arg~min}\limits_\beta \left[ \sum_i^N(y_i -\mathbf{x_i} \mathbf{\beta})^2 + \frac{norm(\mathbf{\beta})}{\tau^2} \right]$$ ノルムとしてよく使われているのは、L0ノルム \( || \mathbf{\beta} ||_0 \)=非ゼロ要素数、L1ノルム\( || \mathbf{\beta} ||_1 = \sum_j^p |\beta_j| \), L2ノルム \( || \mathbf{\beta} ||_2 = \sqrt{\sum_j^p \beta_j^2} \)である。\(norm(\mathbf{\beta})/\tau^2\)を罰則項という。
　\(\hat{\mathbf{\beta}}_R\)の有名な特殊ケースとしてridgeとlassoがある。[…それぞれについて簡単に説明している。めんどくさいので中略…]
　というわけで、頻度主義的な回帰の枠組みにおいて、正則化とは目的変数に罰則をつけることである。

　よく知られている正則化手法は、罰則関数を事前分布として指定したベイジアンモデルとしても捉えられる。ベイジアンの手法についてはMuthen & Asparouhov (2012 Psych.Method)をみてほしいんだけど… ベイジアン正則化は伝統的なベイジアン手法とは異なり、事前分布でもって未知パラメータの値の分布を指定するだけでなく、モデル構造にまで影響を与える。
　回帰モデルの場合、\(\mathbf{\beta}\)の事前分布としてはガウシアンかラプラスが選ばれることが多い。順に$$ p(\mathbf{\beta} | \sigma^2) \propto \exp \left( -\frac{1}{2\sigma^2} \sum_j^p \beta_j^2 \right)$$ $$ p(\mathbf{\beta} | \sigma^2) \propto \exp \left( -\frac{1}{2\sigma^2} \sum_j^p |\beta_j| \right)$$ これが頻度主義でいうridge推定量, lasso推定量と等しいことを示せる。

　[話の進め方が早すぎるので、ridge推定量の場合についてゆっくりやり直す。
　頻度主義側からみると、ridge推定量はL2ノルムを使って$$ \hat{\mathbf{\beta}}_R = \mathop{\rm arg~min}\limits_\beta \left[ \sum_i^N(y_i – \mathbf{x_i} \mathbf{\beta})^2 + \frac{1}{\tau^2} \sum_j^p \beta_j^2 \right]$$ である。
　いっぽうベイジアン側からは、$$ (y_i | \mathbf{\beta}, \mathbf{x}_i) \sim x_t^\top \mathbf{\beta} + e_i, \ \ \epsilon_i \ \sim_{iid} \ N(0, \sigma^2_e I)$$と仮定し、\(\mathbf{\beta}\)の個々の要素について\(\beta_j \sim_{iid} N(0, \sigma^2)\)と仮定する。正規分布の確率密度関数はどんなんでしたっけ？ たしか $$ f(x) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left( – \frac{(x – \mu)^2}{2 \sigma^2}\right) $$ でしたよね。つまり個々の要素について$$ p(\beta_j|\sigma^2) \propto \exp \left( -\frac{1}{2\sigma^2} \beta_j^2 \right)$$だ。すべての要素の同時分布は$$ p(\mathbf{\beta}|\sigma^2) \propto \exp \left( -\frac{1}{2\sigma^2} \sum_j^p \beta_j^2 \right)$$ である。これが上記の式である。
　頻度主義側からみてもベイジアン側からみても結局は同じことである。なぜか。ベイジアン側からみて、\(\mathbf{\beta}\)の事後分布は $$ p(\mathbf{\beta}|\mathbf{y}, \mathbf{x}) \propto p(\mathbf{\beta}) p(\mathbf{y} | \mathbf{x}, \mathbf{\beta})$$ である。再び正規分布の密度関数を思い出すと$$ p(\mathbf{y} | \mathbf{x}, \mathbf{\beta}) \propto \exp \left( -\frac{\sum_i^n (y_i – \mathbf{x}_i \mathbf{\beta})^2}{2\sigma^2_\epsilon} \right)$$ 正規分布と正規分布のかけ算は正規分布になるから、\(\mathbf{\beta}\)の事後分布は結局 $$ p(\mathbf{\beta}|\mathbf{y}, \mathbf{x}) \propto \exp \left( -\frac{\sum_i^N (y_i – \mathbf{x}_i \mathbf{\beta} )^2}{2\sigma^2_\epsilon} – \frac{\sum_j^p \beta_j^2}{2\sigma^2} \right)$$ \(\mathbf{\beta}\)の事後分布の最頻値を\(\hat{\mathbf{\beta}}_{Bayes}\)と呼ぶとして、それはこの密度関数を最大化する\(\mathbf{\beta}\)、つまり、\( \sum_i^N (y_i- \mathbf{x}_i \mathbf{\beta})^2 + \frac{\sigma^2_\epsilon}{\sigma^2} \sum_j^p \beta_j^2 \) を最小化する\(\mathbf{\beta}\)、つまり$$ \hat{\mathbf{\beta}}_{Bayes} = \mathop{\rm arg~min}\limits_\beta \left[ \sum_i^N(y_i -\mathbf{x}_i \mathbf{\beta})^2 + \frac{\sigma^2_\epsilon}{\sigma^2} \sum_j^p \beta_j^2 \right]$$ である。なるほど等しいですね]

　ベイジアン正則化には頻度主義にくらべて利点がある。頻度主義lassoではSEの算出が難しいがベイジアンなら事後分布からサンプリングすればよい。モデルの不確実性もモデルの事後確率として示せる。推定できるモデルの幅もより広い。

ベイジアンCFAにおける交差負荷の正則化
　では、因子分析の負荷の事前分布はどうするか。ふたつのアプローチを比べよう。

　その1，Muthen & Asparouhov(2012)のアプローチ。負荷行列のゼロのかわりに分散の小さい情報分布を埋める。
　CFAの場合、\(\Lambda\)の行\(j\), 列\(k\)を\(\lambda_{jk}\)として、$$ p(\Lambda | \Sigma_0) \propto \exp \left( \sum_j^p \sum_k^q – \frac{1}{2\sigma^2_{jk}} \lambda^2_{jk} \right) $$ ただし\(\Sigma_0\)は\(\sigma^2_{jk}\)の行列。
　標準的なCFAでは、\(\lambda_{jk}\)を主負荷と交差負荷に分けて考え、交差負荷をゼロにする。ここでは、主負荷の\(\sigma^2_{jk}\)を十分に大きくし、交差負荷の\(\sigma^2_{jk}\)を0のまわりですごく小さくする(0.01とか)。他の識別のための制約はかけない。
　[ちょっと戸惑ったんだけど、この論文では基本的に、ユーザがゼロだと思っている\(\lambda_{jk}\)が交差負荷、ゼロじゃないと思っている\(\lambda_{jk}\)が主負荷である。単純構造に合致するのが主負荷、単純構造からはみ出すのが交差負荷、という話ではない。うーん、cross loadingという言葉からは「この項目は因子Aの指標なんだけどBもちょっと反映しちゃっている」というような、ある項目のふたつめの負荷というようなニュアンスを感じるんだけど、そういう意味合いはないわけね]

　Muthen導師らはそうは言ってないけど、これはベイジアンridge回帰事前分布(RP)の特殊ケースである[あ、そうか、確かに…]。ちがいは、ベイジアンridgeでは事前分散もパラメータであってそのまた事前分布を持つけど、Muthen-Asparouhovでは固定だという点である。以下、このアプローチをBSEM-RPと呼ぼう。

　どうやって実装するのかというと… まず\(\mathbf{\Phi}\)[因子得点の共分散行列]、\(\mathbf{\Psi}\)[独自因子の共分散行列。対角行列]、\(\mathbf{\Lambda}\)[因子負荷行列] に事前分布を与える。$$ p(\mathbf{\Phi}) = IW(\rho_0, \mathbf{\Phi}_0)$$ $$ p(\Psi_j) = IG(\alpha_{1j}, \alpha_{2j})$$ $$ p(\mu_j | \Psi_j) = N(\mu_{0j}, \Phi_j \sigma^2_{\mu 0j}) $$ $$ p(\lambda_{jk} | \Phi_j) = N(\lambda_{0jk}, \Psi_j \sigma^2_{\lambda 0jk}) $$という段取りである。ハイパーパラメータは、\(\mathbf{\Phi}_0\)(正定値行列)、\(\rho_0 (> q-1) \)、\(\alpha_{1j} (> 0), \alpha_{2j} (> 0) \)、\(\mu_{0j}\), \(\sigma^2_{\mu 0j}(> 0)\), \(\sigma^2_{\lambda 0jk}(> 0)\)である。
　なおMplusのデフォルトでは、最初の2本が$$ p(\mathbf{\Phi}) \propto 1$$ $$ p(\Psi_j) \propto 1$$ である。非正則なわけです。
　でもって、Gibbsサンプリングで\(\mathbf{\Lambda}, \mathbf{\Phi}, \mathbf{\Psi}\)の事後分布を得る。[原文ではここに\(\mathbf{R}\)というのも入っているんだけど、これは誤記じゃないかしらん？ BSEM-RPには\(\mathbf{R}\)ってのは出てこないと思う]
　こうやって因子負荷とその信用区間を調べてから、CFAモデルを修正して再度推定せよ、というのが導師Muthenらの主張である。修正指標を使うやり方のような逐次的なモデル選択はやらないけど、二段階の手順ではある。

　その2, 我々の提案。spike-and-slab事前分布(SSP)を使う。交差負荷について、識別のために一部を0に固定し、残りについては$$ p(\lambda_{jk}|r_{jk}) = (1-r_{jk})\delta_0 + r_{jk} N(0, c^2_{jk})$$ $$ p(r_{jk} = Bernoulli(p_jk) $$ \(\delta_0\)は0におけるpoint mass function[0となる確率が1ですという確率密度関数のことだろう]。\(r_{jk}\)はなにをしているのだとお嘆きのみなさん、これは混合モデルで使うような二値潜在変数でありまして、そこの負荷がゼロなのかはたまた自由推定なのかを決めているのです。つまり、\(\lambda_{jk}=0\)ですというspikeと、\(\lambda\)はなんだかわかりませんというslabの混合である。\(p_{jk}\)はハイパーパラメータで、ユーザの事前知識を表現する。
　SSPは階層正規混合事前分布$$ p(\lambda_{jk}|r_{jk}) = (1-r_{jk}) N(0, \sigma^2_{jk}) + r_{jk} N(0, c^2_{jk})$$ $$ p(r_{jk} = Bernoulli(p_jk) $$ の特殊ケースである。[…中略…]
　以下、このアプローチをBSEM-SSPと呼ぼう。

　どうやって実装するのかというと… 事前分布の与え方はBSEM-RPと同じ[因子負荷以外は、ってことだろう]。\(r_
{jk}\)を要素として持つ\(p \times q\)の二値行列を\(\mathbf{R}\)とする。
　仮定されているモデルの交差負荷構造を\(\tilde{\mathbf{R}} = (\tilde{r}_{jk})\)とする。因子負荷の点推定は次の3種類考えられる。

• SSP確認的推定量: \(E(\lambda_{jk}|\mathbf{Y}, \mathbf{R} = \tilde{\mathbf{R}}) \)
• SSP条件付き推定量: \( E(\lambda_{jk} | \mathbf{Y}, r_{jk} = 1)\)
• SSP周辺推定量: \( E(\lambda_{jk} | \mathbf{Y}) \)

　CFAに対応するのは確認的推定量である。しかしMCMCで推定する際、MCMCサンプルのうち\(\mathbf{R}\)が\(\tilde{\mathbf{R}}\)である奴だけをとってきて平均したのが確認的推定量なんだけど、そのサンプルはすごく小さくなるだろう。むしろ条件付き推定量を使ったほうがいい[つまり、他の負荷の\(r_{jk}\)を無視して、今見ている負荷の\(r_{jk}\)だけを1と条件づけたときの\(\lambda_{jk}\)の期待値をみればいいじゃんってことね]。
　リサーチャーがモデルに確信を持っておらず、モデルに不確実性を組み込みたいときはSSP周辺推定量がよろしい。

　BSEM-SSPは一段階の手順である。もっとも、モデル比較指標をMCMCサンプルから求めることもできる。[ここからよくわからんので逐語訳]

たくさんの候補モデルから最終モデルを選択するためのひとつの方法は、\(　p(r_{jk} = 1 | \mathbf{Y})\ > 0.5 \)を持つ交差負荷だけを保持した最終モデルを選ぶ to select a final model in which only cross loadings with \(　p(r_{jk} = 1 | \mathbf{Y})\ > 0.5 \) are retainedという方法である。これはmedian probability model (Berger & Barbieri, 2004)となる。このモデルは、独立同分布のガウシアン誤差を持つ回帰の最適予測モデルである。

[わ、わからん… リサーチャーが因子負荷0だと思ってる場所(この論文でいう交差負荷)と自由推定したいと思ってる箇所(主負荷)がありますよね。BSEM-SSPではこの信念を\(p_{jk}\)の行列として与えるわけですよね。で、MCMCサンプルをたくさん得ますよね。二値行列\(\mathbf{R}\)のサンプルを平均しますよね。すると、交差負荷なのに1が5割超えている箇所や、主負荷なのに1が5割切っている箇所が出てきますよね。ここで述べているのは、\(p_{jk}\)をどう与えたかは忘れて、MCMCサンプルで1が5割超えている箇所は主負荷、そうでない箇所を交差負荷にした構造を仮定しろってこと？ それとも、主負荷だと思っていた場所は1が何割であろうが主負荷のままで、交差負荷のうち1が5割超えてるやつを交差負荷に格上げした構造を仮定しろってこと？
　それからもう一つ疑問なのは… そうやって新しいモデルを仮定したとして、そのパラメータ推定はどうするの？ 著者らの主張では、CFAをやりなおすわけじゃなくて、このMCMCサンプルから求めるわけですよね？ ひょっとして、主負荷の推定値としては \( E(\lambda_{jk} | \mathbf{Y}, r_{jk} = 1)\)を使え、交差負荷の推定値は0だと思え、ってこと？ でもそれおかしくないですか、新たに仮定したモデルを\(\tilde{\mathbf{R}}_{final} = (\tilde{r}_{final,jk})\)と書くとして、推定したいのは原理的には\( E(\lambda_{jk} | \mathbf{Y}, r_{jk} = \tilde{r}_{final,jk})\)ではなくて\(E(\mathbf{\Lambda}|\mathbf{Y}, \mathbf{R} = \tilde{\mathbf{R}}_{final}) \)でしょう？ MCMCサンプルからそういうサブセットをとってくるのは難しいというのはよくわかるけど、それはあなたの都合であって、難しいってんならCFAをやり直すべきなのでは？]

モデル比較は、すべての候補モデルの事後確率を評価することで行うこともできるだろう。いま関心ある仮説的な因子分析モデルを\(M\)とし、それに対応する交差負荷構造を\(\tilde{\mathbf{R}} = (\tilde{r}_{jk})\)とする。\(M\)も事後確率\(p(\mathbf{R} = \tilde{\mathbf{R}} | \mathbf{Y})\)を、BSEM-SSPをフィッティングして得たMCMCサンプルから計算できるし、ベイズ・ファクターも求めることができる。

[なるほど、このくだりはわかる。きっとここでは、「負荷行列のどこにゼロを埋めるか」というモデル選択の話だけをしてるんだろうな。でもモデルを選択したら当然パラメータを推定したくなるわけで、その推定量はどうするのか…やっぱりぴんととこないなあ]

　BSEM-SSPには次の理論的利点がある。

• BSEM-CFAのハイブリッドという枠組みでうまく変数選択できるかどうかは、ゼロ負荷と非ゼロ負荷を正しく区別できるかどうかにかかっている。BSEM-SSPはBSEM-RPやEFAと比べてより倹約的なモデルになりやすい。非ゼロ負荷については、SSPでは\(r_{ij} = 1\)の事後分布を学習するので検出力が高く、大きな負荷についても縮小が小さい。いっぽうBSEM-RPは大きな負荷に大きな罰則をかけるので負荷が大きいときに縮小しすぎてしまう。リッジ回帰と同様に、すべての交差負荷が比例的に0に縮小してしまうのである。
• BSEM-SSPは一段階アプローチであり、CFAをやりなおす必要がない。
• BSEM-SSPの事前分布は情報的でないので、事後分布への影響が小さい。いっぽう導師Muthenらは事前分布の分散を動かしてsensitivityをチェックせよといっている。[そうなのか…　それはまあ、ユーザとしては素直に魅力的だな]
• BSEM-SSPならベイズ・ファクターを簡単に計算できる。負荷が0である(ない)確率も定量評価できるしそれを推定値に組み込める。

シミュレーション
　さあシミュレーションだ。がんばっていきましょう！！ [←原文には書いてないけど、私には聞こえるぞ、世界の読者の声が]

　ふたつの因子負荷構造を考える。どちらも3因子, 9変数。[以下、原文から離れてMplusの文法でメモする。まず、どちらの因子負荷構造でも V1-V3 by F1; V4-V6 by F2; V7-V9 by F3; である(V1, V4, V7の負荷が1に固定されている点に注意)。以上を主負荷と呼ぶ。
　因子負荷構造その１、単一交差負荷条件。上記に V2 by F2 (x); を追加する。
　因子負荷構造その２、多重交差負荷条件。V4 by F1 (x); V7 by F2 (x); V1 by F3 (x); を追加する。なんと、xは3カ所で共通なのね…]

　シミュレーションで動かす要因は以下の通り。(1)因子負荷構造、上記の2水準。(2)標本サイズ。{100, 200, 300}。(3)交差負荷xの強さ。{0.1, 0.2, 0.3}。2x3x3=18条件となる。
　因子負荷以外のパラメータは、\(\mu_j = 0, \Psi_j = 0.3\), \(\Psi\)は対角1, 非対角0.5に固定。各条件あたりMCMCを4000回転バーンインしてから4000回転。

　選手入場！

• BSEM-SSP。すべての0とxにspile-and-slab事前分布を与えます。\(p_{jk} = 0.5\), \(c^2_{jk}\)(slabの分散ね)は\(10^4\)とします。\(\mu_{0j}=\lambda_{0jk}=0, \sigma^2_{\lambda 0jk}=\sigma^2_{\mu 0j} = 10^10\), \(\mathbf{\Phi}, \mathbf{\Psi}\)には無情報事前分布を与えます。識別条件についてひとこと。本来なら\(q^2\)箇所を固定したいけど、このシミュレーションでは主負荷の\(q\)個を固定するだけでよさそうだったので、そうします。
• BSEM-RP。すべての0とxに\(N(0,0.01)\)を与えます。あとはBSEM-SSPと同じ。
• FBS-MI。FBSってのは前向き後向きステップワイズの略。交差負荷がないモデルから出発し、MIが有意な交差負荷をひとつ追加して再推定するのを繰り返し、止まったらすべての交差負荷についてWald検定して有意じゃないやつを一つ削って再推定するのを繰り返す。RのMplusAutomationパッケージで実装しましたー。
• FS。Scheffe型テストで前向き選択の略。[面倒くさいのでメモ省略]
• EFA。主負荷のみの行列をターゲットにして回転する。

　感度分析について… [力尽きたので丸ごとスキップ]

　結果。
　交差負荷のうち0とxを区別できたかどうか、ROC曲線を描いてみると、おおまかにいってBSEM-SSPとFBI-MIの成績がよい。FAとEFAはちょっと劣る。
　パラメータ推定値は…感度分析は…[嗚呼、面倒くさくなってきた… 申し訳ないけどスキップ]

実例
[教員を対象にした横断質問紙で6因子モデルを推定するという話。パス]

考察
　というわけで、CFAにおけるモデル選択を、spike-and-slab事前分布を使った一段階のベイジアン変数選択として解くという方法をご提案しました。ついでに導師MuthenらによるBSEM-RPを拡張したBSEM-RP-VPをご提案しました[ここのくだりはメモをスキップした]。
　BSEM-SSPは全体として、標本サイズと効果量が比較的に大きいときには他のアプローチより成績が良く、そうでないときには倹約的な負荷構造を提供しました。
　重要な交差負荷を検出するという成績でいえばFBS-MIのほうが高いんだけど、false positiveも高い。両方考慮するとBSEM-SSPが優れている。モデル選択の不確実性を定量的に表現できるという利点もある。
　[…FSやEFAについての批評。中略…]
　
　今後の拡張:

• カテゴリカル変数への拡張。
• 独自因子に共分散をいれる。
• spike-and-slabを別の状況で使う。潜在変数間の相関にいれるとか。
• もっと広範なシミュレーション。因子数とか、モデルの誤指定があるときとか。[←これはホントにそうです。正直、Muthen先生たちのBSEMに感じている魅力は、CFAでモデル指定をしくじったときにそれに気づけるだろうという期待感とか、実質的知識を直接表現したCFAモデルが全然推定できない、かといってEFAに知識を入れるのはハードルが高すぎる、でもBSEMは弱めの事前分布を指定すればEFAの代わりに使えるよね… という点にあります。だから、著者らのアプローチでは識別用に\(q^2\)個の負荷を固定しないといけないという点にちょっと引いてしまう。そこの箇所でしくじったら一巻の終わりなんじゃないかな]

　[中略…] 　

　本研究で解決できたわけじゃないけど、BSEM-RPにはいくつか問題点がある。(1)BSEM-RPでは同じデータを使ってBSEMとCFAをやるわけで、そのぶんover-fittingが起きている(SEを過小評価している)はずである。データを分割してCVしたほうがいいのでは。もっとも標本サイズが下がって正確性が下がるけど。(2)一段階目(BSEM)で得たモデル選択の不確実性を二段階目(CFA)で使っていない。
　本研究では、BSEM-RMとBSEM-SSPではちょっと違った解が得られた。どちらを採用すべきなのだろうか？ 思うに、事前分布のちがいはデータに対する信念のちがいなのであり、どちらがよいかは状況による。自分の事前信念にあった事前分布を選びましょう。
　云々。

　。。。面倒くさくて途中から飛ばし読みになっちゃったけど、面白かった。

　ちゃんと読んでないのに申し訳ないんだけど、記録のため、疑問点を３点メモしておく。

　その1, BSEM-SSPでは「どこに非ゼロの負荷があるか」の二値行列\(\mathbf{R}\)の事後分布が得られる。素晴らしい。また\(\mathbf{R}\)の不確実性を考慮した、負荷行列の\(\mathbf{\Lambda}\)の周辺的な事後分布が得られる。これも素晴らしい。でも、分析が探索的なフェイズから確認的なフェイズに移行するにつれ、関心の対象は、(\(\mathbf{P}\)をハイパーパラメータ\(p_{jk}\)の行列として) 負荷行列の周辺的な事後分布 \(P(\mathbf{\Lambda}|\mathbf{Y}, \mathbf{P})\)から 、その分析の結果として得られた最終モデル\(\tilde{\mathbf{R}}\)の下での事後分布\(P(\mathbf{\Lambda}|\mathbf{Y}, \mathbf{P}, \tilde{\mathbf{R}})\)へと移行していくんじゃないかと思う。後者をMCMCサンプルからどうやって推定するのかがよく分からなかった。ほんとに\( P(\lambda_{jk} | \mathbf{Y}, \mathbf{P}, r_{jk} = 1)\)でいいの？
　いや、それは私がきちんと理解できていないだけかもしれないけど、でもやっぱり話の筋道としては、モデルを選んだんだから、選んだ\(\tilde{\mathbf{R}}\)のもとでCFAをやり直すべきなのではないかと思う。そうするとモデル選択の不確実性が表現できないというのはわかりますよ。でも、そもそも分析が探索から検証へと進んでいくプロセスとは、不確実性をモデルの内側から外側(リサーチャーの責任範囲にある実質的推論)へと引き剥がしていくプロセスというか、リサーチャーがより多くの不確実性を実質科学的に抱きしめていくプロセスなんじゃないかと思うわけである。

　その2, BSEM-SSPって、識別のために負荷をいくつか固定しないといけないので、やっぱしモデルを誤指定したときが怖いと思う。その点はやっぱりMuthen的アプローチのほうが有利なんじゃないかしらん

　その3。これは表現するのが難しいんだけど…
　著者のいうように事前分布とは信念の表現なので、spike-and-slabに対応した実質的信念がある場面ではそれが適切だし、そうでない場面では不適切なんだろうと思う。気になるのは、因子分析の文脈でspike-and-slabに対応した実質的信念を持つことがあるのか？という点である。
　回帰の文脈であれば、spike-and-slabで変数選択するというのは自然に感じられる。「この目的変数を説明するために必要な変数はきっと少数しかないはずだ、それがどれかはわからないけど」って思う場面は、確かにありますよね。いっぽう、「いずれの説明変数についても偏回帰係数は0の回りを散らばっているはず」と思うことはそうそうないですよね。だから、ridge回帰のパラメータ推定値を実質的に解釈することは少ないんだと思う。ふつうridge正則化は推定上のテクニカルな事情のために行っているだけだ。
　ではCFAの場合はどうか。spike-and-slabってのは、「因子負荷行列にはきっとたくさんのゼロが埋まっている、それがどこかはわからないけど」という信念に対応する事前分布だと思う。そういう実質的信念を持つ場面が、ちょっと想像できないのである。ある項目が少数の因子にしか負荷を持たないというのは、因子回転の背後にある最適化の原理(「わかりやすさ」の数理的表現)に過ぎず、ユーザが世界について持っている実質的な信念ではないですよね。いっぽうMuthenらの事前分布は、回帰でいえばたしかにridge事前分布だけれども、「因子負荷行列のこことここはゼロだと思うんだけど、でも俺は時々間違えるんだだよね」という実質的信念に対応していると思う。それってすごく自然だよなあ、と感じられる次第である。