構造方程式モデリングのソフトウェア Mplus を使っていて、よくわからないことが起きて解説を辿っていくと、最終的にすごく基礎的な説明までさかのぼることになり、難しすぎて断念することがある。今回も仕事の都合で、ちょっとそういうことが起きかけて、大変イライラしたので仕事を放りだし、初期の技術文書をメモを取りながら読んだ。
　以下は Mplus Technical Appendix 2. The general modeling framework からのメモである。Mplus ver.3の頃の古い文書なので(Mplusの最新版は8.5)、当然ながらベイジアンな話は出てこない。

変数
　Mplusでは変数を従属変数\(y\)と独立変数\(x\)に分ける。[\(x\)は完全に外生的な変数だけ、つまり、そこからパスは伸びるが刺さらないような変数だけ。あとはみな\(y\)である]
　\(y\)は連続変数かも知れないし、二値変数もしれないし、カテゴリ数\(C\)の順序変数かも知れない。個々の\(y\)の裏には、それに対応する連続潜在変数\(y^{*}\)があると考える。

• \(y\)が連続変数なら\(y = y^{*}\)である。
• \(y\)が二値変数の場合、\(y^{*} > \tau\) のとき \(y=1\)、そうでないとき0とする。
• \(y\)が順序変数の場合、カテゴリ\(0, 1, \ldots, C-1\)の閾値を\(\tau_0, \tau_1, \ldots, \tau_C\)として、\(\tau_0 = -\infty, \tau_C = \infty\)とし、\(\tau_c < y^{*} \leq \tau_{c+1}\)のとき\(y = c\)とする。

　以下、\(i\)を個体のインデクス, \(j\)を変数のインデクスとする。

統計モデル
　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}$$　一本目が測定方程式。左辺の\(\mathbf{y}^*_i\)は\(\mathbf{y}_i\)の裏にある潜在変数ベクトル。\(\mathbf{\eta}_i\)は\(m\)次元の潜在変数ベクトル、\(\mathbf{x}_i\)は\(q\)次元の独立変数ベクトル。測定誤差ベクトル\(\mathbf{\varepsilon}_i\)は他の変数とは独立。[一般化線型モデルであろうがなんだろうが測定誤差項はつくわけだ]
　二本目が構造方程式。\(\mathbf{B}\)の対角要素は0。\(\mathbf{I} – \mathbf{B}\)が正則だと仮定する[下でその理由がわかる]。

　構造方程式では両辺に\(\mathbf{\eta}\)が入っているので、左辺に移そう。$$ (\mathbf{I} – \mathbf{B}) \mathbf{\eta}_i = \mathbf{\alpha} + \mathbf{\Gamma} \mathbf{x}_i + \mathbf{\zeta}_i$$ \(\mathbf{I} – \mathbf{B}\)が正則だから $$ \mathbf{\eta}_i = (\mathbf{I} – \mathbf{B})^{-1} (\mathbf{\alpha} + \mathbf{\Gamma} \mathbf{x}_i + \mathbf{\zeta}_i ) $$ である。

　Mplusの枠組みでは、\(\mathbf{x}\)のもとで\(\mathbf{y}^{*}\)が条件つき正規性を持つ(周辺分布が正規なわけではない)。期待値は [測定方程式の期待値に\(\mathbf{\eta}_i\)の期待値を代入して] $$ E(\mathbf{y}^{*}_i | \mathbf{x}_i ) = \mathbf{\nu} + \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\alpha} + \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\Gamma} \mathbf{x}_i + \mathbf{K} \mathbf{x}_i$$ 共分散は [一般に\(V(\mathbf{Ax}) = \mathbf{A} V(\mathbf{x}) \mathbf{A}^\top\)だから、\(\mathbf{A} = (\mathbf{I} – \mathbf{B})^{-1}\)と略記して$$ V(\mathbf{\eta}_i|\mathbf{x}_i) = \mathbf{A} V(\mathbf{\eta}) \mathbf{A}^\top $$ $$ V(\mathbf{\Lambda} \mathbf{\eta}_i|\mathbf{x}_i) = \mathbf{\Lambda} \mathbf{A} V(\mathbf{\eta}) \mathbf{A}^\top \mathbf{\Lambda}^\top $$ $$ V(\mathbf{\Lambda} \mathbf{\eta}_i + \mathbf{\varepsilon}_i|\mathbf{x}_i) = \mathbf{\Lambda} \mathbf{A} V(\mathbf{\eta}) \mathbf{A}^\top \mathbf{\Lambda}^\top + V(\mathbf{\varepsilon})$$ というわけで] $$ V(\mathbf{y}^{*}_i | \mathbf{x}_i) = \mathbf{\Lambda}(\mathbf{I}-\mathbf{B})^{-1} \mathbf{\Psi} (\mathbf{I} – \mathbf{B})^{\top-1}\mathbf{\Lambda}^{\top} + \mathbf{\Theta}$$

