カテゴリカル変数について因子分析をしたとき、その因子得点ってどうやって求めるの？ いや、ベイズ推定のときにシミュレーションで求めるのは知っているけど(plausible valuesっていうんでしたっけ)、そうじゃなくて頻度主義的に求めたときには…
　仕事の都合で上記の疑問を抱き、構造方程式モデリングのソフトウェア Mplus の技術文書を読んでみたんだけど、途中でさらに別の文書に遡る羽目になり、ずいぶん時間がかかってしまった。
　以下は Mplus Technical Appendix 11. Estimation of Factor Scoresのメモである。Mplus ver.3の頃の古い文書であり、パラメータのベイズ推定量の話は出てこない。

　おさらいしておくと、Mplusの統計モデルは
$$ \mathbf{y}^*_i = \mathbf{\nu} + \mathbf{\Lambda} \mathbf{\eta}_i + \mathbf{K} \mathbf{x}_i + \mathbf{\varepsilon}_i, \ \ V(\mathbf{\varepsilon}) = \mathbf{\Theta}$$ $$ \mathbf{\eta}_i = \mathbf{\alpha} + \mathbf{B} \mathbf{\eta}_i + \mathbf{\Gamma} \mathbf{x}_i + \mathbf{\zeta}_i, \ \ V(\mathbf{\zeta}) = \mathbf{\Psi}$$ 2本目の構造方程式は次のように書き換えられる。$$ \mathbf{\eta}_i = (\mathbf{I} – \mathbf{B})^{-1} (\mathbf{\alpha} + \mathbf{\Gamma} \mathbf{x}_i + \mathbf{\zeta}_i ) $$ \(\mathbf{\eta}_i\)の期待値と分散は$$ E(\mathbf{\eta}_i | \mathbf{x}_i) = (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\alpha} + (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\Gamma} \mathbf{x}_i $$ $$ V(\mathbf{\eta}_i | \mathbf{x}_i) = (\mathbf{I}-\mathbf{B})^{-1} \mathbf{\Psi} (\mathbf{I}-\mathbf{B})^{\top -1} $$ これを\(\mathbf{\mu}_i, \mathbf{\Sigma}\)と書く。

　\(\mathbf{\eta}_i\)の事後分布は$$ g(\mathbf{\eta}_i | \mathbf{y}_i, \mathbf{x}_i) \propto \phi(\mathbf{\eta}_i | \mathbf{x}_i) f(\mathbf{y}_i | \mathbf{\eta}_i, \mathbf{x}_i)$$ \(\phi(\mathbf{\eta}_i | \mathbf{x}_i)\)は事前分布で、平均\(\mathbf{\mu}_i\), 共分散\(\mathbf{\Sigma}\)のMVN。

　ここまではよろしいですね？
　[あとで気が付いたんだけど、ここでいっている事前分布というのは、モデルのパラメータを推定する前にユーザだかソフトだかが設定した事前分布という意味ではない。モデルのパラメータは推定しました、\(\mathbf{x}_i\)のもとでの\(\mathbf{\eta}_i\)の条件つき期待値と条件つき分散も推定しました、でも\(\mathbf{y}_i\)はまだみてませんし個々の\(\mathbf{\eta}_i\)はまだ推定してません、という意味での事前分布である。\(\mathbf{x}_i\)はわかってるのに\(\mathbf{y}_i\)はわかりませんってところが、いささかシュールというか、混乱を招く点である]
　
　因子得点 \(\mathbf{\eta}_i\)の計算方法は、次の3つのケースに分けて考える必要がある。

