Asparouhov, T., Muthen, B. (2021) Residual Structural Equation Models. Mplus.
泣く子も黙る構造方程式モデリング用ソフトウェアMplusは、バージョンアップのたびになんらかの先進的すぎる謎機能を搭載してくることで有名である(私のなかで)。
今月リリースされたVersion 8.7では、えーと、従来は残差動的構造方程式モデルのベイズ推定のみが可能であったラグ変数残差間回帰が単一レベルモデルの最尤推定・重みつき最小二乗推定・ベイズ推定へと拡張され、パネルデータのランダム切片クロスラグモデルならびにランダム切片自己回帰移動平均モデルの推定が可能となった、のだそうだ。
はあ、そうですか、と虚ろな目でディスプレイに相槌を打つ私である。わからない。もうなにもわからない。
いちいち腑抜けているのもしゃくに触るので、リリースと同時に発表されたPDFファイルを読んでみた。Technical Paperという位置づけだが、いずれ論文化されるのだろう。
背景としては… Mplusの中の人であるMuthen導師とその仲間たちは、Version 8の頃から動的SEM(DSEM), 8.1の頃から残差動的SEM(RDSEM)という枠組みを提唱している。今回提案されるRSEMはその拡張版、ということらしい。
1. イントロダクション
[略]
2. 一般的なRSEMモデル
基本的なSEMモデルはこう書ける。\(Y\)を連続観察従属変数のベクトル、\(\eta\)を連続潜在変数のベクトル、\(X\)を共変量ベクトルとして$$ Y = \nu + \Lambda \eta + K X + \varepsilon $$ $$ \eta = \alpha + B \eta + \Gamma X + \xi$$ さて、RSEM(残差SEM)モデルとはこういうモデルである。残差変数\(\varepsilon, \xi\)を縦に積んで\(R\)としよう。$$ R = B_r R + \zeta, \ \ \zeta \sim N(0, \Psi)$$ [ここまで読んで、ああああなるほどおおお、と小声で叫んだ次第である。残差に構造を持たせる気なのか！ これ、クロス・セクショナル・データのSEMではなく、たとえばパネルデータの分析で、\(Y\)がある個体の複数時点の測定値からなるベクトルだと思うとわかりやすい。\(B_r\)の非対角要素をうまいこと使えばさあ！自己相関を持つ誤差項なんかをモデル化できるって寸法よ！そうですよね導師！]
カテゴリカル観察変数の場合は、単に観察変数を、その背後に仮定される連続変数に置き換えればよい。この場合は閾値パラメータもモデルに入る。
RSEMモデルとは、\(R\)の共分散行列を$$ ((I-B_r)^{-1})^\top \Psi (I-B_r)^{-1}$$ とした標準的SEMモデルだと考えてもよい。
ついでに拡張RSEMモデルというのも定義しておこう。$$ Y = \nu + \Lambda \eta + K X + B_{1r} R + \epsilon$$ $$ \eta = \alpha + B \eta + \Gamma X + B_{2r} R + \zeta$$ $$ R = B_r R + \zeta, \ \ \zeta \sim N(0, \Psi)$$ Mplusではこのモデルも推定できる(ML, WLSMV推定が使えるがBayes推定はできない)。
RSEMに対応するため、Mplus言語も拡張したぞ。Y の残差変数を Y^ と書く。残差変数同士を ON で回帰したり WITH で相関を持たせたりできる。
推定にあたっては、RSEMモデルは標準的SEMモデルに変換される。残差変数が潜在変数(負荷は1)と新しい残差変数にかわり、後者の分散が0に固定されるわけだ。なので、ML推定とWLSMV推定についていえば、推定プロセスは変わらない。
ML推定は連続従属変数の時しか使えない。なぜなら、Mplusではカテゴリカル変数の残差が直接に相関を持てるのはWLMSVかBayesだけだからだ。
カテゴリカル従属変数の時はthetaパラメータ化を使うことになる。新しい残差変数の残差分散を0に固定したいから。[カテゴリカル従属変数の場合、その後ろにある潜在変数の分散をパラメータとみることもできるし、残差分散をパラメータとみることもできるわけで、どっちにするかというややこしい話が生じるのである。ぶひー]
Bayes推定の場合は俄然話がややこしくなる。モデルを標準的SEMモデルに変換したとして、新しい残差変数の分散をどうやって0に固定するかという問題が生じるからである。
そこで新しいやりかたを考えた。詳しくは2010年のTech. Appendixをみてほしいんだけど…
標準的なSEMモデルであれば、すべての構造パラメータと切片は、他のすべてのモデルパラメータと潜在変数で条件づけた正規分布を使ってMCMC更新される。しかしRSEMの場合は次の手順をとる。
- \(Y, \eta\)の式を、 誤差共分散行列が\( ((I-B_r)^{-1})^\top \Psi (I-B_r)^{-1} \)である標準的SEMモデルとみなして、パラメータをMCMC更新する。
- \(R\)を求める。
- \(R\)の式を推定する。
[恥ずかしながら、これをMCMCのiterationごとにやるという話なのか、それとも順に片付けていくという話なのかわからなかった。後者かなあ…]
3.1 自己回帰誤差構造を持つ成長モデル
個体を\(i\), 時点を\(t\)とする。線形成長モデルは$$ Y_{it} = I_i + S_i t + \epsilon_{it} $$ ランダム切片\(I_i, S_i\)は正規分布するとします。
ここに残差の自己回帰構造をいれるやりかたを考えてみよう。次の3つがある。
- \(\epsilon_{it} \sim N(0, \theta)\)なんだけど、\(Cor(\epsilon_{i,t_1}, \epsilon_{i, t_2}) = \rho ^{|t_1 – t_2|}\)とするやり方。\(\theta, \rho\)はモデルパラメータである。Mplusでは
MODEL CONSTRAINTで書ける。[書けるって、導師… ローデータの行を個体にしておき、すべての\(t\)のペアについてPWITHを書けってことですよね…それ結構面倒くさいですよね]
このやりかたは変数がひとつでせいぜい50時点くらいまでの場合に限られるし、残差分散や自己回帰パラメータを時変させられない。 - 以前提案した、DSEM(動的SEM), RDSEM(残差動的SEM)というフレームワーク。$$ \epsilon_{it} = \rho_i \epsilon_{i, t-1} + \zeta_{it} $$ とし、\(\rho_i\)と\(\log(Var(\zeta_{it}))\)が個人レベルのランダム効果でそれぞれ正規分布すると考える。cross-classified DSEMならばランダム効果を時変させることもできる。ただし、1時点めの予測子となる\(\epsilon_{i,0}\)が指標のない潜在変数となるというのが欠点で、今度はある程度の時点数がないと困ることになる。
- 今回ご提案するRSEM(残差SEM)。2時点目以降のみ$$ \epsilon_{it} = \rho_t \epsilon_{i, t-1} + \zeta_{it} $$ として $$ \zeta_{it} \sim N(0, \theta_t)$$ とする[1時点目は\(\epsilon_{i,1} = \zeta_{i,1}\)ってことかなあ]。これなら、パラメータの部分的な時変を表現したりできるし、時点数は少なくても大丈夫だし、DSEMやRDSEMと違ってカイ二乗検定やBayesian PPPが使える。ただし、自己回帰構造は個人間で等質。
[言いたいことはわかるけど、Mplusコードのイメージがつきにくいので、示されている事例をメモしておく。7時点の線形成長モデルで誤差がAR(1)に従っているというモデルをRSEMで書くとこうなるそうである。行を個人として
MODEL:
y1-y7;
i s | y1@-3 y2@-2 y3@-1 y4@0 y5@1 y6@2 y7@3;
i s; i WITH s; [i s];
y2^-y7^ PON y1^-y6^;なるほどね。切片、傾きは個人間で変動し時変しない。誤差分散と誤差の自己回帰係数は個人間で変動せず時変する。
試しにこれをRDSEMで書いてみると、たぶんこうだと思う。行を個人x時点として
VARIABLE:
...
CLUSTER = subject;
LAGGED = y(1);
ANALYSIS:
TYPE = TWOLEVEL RANDOM:
MODEL:
%WITHIN%
sx | y on time;
sy | y^ ON y^1;
logv | y;
%BETWEEN%
y sx sy logv WITH y sx sy logv ; 今度は、切片、傾き、誤差分散、誤差の自己回帰係数のいずれも、個人間で変動し時変しない。どれかを時変させるとなるとTYPE = CROSSCLASSIFIED RANDOMで書かなければならなくなり、死ぬほどめんどくさい話になるわけだ]
3.2 ランダム切片クロスラグパネルモデル(RI-CLPM)
まずはランダム切片自己回帰(RI-AR)モデルから。$$ Y_t = \nu_t + \eta + \hat{Y}_t$$ $$ \hat{Y}_t = \rho_t \hat{Y}_{t-1} + \epsilon_t$$ [添字が省略されているが、個々の個人についての表現を意図していると思う。意図としては $$ Y_{i,t} = \nu_t + \eta_i + \hat{Y}_{i,t}$$ $$ \hat{Y}_{i,t} = \rho_t \hat{Y}_{i,t-1} + \epsilon_{i,t}$$ であろう]
これを2変量に拡張したのがランダム切片クロスラグパネルモデル(RI-CLPM)である。$$ Y_t = \nu_{1,t} + \eta_Y + \hat{Y}_t$$ $$ Z_t = \nu_{2,t} + \eta_Z + \hat{Z}_t$$ $$ \hat{Y}_t = \beta_{1,t} \hat{Y}_{t-1} + \beta_{2,t} \hat{Z}_{t-1} + \epsilon_{1,t}$$ $$ \hat{Z}_t = \beta_{3,t} \hat{Y}_{t-1} + \beta_{4,t} \hat{Z}_{t-1} + \epsilon_{2,t}$$ [なるほど、個人ごとにラグ1のVARモデルを考えるのね]
\(Var(\epsilon_{j,t}), Cov(\epsilon_{1,t}, \epsilon_{2,t})\)はなんなら時変させてもよい。
カテゴリカル従属変数の場合、残差分散の指定には4つの方略があって…[略]
カテゴリカル従属変数でベイズ推定の場合には他にも注意点があって…[略。推定上の理由で閾値は時変させるべきなのだそうな。へー]
3.3 ARMAモデル
まずはARMA(1,1)モデルから。$$ Y_t= \alpha Y_{t-1} + \beta \epsilon_{t-1} + \epsilon_t$$ 動的ランダム切片ARMAモデル(D-RI-ARMAモデル)を次のように定義する。\(t=1\)のとき $$ Y_{i1} = \alpha_1 + \eta_i + \epsilon_{i1}$$ \(t > 1\)のとき $$ Y_{it} = \alpha_t + \eta_i + \rho_t Y_{i,t-1} + \beta_t \epsilon_{i,t-1} + \epsilon_{it}$$ ただし\(\epsilon_{it} \sim N(0, \theta_t), \eta \sim N(0, \psi)\)とする。[\(\alpha_t\)が時点効果、\(\eta_i\)が個人効果で、個人レベルでARMA(1,1)誤差を持つがAR係数もMA係数も撹乱項分散も時変する、ってわけね]
D-RI-ARMAモデルは、パラメータ数からみて4時点以下では識別できないし、5時点を超えてもSEが大きくなりすぎるので制約が必要。7から50時点くらいが実用的である。
[負荷の制約の付け方とかなんとか。略]
RI-CLPMにしてもD-RI-ARMAにしても、(ふつうのSEMとはちがって)\(T\)以降の時点について推論したいという側面がある。そうした縦断モデリングのためには十分な時点数による縦断的証拠が必要であって、\(N\)が大きくて適合度が高ければよいということにはならない。
さて、ARMAモデルは測定誤差つきARモデル(MEARモデル)と等価である。
[Granger & Morris(1976 JASA)というのがreferされている。おまえは何を言っているんだという感じだが、おそらくこういうことであろう。いま、$$f_t = \phi f_{t-1} + \epsilon_t$$ という潜在的なAR(1)時系列があり、観察時系列はそれに定数とノイズを足した $$Y_t = \mu + f_t + \xi_t$$ だとする。代入すると $$Y_t = \mu + \phi f_{t-1} + \epsilon_t + \xi_t$$ 潜在時系列の前期の値を観察から逆算すると $$f_{t-1} = Y_{t-1} – \mu – \xi_{t-1}$$ これを代入して $$Y_t = \mu + \phi (Y_{t-1} – \mu – \xi_{t-1}) + \epsilon_t + \xi_t$$ $$Y_t – \phi Y_{t-1} = (1-\phi)\mu + \xi_t – \phi \xi_{t-1} + \epsilon_t$$ なので、ARMA(1,1)の形になっている。ARMAモデルはふつう係数を解釈しないが、この形式では解釈できるわけだ]
D-RI-ARMAモデルの係数を時間不変にすると、2レベルのDSEM-MEARモデルと似てくるが、以下の点が異なる。
- DSEM-MEARがARMAと等価なのはMAパラメータが負のときだけである。
- DSEM-MEARはlong formatだがD-RI-ARMAはwide formatだ[前者の行は個体x時点だが後者の行は個体だ]。
- \(t=1\)の扱いが異なる。
- ランダム切片の扱いが異なる。DSEM-ARMAではARMAが残差変数について定義され、そこにはランダム切片がない。
最後の点に関して、DSEM-ARMAやRI-CLPMと同様に、ARMAのなかにランダム切片を含まないモデルを考えよう。$$ Y_{it} = \alpha_t + \eta_i + \epsilon_{it}$$ $$ \epsilon_{i1} = \hat{\epsilon}_{i1}$$ $$ \epsilon_{i2} = \rho_2 \epsilon_{i1} + \hat{\epsilon}_{i2}$$ $$ \epsilon_{it} = \rho_t \epsilon_{i,t-1} + \beta_t \hat{\epsilon}_{i,t-1} + \hat{\epsilon}_{it} \ \ for \ \ t > 2$$
これをRI-ARMAモデルと呼ぶことにしよう。RI-ARMAは新しいMplus言語でも直接は書けず、いったん\(\epsilon_{it}\)を潜在変数として書いてやる必要がある。
[はあ？どういうこと？ と思ったら、たとえば7時点ならこう書けとのことである。
MODEL:
i by y1-y7@1;
f1 BY y1@1; y1@0;
f2 BY y2@1; y2@0;
f3 BY y3@1; y3@0;
f4 BY y4@1; y4@0;
f5 BY y5@1; y5@0;
f6 BY y6@1; y6@0;
f7 BY y7@1; y7@0;
f1 WITH i@0;
f2-f7 PON f1-f6;
f3-f7 PON f2^-f6^; うわあ。理屈はわかるけど、こんなモデルを組むぐらいなら諦めた方がマシだなあ…]
同様にRI-MEARモデルというのも考えよう。$$ Y_{it} = \alpha_t + \eta_i + \hat{Y}_{it}$$ $$ \hat{Y}_{it} = f_{it} + e_{it} $$ $$ f_{i1} = \xi_{i1} $$ $$ f_{it} = \rho_t f_{i,t-1} + \xi_{it} \ \ for \ \ t > 1$$ 識別にはいろいろと制約が必要で[…面倒くさい話がひとくさり続くがスキップ…]、制約すると結局RI-ARMAと等価になる。どちちがいいかというと…[略]。
RI-MEARモデルのMplusでの書き方は2つあって…[略]。
[ここから、カテゴリカルデータの場合の注意点、シミュレーション、2変量モデルへの拡張。なんか面倒くさくなってスキップ。必要になったら読もうと言い訳し…]
3.4 残差変数を予測子とするモデル
たとえば因子分析モデル$$ Y_p = \nu_p + \lambda_p \eta + \epsilon_p, \ \ \epsilon_p \sim N(0, \theta_p), \ \eta \sim N(0,1)$$ で、因子\(\eta\)と残差\(\epsilon_p\)を使ってアウトカム\(Z\)を予測する$$ Z = \alpha + \beta_0 \eta + \sum_{p=1}^{P} \beta_p \epsilon_p + \zeta$$というモデルを組みたくなることがあるだろう。Mplusでは Z ON Y1^ - YP^と書けば良い。
[… そんなモデルを組みたくなることがあるのか？と呆気にとられたが、考えてみると、\(P\)項目への調査回答の背後に共通手法因子\(\eta\)があると考えられるとき、\(\eta\)を抜いた独自因子だけでアウトカムを説明したい、ということならあるね。いっぽう、\(\eta\)コミでいいんなら分解せずに\(Y_p\)を使えばイイじゃんと思うのだが、著者らいわく、いや\(Y_p\)ではなく\(\epsilon_p, \eta\)を使いたいこともあるだろうとのこと。ここはよくわからなかった]
最後の式を$$ Z = \alpha + \beta_0 \eta + \hat{Z}, \ \ \hat{Z} = \sum_p \beta_p \epsilon_p + \zeta$$ と分解して, Z^ ON Y1^ - YP^と書いてもよさそうなものだが、この2つの書き方は以下の点で異なる。
- \(Z\)が他の式で予測子になっているとき、モデルが異なる。さっきのモデルでは\(Z\)の残差は\(\zeta\)だったが今度は\(\hat{Z}\)だから。
- 現在Mplusでは、Bayes推定量は残差変数が互いに回帰しているときにしか使えない。なので後者のモデルしかBayes推定できない。(MLは両方いける)
3.5 不完全な分散共分散ブロックのベイズ推定
分散共分散行列はふつうは逆ウィシャート分布が共役事前分布になるんだけど、ブロック対角でないときには共役分布がないので、MplusはMHアルゴリズムでなんとかしようとする(ALGO = GIBBS(RW); )。ところが複雑なモデルではうまく収束しない。
たとえば\(Z\)が\(Y_1, Y_2\)と相関しているとしよう。おそらく\(Y_1\)と\(Y_2\)も相関しているが、それを推定できないこともある。こういうときはMHに頼らざるを得ない。
この問題はRSEMで解決できる。Z^ on Y1^ Y2^と指定すればよい。
[ここから短い事例を2つ紹介…
残念ながら、上記の説明は知識不足により全く理解できなかった。悔しいので、2つの事例のうち2つめをほぼ逐語訳した]
2つめの事例は成長モデルである。$$ Y_{it} = I_i + S_i t + \epsilon_{it}$$ というようなモデルについて、隣接する観察のあいだについてすべての相関を含めることがよくある( \(Cov(\epsilon_{i1}, \epsilon_{i2}), Cov(\epsilon_{i2}, \epsilon_{i3}), \ldots\) )。こうしたモデルは、技術的にみると自己回帰モデルではないのだが、隣り合う観察の間の相関がもっとも高いから、自己相関の大部分はこれらのパラメータを含めることでモデル化できるだろう。その結果[隣接時点の残差相関をモデル化した結果、ということかな]、共分散行列においては対角要素の隣にあるパラメータだけが非ゼロとなる。つまり分散共分散行列はブロック対角でなくなる。こうしたモデルの推定にはMHアルゴリズムが必要になり、アルゴリズムが非効率なので推定できなくなることも多い。いっぽう残差の完全な分散共分散行列を推定できないことはあきらかである(成長モデルが識別できなくなってしまう)。
この問題はRSEMによって解決できる。隣接相関モデルを完全な自己回帰モデルに置き換えるのである。2つのモデルが同一でないことはあきらかだが(自己回帰モデルの場合は非隣接の残差間でも相関がゼロでない)、実務的には問題を完全に解決できる。
[なんとなくわかってきた… たとえば上の成長モデルで、\(\rho = Cov(\epsilon_{i1}, \epsilon_{i2}) = Cov(\epsilon_{i2}, \epsilon_{i3}) = \ldots\)をモデルのパラメータとして含めることもできるが、分散共分散行列の共役事前分布としてうまい分布が見当たらないので、推定が難しくなる。むしろ$$ \epsilon_{it} = \rho \epsilon_{i,t-1} + \xi_{it}, \ \ \xi_{it} \sim N(0, \theta_t) $$ としてくれたほうがありがたい。っていうことかなあ。パラメータの数はむしろ増えているような気がするのだが…]
3.6 残差の分散共分散を制約しないベイズ推定
未知の残差相関を発見する方法として、我らが提案しているBSEMや、ML推定量における修正指標がある。
BSEMの場合、残差分散共分散行列について非常に制約的な事前分布を指定するわけだが、その指定が難しい。そこでRSEMを使うという手がある。たとえば因子分析モデルで、$$ Y_p = \nu_p + \lambda_p \eta + \hat{Y}_p, \ \ \hat{Y}_p \sim N(0, 1)$$ とするところ、$$ \hat{Y} = \sum_{i=1}^{p-1} \beta_{pi} \hat{Y}_{i} + \zeta_p, \ \ \beta_{pi} \sim N(0, \sigma)$$ とし、\(\sigma\)を小さな値にしておくのである。
[なるほど… これは便利かも。BSEMの場合は共分散につられて分散も大きくなっちゃうしね]
—–
… 全部で84p, のこりの本文は13p。なんだか疲れちゃったのでこのへんで打ち止めにしておく。
なお残りの目次は:
3.7 条件付き依存性のある潜在クラス分析 [読んでないけど、LCAで残差同士を回帰するらしい]
3.8 Pearson事後予測P値 [新しく定義したPPPの紹介。ややこしそう…]
3.9 カテゴリカル変数の混合RI-ARモデル [ややこしい事例の紹介]
3.10 RDSEMにおける構造残差モデル [残差にラグ構造と同時構造の両方をいれるという事例]
4. 結論 [ここはたいしたこと書いてない]