スケーリング
　Mplusでは、潜在反応変数\(\mathbf{y}^{*}_i\)をスケーリングする対角行列\(\mathbf{\Delta}\)というのを考える。\( \mathbf{y}^{*}_{si} = \Delta \mathbf{y}^{*}_{i} \)と書く。

• \(y\)のすべてがカテゴリカルだったら、\(\mathbf{\Delta}\)の対角には、\(\mathbf{x}\)のもとでの\(\mathbf{y}^{*}\)の条件付きSDの逆数をいれる。式で書くと $$ diag[\Delta] = diag[\mathbf{\Lambda}(\mathbf{I}-\mathbf{B})^{-1} \mathbf{\Psi} (\mathbf{I} – \mathbf{B})^{\top-1}\mathbf{\Lambda}^{\top} + \mathbf{\Theta}]^{-1/2}$$ということになる。\( \mathbf{y}^{*}_{si} \)の分散は1になる。\(\mathbf{\Theta}\)は識別できなくなるので、\(\mathbf{\Delta} = \mathbf{I}\)だということにして、残りの部分が\(\mathbf{\Theta}\)だということにする。つまり $$ diag[\mathbf{\Theta}] = \mathbf{I} – diag[\mathbf{\Lambda}(\mathbf{I}-\mathbf{B})^{-1} \mathbf{\Psi} (\mathbf{I} – \mathbf{B})^{\top-1}\mathbf{\Lambda}^{\top}]$$　同じ\(y\)変数を時間を通じて比べたり、群を通じて比べたりするときには、\(\mathbf{\Delta}\)の要素は1時点目ないし1群目で1に標準化しておき、残りの時点ないし群では推定する。このように\(\mathbf{\Delta}\)を\(\mathbf{I}\)に固定しない場合は、$$ diag[\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{\Delta}^{-2}\)の対角要素は\(V(\mathbf{y}^{*}|\mathbf{x})\)の分散である。
• \(y\)の一部がカテゴリカルだったら、\(\mathbf{\Delta}\)の対角要素は、カテゴリカルな\(y\)については上記の通りとし、連続的な\(y\)については1とする。
• \(y\)のすべてが連続的だったら、\(\mathbf{\Delta}\)の対角要素に\(y\)のSDをいれれば、\(\mathbf{y}^{*}_i\)は分散1に標準化される。

　[ここでのスケーリングの説明、ものすごくわかりにくいと思う… まあとにかく、\(\mathbf{y}^{*}_i\)を\(\Delta\)でスケーリングしたりしなかったりするわけだ]

プロビット傾きとプロビット残差相関
　測定方程式に\(\mathbf{\eta}_i\)を代入して $$ \mathbf{y}^{*}_i = \mathbf{\nu} + \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\alpha} + \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\Gamma} \mathbf{x}_i + \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\eta}_i + \mathbf{K} \mathbf{x}_i + \mathbf{\varepsilon}_i$$ これを $$ \mathbf{y}^{*}_i = \mathbf{\pi}_0 + \mathbf{\Pi} \mathbf{x}_i + \mathbf{\delta}_i$$と書くことにする。
　\(\mathbf{\nu}\)は閾値パラメータと切り離して識別できないから、通常は1に固定する。\(\mathbf{\delta}_i\)の分散を$$ V(\mathbf{\delta}_i) = \mathbf{\Lambda}(\mathbf{I}-\mathbf{B})^{-1} \mathbf{\Psi} (\mathbf{I} – \mathbf{B})^{\top-1}\mathbf{\Lambda}^{\top} + \mathbf{\Theta} = \mathbf{\Omega}$$ と書く。
　\(\mathbf{y}_i^{*}\)をスケーリングしたあとの期待値と分散を $$ E(\mathbf{y}^{*}_{s} | \mathbf{x}) = \mathbf{\Delta} [\mathbf{\pi}_0 + \mathbf{\Pi} \mathbf{x}] = \mathbf{\mu}^{*}(\mathbf{x}) $$ $$ V(\mathbf{y}^{*}_s | \mathbf{x}) = \mathbf{\Delta} \mathbf{\Omega} \mathbf{\Delta} = \mathbf{\Sigma}^{*}$$と書く。\(\mathbf{\Sigma}^{*}\)の対角は1である。
　\(\mathbf{\Delta} \mathbf{\Pi}\)のことをプロビット傾きという。\(\mathbf{\Sigma}^{*}\)の非対角要素(つまり相関)のことをプロビット残差相関という。