ケース1. \(\mathbf{y}\)のすべてが連続変数
　Appendix 2 で説明されているが、このときMplusはモデルを次のように勝手に書き換えるのでした。$$ \mathbf{v}_i = \mathbf{\nu}_v + \mathbf{\Lambda}_v \mathbf{\eta}_{vi} + \mathbf{\varepsilon}_{vi} $$ $$ \mathbf{\eta}_{vi} = \mathbf{\alpha}_v + \mathbf{B}_v \mathbf{\eta}_{vi} + \mathbf{\zeta}_{vi}$$ いまや\(\mathbf{\eta}_v\)のなかには、本来の潜在変数\(\mathbf{\eta}\)だけでなく、一部の\(y\)変数やすべての\(x\)変数と同じ値を持つ謎の変数が含まれてしまっているんだけど、そういうことは気にせず、式の形だけみましょう、形だけ。\(\mathbf{\eta}_v\)の期待値と分散を$$ E(\mathbf{\eta}_v) = (\mathbf{I} – \mathbf{B}_v)^{-1} \mathbf{\alpha}_v = \mathbf{\mu}_v$$ $$ V(\mathbf{\eta}_v) = (\mathbf{I}-\mathbf{B}_v)^{-1} \mathbf{\Psi}_v (\mathbf{I}-\mathbf{B}_v)^{\top -1} = \mathbf{\Sigma}_v$$ と書く。
　このとき、\(\mathbf{\eta}_i\)の事後確率の対数を最大化するのは、通常の回帰法による因子得点推定値と等しい。[\(\mathbf{\eta}_i\)の事後分布が条件つきMVNである以上、そのモードであろうが期待値であろうが、はたまた回帰法による推定値だろうが同じことだ、ということでしょうか]
　すなわち$$ \hat{\mathbf{\eta}}_i = \mathbf{\mu}_i + \mathbf{C} (\mathbf{v}_i -\mathbf{\nu}_v – \mathbf{\Lambda}_v \mathbf{\mu}_v)$$ 因子負荷係数行列は$$ \mathbf{C} = \mathbf{\Sigma}_v \mathbf{\Lambda}^\top_v (\mathbf{\Lambda}_v \mathbf{\Sigma}_v \mathbf{\Lambda}^\top_v + \mathbf{\Theta}_v)^{-1}$$　[さあ困った、さっぱりわからない。落ち着いて考えよう。
　まず第1項の\(\mathbf{\mu}_v\)ってのは\(E(\mathbf{\eta}_v)\)のことですよね。これはわかる。ふつうの因子分析モデルでは因子得点の期待値はゼロだが、Mplusモデルでは必ずしもそうでないので、まず期待値を足している。
　第2項の\(\mathbf{C}\)ってのは因子得点係数行列なんですよね？ ふつうの因子分析でいうところの回帰法ってのは、標準化した観察ベクトルを\(\mathbf{Z}\)、観察変数の相関行列を\(\mathbf{R}\), 因子負荷行列を\(\mathbf{\Lambda}\)、因子間相関行列を\(\mathbf{Q}\)として、因子得点を$$ (\mathbf{Q} \mathbf{\Lambda}^\top \mathbf{R}^{-1}) \mathbf{Z} $$ とする、ってことでしょう？ ってことは、\(\mathbf{C}\)は\(\mathbf{Q} \mathbf{\Lambda}^\top \mathbf{R}^{-1}\)にあたるものでないといけないですよね？ えーと、まず\(\mathbf{\Sigma}_v = V(\mathbf{\eta}_v)\)だから、これは因子間相関行列\(\mathbf{Q}\)。\(\mathbf{\Lambda}^{\top}_v\)は因子負荷の転置行列\(\mathbf{\Lambda}^{\top} \)。残る\(\mathbf{\Lambda}_v \mathbf{\Sigma}_v \mathbf{\Lambda}^\top_v + \mathbf{\Theta}_v\)は… あ、観察変数の相関行列\(\mathbf{R}\)だ。なるほど、これであっているのか…
　ふつうの因子分析モデルなら、因子得点係数行列を(標準化された)観察ベクトルに掛けて、因子得点の推定値を得る。ってことは、\(\mathbf{C}\)の後ろに来るのはきっと観察ベクトル\(\mathbf{v}_i\)だろう、いやまて切片を引いて\(\mathbf{v}_i – \mathbf{\nu}_v\)かな、って思うじゃないですか。ところがここに書かれているのは、\(\mathbf{v}_i – (\mathbf{\nu}_v + \mathbf{\Lambda}_v E(\mathbf{\eta}_v) ) \)なのである。えええ？ 因子得点の期待値に負荷をかけた奴も引くの？ なんで？？ 頭が混乱してきた…]

ケース2. \(\mathbf{y}\)のなかに二値変数ないし順序カテゴリカル変数が含まれている場合
　[この節は逐語訳]
　すくなくともひとつの\(y\)が二値変数ないし順序カテゴリカル変数の場合、次の条件付き独立性が仮定される。$$ f(\mathbf{y}_i | \mathbf{\eta}_i, \mathbf{x}_i) = \prod_{j=1}^p f_j (y_{ij}|\mathbf{\eta}_i, \mathbf{x}_i)$$ その結果\(\mathbf{\Theta}\)は対角行列であると仮定される。アウトカムがカテゴリカルなので、\(\mathbf{\Theta}\)はさらに次のように制約される。$$ \mathbf{\Theta} = \mathbf{\Delta}^{-2} – diag[\mathbf{\Lambda}(\mathbf{I}-\mathbf{B})^{-1} \mathbf{\Psi} (\mathbf{I}-\mathbf{B})^{\top-1} \mathbf{\Lambda}^\top]$$ この制約は、潜在反応変数\(\mathbf{y}^{*}\)がカテゴリカルな\(y\)変数で測定されていることの自然な帰結である。いっぽう、残差共分散がゼロであるという仮定は、因子得点の推定を単純化するための仮定である。
　[はいそこ！くじけない！あきらめない！寝ない！
　右辺第二項は、\(V(\mathbf{\Lambda} \mathbf{\eta}_i | \mathbf{x}_i)\)の対角要素である。つまり、測定方程式の残差分散\(\mathbf{\Theta}\)と、 潜在変数に由来する分散\(V( \mathbf{\Lambda} \mathbf{\eta}_i | \mathbf{x}_i ) \)の和が、\( V(\mathbf{y}^{*}_i | \mathbf{x}_i) \)と一致しなきゃいけないよね、といっているのだ。ではなぜ、\(V(\mathbf{y}^{*}_i | \mathbf{x}_i)\)を\(\mathbf{\Delta}^{-2}\)と書けるのかというと…ええと…それはつまりそういう風にスケーリングしているからなんだろうけど… ええい、よくわからん！]

　カテゴリ\(s = 0,1,2,\ldots, S_j-1\)を持つカテゴリカル変数\(y_j\)について考えよう。閾値を\(\tau_{j,k,0} = -\infty, \tau_{j,k,S_j}= \infty\)とする。\(y_j\)がカテゴリ\(s\)として観察される確率を次のように定義する。$$ f_j(y_{ij}|\mathbf{\eta}_i, \mathbf{x}_i) = \Phi[(\tau_{j,s+1} – \mathbf{\lambda}^\top_j \mathbf{\eta}_i – \mathbf{k}^\top_j \mathbf{x}_i) \theta^{-1/2}_{jj}] – \Phi[(\tau_{j,s} – \mathbf{\lambda}^\top_j \mathbf{\eta}_i – \mathbf{k}^\top_j \mathbf{x}_i) \theta^{-1/2}_{jj}]$$ ここで、\(\mathbf{\lambda}^\top_j\)は\(\mathbf{\Lambda}\)の\(j\)行目、\(\mathbf{k}^\top_j\)は\(\mathbf{K}\)の\(j\)行目、\(\theta_{jj}\)は\(\mathbf{\Theta}\)の\(j\)番目の対角要素である。
　[めんどくさいので先に進むが、たぶん\(y^{*}_j\)が閾値と閾値の間に落ちる確率を求めているのではないかと思う]

　いっぽう連続的な\(y_j\)の場合は$$ f_j(y_{ij}|\mathbf{\eta}_i, \mathbf{x}_i) = \exp(-\frac{1}{2} (y_{ij} – \nu_j – \mathbf{\lambda}^\top_j \mathbf{\eta}_i – \mathbf{k}^\top_j \mathbf{x}_i)^2 \theta^{-1}_{jj})$$である。
　[落ち着け、そんなに難しいことは言っていないはずである。
　\(y_j\)が連続変数の場合、それは条件つき正規分布に従う。その期待値は\(\nu_j + \mathbf{\lambda}^\top_j \mathbf{\eta}_j + \mathbf{k}^\top_j \mathbf{x}_j \)でその分散は\(\theta_{jj}\)である。つまり\(z_{ij} = (y_{ij} – (\nu_j + \mathbf{\lambda}^\top_j \mathbf{\eta}_j + \mathbf{k}^\top_j \mathbf{x}_j)) \theta_{jj}^{-1/2} \)が標準正規分布に従う。その確率密度は$$ \frac{1}{\sqrt{2\pi}} \exp(-\frac{1}{2} z_{ij}^2) = \frac{1}{\sqrt{2\pi}} \exp( -\frac{1}{2} (y_{ij} – \nu_j – \mathbf{\lambda}^\top_j \mathbf{\eta}_j – \mathbf{k}^\top_j \mathbf{x}_j)^2 \theta_{jj}^{-1} ) $$ である。上の\(f_j(y_{ij}|\mathbf{\eta}_i, \mathbf{x}_i)\)はこれに比例している。なるほど？]