\(y\)の条件付き確率
　\(\mathbf{y}^{*}\)の条件付き正規性の仮定があるから、\(y\)の条件付き確率は下式で表現できる。閾値\(\tau_j\)に\(\mathbf{\Delta}\)の\(j\)番目の対角要素を掛けたのを\(\tau^{*}_j\)として、$$ P(y_j = 1 | \mathbf{x}) = \int^{\infty}_{\tau^{*}_j-\mu^{*}_j(\mathbf{x})} \phi_1(y^{*}_j | \mathbf{x}) d y^{*}_j$$ \(\phi_1\)は単変量の標準正規密度関数である。
　2変量の同時確率はどうなるかというと $$ P(y_j = 1, y_k = 1| \mathbf{x}) = \int^{\infty}_{\tau^{*}_j-\mu^{*}_j(\mathbf{x})} \int^{\infty}_{\tau^{*}_k-\mu^{*}_k(\mathbf{x})} \phi_2(y^{*}_j,y^{*}_k | \mathbf{x}) d y^{*}_k d y^{*}_j $$ってことになりますね。\(\phi_2\)は2変量の標準正規密度関数である。

推定するパラメータ
　推定するパラメータをまとめておく。

• \(\mathbf{\tau}\): モデルに出てくる閾値のベクトル。二値変数なら閾値を1個、\(C\)カテゴリのカテゴリカル変数なら閾値を\(C-1\)個持つ。
• \(\mathbf{\nu}\): 切片ベクトル。長さは変数の数\(p\)。
• \(\mathbf{\Lambda}\): 測定方程式における傾き(負荷)の行列。\(p \times m\)。
• \(\mathbf{\Theta}\): 測定方程式の残差の共分散行列。\(p \times p\)。
• \(\mathbf{\alpha}\): 潜在変数の平均・切片のベクトル。長さ\(m\)。
• \(\mathbf{B}\): 潜在変数の間の回帰の傾きの行列。対角はゼロ。\(m \times m\)。
• \(\mathbf{\Gamma}\): 潜在変数の\(x\)変数に対する回帰係数の行列。\(m \times q\)。
• \(\mathbf{\Psi}\): 潜在変数の残差の共分散行列。\(m \times m\)。
• \(\mathbf{\Delta}\): スケーリング・ファクター。対角行列、\(p \times p\)。

　[もう一度書いておくけど、いまやモデルは測定方程式も構造方程式も一緒くたにして以下のように再定式化されている。
$$ \mathbf{y}^{*}_i = \mathbf{\pi}_0 + \mathbf{\Pi} \mathbf{x}_i + \mathbf{\delta}_i$$ ここで

• \( \mathbf{\pi}_0 \equiv \mathbf{\nu} + \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\alpha}\)
• \(\mathbf{\Pi} \equiv \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\Gamma} + \mathbf{K}\)
• \(\mathbf{\delta}_i \equiv \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\eta}_i + \mathbf{\varepsilon}_i\)