　因子得点\(\hat{\mathbf{\eta}}_i \)は事後分布のモードとする。[はいはい、MAP推定量なわけね]
　事後分布のモードは、次の関数\(F\)を\(\mathbf{\eta}_i\)について最小化すると得られる。$$F = 1/2 (\mathbf{\eta}_i – \mathbf{\mu}_i)^\top \mathbf{\Sigma}^{-1}(\mathbf{\eta}_i – \mathbf{\mu}_i) – \sum_j^p \log f_j(y_{ij}|\mathbf{\eta}_i, \mathbf{x}_i)$$ 最小化を行うには反復法が必要である。Mplusは準ニュートン法を使う(\(F\)の1次導関数だけが必要になる)。
　[ううむ、わからん…
　第2項は\(y_{ij}, \mathbf{x}_i\)のもとでの\(\mathbf{\eta}_i\)の対数尤度だと思う。第1項\(1/2 (\mathbf{\eta}_i – \mathbf{\mu}_i)^\top \mathbf{\Sigma}^{-1}(\mathbf{\eta}_i – \mathbf{\mu}_i)\)ってなんなんだろう？ 事前分布 \(\phi(\mathbf{\eta}_i | \mathbf{x})\) の密度の対数に-1を掛けたものだろうか…???]

　[あとで気が付いたのだが、こういうことである。
　最初に書いたように、\(\eta\)の事後分布は事前分布と尤度の積、つまり $$ g(\mathbf{\eta}_i | \mathbf{y}_i, \mathbf{x}_i) \propto \phi(\mathbf{\eta}_i | \mathbf{x}_i) f(\mathbf{y}_i | \mathbf{\eta}_i, \mathbf{x}_i)$$ である。Muthen導師は、この確率密度関数を最大化する\(\eta\)の値を数値的に求めるよ、と仰っているだけなのである。
　事前分布は平均\(\mathbf{\mu}_i\), 分散\(\mathbf{\Sigma}\)の多変量正規分布だから、$$ \phi( \mathbf{\eta}_i | \mathbf{x}_i ) = \frac{1}{ (\sqrt{2 \pi})^m \sqrt{ | \mathbf{\Sigma} | }} \exp \left( – \frac{1}{2} (\mathbf{\eta}_i – \mathbf{\mu}_i)^\top \mathbf{\Sigma}^{-1} (\mathbf{\eta}_i – \mathbf{\mu}_i) \right) $$ 尤度は $$ f(\mathbf{y}_i | \mathbf{\eta}_i, \mathbf{x}_i) = \prod_j^p f(y_{ij} | \mathbf{\eta}_i, \mathbf{x}_i) $$ この2つの積をとって、最初の正規化項を捨てて対数をとると $$ – \frac{1}{2} (\mathbf{\eta}_i – \mathbf{\mu}_i)^\top \mathbf{\Sigma}^{-1} (\mathbf{\eta}_i – \mathbf{\mu}_i) + \sum_j^p \log f(y_{ij} | \mathbf{\eta}_i, \mathbf{x}_i) $$ こいつを最大化するよ、という話である。
　では\(\mathbf{\mu}_i, \mathbf{\Sigma}\)はどこから得るのかというと、$$ \mathbf{\mu}_i = E(\mathbf{\eta}_i | \mathbf{x}_i) = (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\alpha} + (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\Gamma} \mathbf{x}_i $$ $$ \mathbf{\Sigma} = V(\mathbf{\eta}_i | \mathbf{x}_i) = (\mathbf{I}-\mathbf{B})^{-1} \mathbf{\Psi} (\mathbf{I}-\mathbf{B})^{\top -1} $$ に\(\hat{\mathbf{\alpha}}, \hat{\mathbf{B}}, \hat{\mathbf{\Gamma}}, \hat{\mathbf{\Psi}}\)を放り込んで得るのだと思う]