この書き方に従って、\(\mathbf{y}^{*}_i\)の期待値と分散を求めると$$E(\mathbf{y}^{*}_i | \mathbf{x}_i) = \mathbf{\pi}_0 + \mathbf{\Pi} \mathbf{x}_i $$ $$ V(\mathbf{y}^{*}_i | \mathbf{x}_i) = \mathbf{\Lambda}(\mathbf{I}-\mathbf{B})^{-1} \mathbf{\Psi} (\mathbf{I} – \mathbf{B})^{\top-1}\mathbf{\Lambda}^{\top} + \mathbf{\Theta} = \mathbf{\Omega}$$ である。スケーリングした後は $$ E(\mathbf{y}^{*}_{si} | \mathbf{x}_i) = \mathbf{\Delta} [\mathbf{\pi}_0 + \mathbf{\Pi} \mathbf{x}_i] = \mathbf{\mu}^{*}(\mathbf{x}_i) $$ $$ V(\mathbf{y}^{*}_{si} | \mathbf{x}_i) = \mathbf{\Delta} \mathbf{\Omega} \mathbf{\Delta} = \mathbf{\Sigma}^{*}$$ である。]

実装(1) 二値ないしカテゴリカルな\(y\)があるとき
　いよいよ計算上の実装について。まず、\(y\)には二値ないしカテゴリカル変数が含まれているものとする。
　計算上は、モデルを次の3つのパートに分ける。

　パート1、平均(閾値・切片)構造。$$\mathbf{\Delta}^{*} [ \mathbf{K}_\tau \mathbf{\tau} – \mathbf{K}_\nu [\mathbf{\nu} + \mathbf{\Lambda}(\mathbf{I}-\mathbf{B})^{-1} \mathbf{\alpha}]]$$ 　\(\mathbf{K}_\tau, \mathbf{K}_\nu\)は選択行列。
　[式をよく見ると、これは \( \mathbf{\Delta}^{*} [ \mathbf{K}_\tau \mathbf{\tau} – \mathbf{K}_\nu \mathbf{\pi}_0] \)である。
　想像するにこういうことだと思う。たとえば\(y\)変数が3つあり、\(y_1\)が連続変数, \(y_2\)が二値変数でその閾値が\(\tau_2\), \(y_3\)が3カテゴリのカテゴリカル変数でその閾値が\(\tau_{31}, \tau_{32}\) だとしよう。\( \mathbf{\tau} = (\tau_2, \tau_{31}, \tau_{32})^{\top} \)である。いっぽう\(\mathbf{\pi}_0 \equiv \mathbf{\nu} + \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\Gamma}\)は長さ3のベクトルで、これを\(\mathbf{\pi}_0 = (\pi_{0,1}, \pi_{0,2}, \pi_{0,3})^{\top} \)としよう。
　\(\mathbf{K}_\tau\)は、かっこ内の第1項が\(\mathbf{K}_\tau \mathbf{\tau} = (0, \tau_2, \tau_{31}, \tau_{32})^\top\)となるような\(4 \times 3\) の0-1行列。
　\(\mathbf{K}_\nu\)は、かっこ内の第2項が\(\mathbf{K}_\nu \mathbf{\pi}_0 = (\pi_{0,1}, \pi_{0,2}, \pi_{0,3}, \pi_{0,3})^{\top} \)となるような\(4 \times 3\) の0-1行列。
　\(\mathbf{\Delta}^{*}\)は、対角に\(\mathbf{\Delta}\)の対角要素の1個目, 2個目, 3個目, 3個目を持つ対角行列。
　たとえば\(y\)がすべて二値変数だとしたら、これは\(\mathbf{\Delta} [\mathbf{\tau} – \mathbf{\pi}_0]\)、つまり「\(\mathbf{y}\)の閾値から\(\mathbf{y}^{*}\)の切片を引いてスケーリングしたやつ」である… ってことで、あってますかね？]

　パート2、傾き構造。
　行列\(\mathbf{A}\)の要素(対称行列だったら下三角の要素)を行順にとってきてベクトルにするオペレータを\(vec[\mathbf{A}]\)として、$$ vec[\mathbf{\Delta} \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\Gamma}] $$　[要するにこれは、\(\mathbf{\Pi} \equiv \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\Gamma} + \mathbf{K}\) の第1項をスケーリングしたものですね。第2項\(\mathbf{K}\)は含んでいない。原文には「\(\mathbf{K}\)は連続的\(y\)変数に関して後述するかたちで扱われる」と書いてあるのだが、どうもよくわからない… まあ先にすすみましょう]

　パート3、共分散構造。$$ \mathbf{K}_\sigma vec[\mathbf{\Delta} [ \mathbf{\Lambda} (\mathbf{I} – \mathbf{B})^{-1} \mathbf{\Psi} (\mathbf{I} – \mathbf{B})^{\top -1} \mathbf{\Lambda}^\top + \mathbf{\Theta}] \mathbf{\Delta}] $$ ここで\(\mathbf{K}_\sigma\)は、それが事前に掛けている(premultiplies)ベクトルから要素を選び、もし対応する\(y\)変数がカテゴリカルだったら対角要素を含まないようにしている。[なにをいってんだかさっぱりわからん… vecオペレータの内側は、\(\mathbf{y}_i^{*}\)をスケーリングしたあとの分散\(V(\mathbf{y}^{*}_{si} | \mathbf{x}_i) = \mathbf{\Delta} \mathbf{\Omega} \mathbf{\Delta} \)ですよね。上の例でいえば\(3 \times 3\)の対称行列である。よくわからんが、これをなんだかイイ感じに\(4 \times 4\)の行列に拡張しているのだろうか]

実装(2) すべての\(y\)が連続のとき
　すべての\(y\)が連続だと、話が変わってくる。
　まず、パート1から閾値\(\mathbf{\tau}\)がなくなる。[どういうこと？ 単に\(\mathbf{\Delta} \mathbf{\pi}_0\)になる、っていうことであってるの？ もっと親切に書いてよ…]
　パート2はまるごとなくなる。プロビット傾きはもう考えなくていいから。
　パート3は共分散構造ないし相関構造を捉える。計算にあたっては、\(\mathbf{K}\)や\(\mathbf{\Gamma}\) [\(\mathbf{x}\)の測定方程式での係数と構造方程式での係数]を求めるのではなく、次の測定方程式・構造方程式に基づいた、より少数のパラメータの配列を用いる。

　まずは測定方程式から。$$ \mathbf{v}_i = \mathbf{\nu}_v + \mathbf{\Lambda}_v \mathbf{\eta}_{vi} + \mathbf{\varepsilon}_{vi} $$ ここで、\(\mathbf{v}_i = (\mathbf{y}^\top_i, \mathbf{x}^\top_i )^\top\)である。[つまり、あらゆる観察変数を縦に積んだベクトルだ。発想をがらっと切り替えないといけない。たとえば、\(y_1\)は\(x\)に回帰し、\(y_2\)は\(y_1\)と\(x\)に回帰する、というパスモデルだったら、左辺は3変数を縦に積んだベクトル\(\mathbf{v}_i = (y_{2i}, y_{1i}, x)^\top \)である]
　第1項は $$ \mathbf{\nu}_v = \left( \begin{array}{c} \mathbf{\nu}_{nd} \\ \mathbf{0} \\ \mathbf{0} \end{array} \right)$$ である。[こいつは縦ベクトル。\(\mathbf{\nu}_{nd}\)はここが初出だけど、おそらくオリジナルの\(\mathbf{\nu}\)そのものではない。\(y\)変数を、完全に従属なタイプの\(y\)変数と、他の\(y\)変数の説明変数になるような\(y\)変数に分け、前者の切片だけを表しているのではないかと思う。後者の切片は構造方程式のほうにいれるのだ]
　第2項。\(\mathbf{\eta}_{vi} = (\mathbf{\eta}^\top_i, \mathbf{\eta}^\top_{yi}, \mathbf{\eta}^\top_{xi})^\top\)である。\(\mathbf{\eta}^\top_{y}\)は、\(y\)変数のうち\(y\)が指標でない他の変数に回帰している変数に対応する新しい潜在変数のベクトルであり、\(\mathbf{\eta}^\top_{x}\)はそれぞれの\(x\)変数に対応する潜在変数ベクトルである。係数は$$ \mathbf{\Lambda}_v = \left( \begin{array}{ccc} \mathbf{\Lambda}_{nd} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{I} \end{array} \right)$$ [つまり、本来の潜在変数\(\mathbf{\eta}_i\)の下に、他の\(y\)変数の説明変数になるタイプの\(y\)変数そのものを表す潜在変数と、\(x\)そのものを表す潜在変数を積んでいる]
　第3項の共分散は$$ \mathbf{\Theta}_v = \left( \begin{array}{ccc} \mathbf{\Theta}_{nd} & & symm. \\ \mathbf{0} & \mathbf{0} & \\ \mathbf{0} & \mathbf{0} & \mathbf{0} \end{array} \right)$$