ケース3. 混合モデルの場合
　混合モデルの場合は、$$E(\mathbf{\eta}_i | \mathbf{x}_i, \mathbf{y}_i) = E_{c|x,y}[E(\mathbf{\eta}_i|\mathbf{x}_i,\mathbf{y}_i,\mathbf{c}_i)]$$を求める。\(E(\mathbf{\eta}_i|\mathbf{x}_i,\mathbf{y}_i,\mathbf{c}_i)\)はクラスごとのパラメータ値を使って、ケース2と同じやり方で求める。\(E_{c|x,y}\)は\(\mathbf{x}, \mathbf{y}\)のもとでの\(\mathbf{c}\)の条件つき分布を通じた期待値という意味で、結局これは、クラスを事後確率を使って合計することに帰着する。

——
　このたび仕事のなかで、「Mplusでカテゴリカル因子分析をしたとき、因子得点係数行列をアウトプットできないのはなぜ？」と疑問に思い、このAppendixを読んでみた次第である。
　なるほど。連続変数だけの因子分析なら伝統的な回帰法で係数行列を求めます。カテゴリカル変数がひとつでも入っていたら、因子得点の事後分布のモードを数値的に求めます。ってことですね。

　疑問その1：カテゴリカル変数入りの因子分析でも、因子共分散行列\( \mathbf{\Sigma}\), 負荷行列\(\mathbf{\Lambda}\)、独自分散の対角行列\(\mathbf{\Theta}\)は出力されるじゃないですか。そこから因子得点係数行列 $$ \mathbf{C} = \mathbf{\Sigma} \mathbf{\Lambda}^\top (\mathbf{\Lambda} \mathbf{\Sigma} \mathbf{\Lambda}^\top + \mathbf{\Theta})^{-1}$$ を求めてはいかんの？
　あっ、そうか、それを求めるのはおまえの勝手だが、それを掛ける相手としては\(\mathbf{y}_i\)じゃなくて\( \hat{\mathbf{y}}^{*}_i \)が必要なのに、いったい何に使うつもりだい？ってことかなあ…

　疑問その2: カテゴリカル変数入りの因子分析では、因子得点を簡単に推定できない、ってことはわかったけど、近似でもいいから簡単に推定したいじゃないっすか。なにかよい方法はないもんだろうか。
　あっ、そうだ、これってIRTでも同じことだ。項目の特性関数が推定できていても、項目の正誤からその人の能力を簡単に推定することはできないけど、でも実務的には、特性関数に基づき項目になんらか配点を与え、合計得点で能力を近似したいじゃないですか。そういうときの配点のしかた、どっかで読んだことがあるような… 勤務先の本棚をみればわかるんだけど… うぐぐぐ…

　疑問その3: 要するに、Mplusが出力する因子得点は、ケース1(すべての従属変数が連続)では回帰法による推定値(いいかえれば事後分布の期待値の推定値)となり、ケース2(従属変数の中に二値変数か順序変数が含まれている)では事後分布のモードを反復推定する、ということだと思う。
　ところで、このAppendixよりおそらく新しいであろう別の文書によれば、Mplusが出力する因子得点は、パラメータをML推定したときには事後分布の期待値つまりEAP推定量、パラメータをWLS推定したときには事後分布のモードつまりMAP推定量になるのだそうだ。Mplus 8では、TYPE=GENERALのときのデフォルトのパラメータ推定量は、すべての\(y\)が連続変数の時はML, 二値変数か順序変数が含まれていたらWLSMVだから、デフォルトでは辻褄があっている。
　では、ケース1でパラメータ推定量を明示的にWLSと指定した場合はどうなるんだろうか。この文書に書かれているように、依然として回帰法なのか？ それとも別の文書のほうが正しくて、ケース1であっても因子得点はMAP推定量となるのか？ あ、でもどのみち因子得点の事後分布は条件付きMVNだから、回帰法(EAP推定量)でもMAP推定量でも同じこと、なのかなあ…
　ケース2で明示的にMLと指定した場合、因子得点の推定量はどうなるんだろう？ 依然としてMAP推定量の反復計算なのだろうか、それとも別の文書のほうが正しくて、EAP推定量になるのだろうか。仮に後者だとすると、\(\mathbf{y}^{*}_i\)がわからない以上回帰法というわけにいかないから、ここに書かれてない別の最適化関数が必要になるのではなかろうか。

2021/01/01追記: このメモをブログにアップして、一晩寝て、先ほど(大晦日の晩)に風呂に浸かっていたら、ケース2の最適化関数の由来についてハタと気がついた。余計なことを考えすぎだ。もっと素直に考えるべし。というわけで、いくつかメモを追記した。