　構造方程式は$$ \mathbf{\eta}_{vi} = \mathbf{\alpha}_v + \mathbf{B}_v \mathbf{\eta}_{vi} + \mathbf{\zeta}_{vi}$$ 第1項は $$ \mathbf{\alpha}_v = \left( \begin{array}{c} \mathbf{\alpha} \\ \mathbf{\nu}_d \\ \mathbf{\mu}_x \end{array} \right)$$ \(\mathbf{\mu}_x\)は\(x\)変数の平均ベクトル。推定の際には標本平均に固定する。[ああなるほど。本来の\(\mathbf{\nu}\)を、完全な従属変数である\(y\)のパート\(\mathbf{\nu}_{nd}\)と、他の\(y\)の説明変数になるタイプの\(y\)のパート\(\mathbf{\nu}_{d}\)にわけて、後者を状態方程式、じゃなかった、構造方程式のほうにいれたわけね。要はあれだ、状態空間モデルみたいなものだと思えばいいのか]
　第2項の係数は$$ \mathbf{B}_v = \left( \begin{array}{ccc} \mathbf{B} & \mathbf{B}_{\eta d} & \mathbf{\Gamma} \\ \mathbf{\Lambda}_d & \mathbf{B}_d & \mathbf{K}_d \\ \mathbf{0} & \mathbf{0} & \mathbf{0} \end{array}\right) $$ [1行目から見ていこう。1列目は本来の\(\mathbf{\eta}\)の遷移行列だからオリジナルの\(\mathbf{B}\)。2列目は、「他の\(y\)の説明変数になる\(y\)」から本来の\(\mathbf{\eta}\)へのパス係数で… 頭が混乱してきちゃうけど、これは元のモデルではどこに書かれているんだろう? やっぱり\(\mathbf{B}\)かな？ 3列目は\(x\)から本来の\(\mathbf{\eta}\)へのパス係数、つまりオリジナルの\(\mathbf{\Gamma}\)。
　2行目。1列目は本来の\(\mathbf{\eta}\)から「他の\(y\)の説明変数になる\(y\)」へのパス係数、つまりオリジナルの\(\mathbf{\Lambda}\)の一部。2列目は「他の\(y\)の説明変数になる\(y\)」同士の遷移行列だから… よくわかんなくなってきたけどやっぱり\(\mathbf{B}\)の一部なんだろうな。3列目は\(x\)から「他の\(y\)の説明変数になる\(y\)」へのパス係数だから、当然\(\mathbf{K}\)の一部である。
　3行目はすべてゼロになる。\(x\)にパスは刺さんないから]
　第3項の共分散は$$ \mathbf{\Psi}_v = \left( \begin{array}{ccc} \mathbf{\Psi} & & symm. \\ \mathbf{\Psi}_{d\eta} & \mathbf{\Theta}_d & \\ \mathbf{\Psi}_{x\eta} & \mathbf{\Psi}_{xd} & \mathbf{\Sigma}_{xx} \end{array} \right)$$ \(\mathbf{\Sigma}_{xx}\)は\(x\)変数の共分散行列で、推定の際には標本共分散に固定する。
　[よくわからん。こうして再定式化したことで、どういういいことがあるんだろうか…]

　… いやあ、むずかしかった…
　Muthen導師の説明は、いつもはすごくわかりやすいのに、ここでの説明は大変にわかりにくかった。2004年の文章だから、導師もまだ修業が足りなかったんですかね。(こういうとき、自分の能力は棚の上にあげないといけない。高ーい高ーい棚の上に)

　結局、最後の再定式化のなかで、二値・カテゴリカルな\(y\)があるとき\(\mathbf{K}\)がどこに行っちゃうのか、よくわからなかった。急いでどうしても知りたいわけじゃないんだけど、自分の無能さを突き付けられたという感じだ。
　もうひとつ疑問に思ったのは、ほんとにこれ、読まないといけないものだったのか、ということ。別の解説で引用されていたので仕方なく目を通したんだけど、内部でどう実装しているかなんて、私のレベルのユーザにはどうでもよいのではなかろうか。$$ \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}$$ で\(\mathbf{y}^{*}_i\)は\(\mathbf{x}_i\)のもとでMVN、というあたりまで知っていればよいのでは、と…
　なんだかがっくりだ。まあいいけどさ。俺には風呂と布団がある。