覚え書き: 咲き誇れ、乱れ飛べ、動的SEMのよくわからないモデルたちよ

　これは完全に、私による私のための覚書です。いや、全部そうなんですけどね。

　仕事の都合もあって、性懲りもなく時系列のマルコフ・スイッチング・モデルについて調べているのだけれど、いろいろと困惑している。手元の具体的問題というのは常になにかしら特殊性を帯びているもので、教科書的な事例がそのまま適応できることはめったにない。
　このたびは、ものすごおおおくたくさんの時系列がある、カレンダー時間は共通、レジームも共通、でもスイッチングのタイミングは時系列ごとに推定したい、という問題を抱えている。こういうのってなんていえばいいんですかね？ そんなに変わった問題ではないと思うんだけど、計量経済学の先生は頼りにならなさそうだ。

　この手の話について検索すると途端に上位に出現するのが、かのMplusの開発元、Muthen一家による動的SEM(DSEM)の論文である。パネルの縦断測定をSEMの枠組みで分析しちゃおう、計量経済学者たちよ、おまえらがやってんのは所詮は我々の分析の特殊ケース(N=1のケース)に過ぎない、そこにひざまずいて心理学者に謝れ、とおっしゃる危ない奴らである(いってません)。
　正直、勘弁してほしい。いや、Muthen先生たちに罪はないのよ。彼らの論文はわかりやすいことが多いし、Mplusは我々いたいけなユーザを泥沼に沈めるいっぽうでとんでもない問題を一発で解決してくれることもある。でも… あの先生たちの論文って往々にして、聞きなれないモデル名称が頻出するじゃないですか。さらにいえば、そもそも論文など読みたかないじゃないですか。疲れているんです、生きているだけで。

　というわけで、以下からメモ。

Asparouhov, T., Muthen, B., Hamaker, E. (2016) Latent class analysis for intensive longitudinal data, Hidden Markov processses, Regime Switching models and Dynamic Structural Equations in Mplus.

　講義のスライドである。調べたところ、U. Connecticutで2011年から隔年開催されているModern Modeling Methods (M3) conferenceの2016年の事前ワークショップとして、Muthen一家による”Mplus, Mplus 8: Dynamic SEM, and Dynamic Structural Equation Models”という一日レクチャーが開かれたようで、その際の資料らしい。朝8時半から昼までぶっつづけでMuthen導師が回帰と媒介分析についてしゃべり、午後はHamakerさんが時系列分析について90分、その後はMuthen一家若頭のAsparouhov兄貴が出てきて動的SEMについて90分を2コマ語りまくったようだ。きっと終了時には死者累々だったことだろう。
　このスライドはAsparouhovさんのパートの2コマめ。強縦断データ(intensive longitudinal data)がテーマである。
　Muthen一家はその後も精力的にセミナーを開いており、もっとよさそうな教材が一杯あるのだが、なぜそういうのではなくてわざわざこのスライドをみるかというと、見つけた範囲では、これがもっとも多くの用語が乱舞する資料だからである。まとまりがないといってもよい。ぶひー。

　ちなみに2024年のM3 conferenceは、事前ワークショップはEpskampさんという人の”Network Psychometrics in R”。基調講演は心理測定論のBorsboomさん、タイトルは”The Nature of the Measurement Game: Psychologyical Constructs as Complex Systems”(聞きてえ… これ動画とかにないのかなあ…)。プログラムをざっと眺めたけれど、結構な数の口頭セッションがパラレルで走り、ポスター発表も多い。日本人らしき方の発表を探したのだけれど、ポスターに一件フィンランドのポスドクの方がいらしたという程度であった。いやー、たまんねえな。私トランプさんにアメリカの学術研究を滅茶苦茶に破壊するよう頼んでおいたんだけど、いささか物足りない、もっと本腰を入れて頑張ってほしい。(すいません冗談です)

Mplusにおける一般的な動的構造方程式モデリング(DSEM)フレームワーク / Mplus general DSEM framework
　DSEMフレームワークとはこういうのである(正確に言うと、フレームワークのうちcross-classifiedモデルを含まない単純なバージョンは以下である)。
　個人\(i\), 時点\(t\)の観察ベクトルを\(Y_{it}\)、潜在変数ベクトルを\(\eta_{it}\), 観察共変量ベクトルを\(X_{it}\)とする。また、時間不変な観察ベクトルを\(Y_i\), 潜在変数ベクトルを\(\eta_i\), 観察共変量ベクトルを\(X_i\)とする。$$ Y_{it}= Y_{1,it} +Y_{2,i}$$ $$ Y_{1,it} = \nu_1 + \sum_{l=0}^L \Lambda_{1,l} \eta_{i,t-1} + \epsilon_{it} $$ $$ \eta_{it} = \alpha_1 + \sum_{l=0}^L B_{1,l} \eta_{i,t-l} + \Gamma_1 X_{it} + \xi_{it}$$ $$ Y_{2,i} + \nu_2 + \Lambda_2 \eta_i + \epsilon_i $$ $$ \eta_i = \alpha_2 + B_2 \eta_i + \Gamma_2 X_i + \xi_t $$ [要するに、人について時系列的に観察しているとき(人は複数だと考えるとわかりやすいが一人でもよい)、観察を時間変動部分と時間不変部分に分け(人をクラスタとする2レベルモデルと考えてもよい)、どちらのレベルでも潜在変数に対する(ラグ付き)回帰を考えて、潜在変数については外生変数つきのARモデルを考えるわけである。指標ベクトルは潜在変数ベクトルより長いかもしれない、つまり観察方程式は動的因子分析モデルになっているかもしれない]

DSEM混合モデル / DSEM Mixture Model
　レベル1(withinレベル)に水準数\(K\)のカテゴリカル潜在変数(\(S_{it}\)をいれよう。これを「状態」と呼ぶ。レベル2(betweenレベル)は状態に効かない。$$ [Y_{1, it} | S_{it} = s] = \nu_{1,s} + \sum_l \Lambda_{1,l.s}\eta_{i,t-l} + \epsilon_{it} $$ $$ [\eta_{i,t}|S_{it} = s] = \alpha_{1,s} + \sum_l B_{1,l,s} \eta_{i,t-l} + \Gamma_{1,s} x_{it} + \xi_{it}$$ \(\epsilon_{it}, \xi_{it}\)も状態ごとに考える。
　で、$$ P(S_{it}=s) = \frac{\exp(\alpha_{is})}{\sum_{s=1}^K \exp(\alpha_{is})} $$ とする。\(\alpha_{is}\)は正規ランダム効果で、\(\eta_i\)の一部である。識別のため\(\alpha_{iK} = 0\)とする。
　[最初に読んだ時、この記述が謎めていて途方に暮れていたんだけど、思い余って冒頭に戻りメモを取りながら読み直して気が付いた。\(\alpha_{is}\)って、構造方程式の右辺にある\(\alpha_{1,s}\)とは別の変数で、ベクトル\(\eta_{it}\)のひとつの要素なのだ。つまり、状態の魅力みたいな時変潜在変数があって、それが状態の確率を決めていると考えればいいよね、という話をしているのである。それを正規ランダム効果だといっているのは、潜在変数の攪乱項\(\xi_{it}\)が正規だという意味であって、これ自体が正規分布に従うとはいっていないのだと思う。Asparouhov兄貴、あってます？ というか、ここはちょっと不親切ですぜ…]

　MCMCにおいてはMHでランダム効果を更新する。バーンイン期間で\(\Sigma = Var(\alpha_{is})\)を推定し、提案分布を\(N(\alpha_{is}, \Sigma)\)とする。受容率は$$ \frac{ Prior(\hat\alpha_{is}) Likelihoood(S_{it} |\hat{\alpha}_{is}) }{ Prior(\alpha_{is}) Likelihoood(S_{it}|\alpha_{is}) }$$となる。MCMCの混合を改善する方法としては…

  • \(c\Sigma\)として、受容率が15-35%になるように\(c\)を調整する
  • バーンイン期間でクラスタごとに\(\Sigma_i = Var(\alpha_{is}|data)\)を得てクラスタごとの提案分布を\(c\Sigma_i\)とする
  • \(\alpha_{is}\)を他のランダム効果から分離する

[いささか唐突なコメントで、いまこの話をする理由がわからない。ま、講義用のスライドだから仕方ないな]

マルチレベル混合モデル / Bayes Multilevel Mixture Model
　ここからはしばらく、強縦断データから離れ、ラグつき潜在変数ベクトル\(\eta_{i,t-l}\)がないDSEM混合モデルを扱う。これはMplusでいう2レベル混合モデルそのもの。
　ベイズ推定する。最尤推定もできるんだけど、その場合は\(\alpha_{is}\)の全てを推定できないので、\(\alpha_{is}\)がそれに比例するような因子を考える形で制約をかけるしかない。詳細はAsparouhov & Muthen(2008, in Hancock & Samuelsen)をみよ。

　LCAでデータがクラスタ化されてたらどうするか。ふつうは1レベルでML推定し、SEとしてはロバストSEないしサンドイッチSEを使う[このへん不勉強でいまいちわかってないんだけど、なんかそういうのありますね]。しかし問題がある:

  • クラスの分布をクラスタごとに変えることができない。[そりゃそうだわな]
  • 潜在クラス変数の測定モデルについて完全測定不変を仮定していることになる
  • これらの仮定が成り立たない場合には擬似クラスが生まれる

測定不変でないLCA / LCA with measurement non-invariance
　OK、では測定不変じゃないLCAを考えよう。クラスタ\(j\)の個人\(i\)の指標\(p\)について[説明がないけど二値変数\(U_{pij}\)があるのだろう]、クラス別の閾値パラメータ\(\tau_{pk}\)だけじゃなくて測定変動を表すパラメータ\(\tau_{pj}\)を考えて $$ P(U_{pij} = 1 | C_{ij} = k) \Phi(\tau_{pk} + \tau_{pj}) $$ とする。\(\tau_{pk}\)は固定効果で\(\tau_{pj}\)はランダム効果である(平均ゼロ)。で、クラスタの確率は$$ P(C_{ij}=k) = \frac{\exp(\alpha_j + \alpha_{jk})}{\sum_{s=1}^K \exp(\alpha_j + \alpha_{jk}) } $$ とする。\(\alpha_j\)は固定効果で\(\alpha_{jk}\)はランダム効果である(平均ゼロ)。
　[クラスタによってクラスの確率を変えているわけね。学校によって生徒のセグメントの割合が違うというような話だ。横断データの場合は正直やり過ぎじゃないかと思うけど、あとで隠れマルコフモデルの話をするための準備なのであろう]
　Mplusのコードでいうとこんな感じ。指標はu1-u8である。ANALYSIS: TYPE = mixture twolevel; としておいて、

MODEL:
%WITHIN%
%overall%
%BETWEEN%
%overall%
C#1-C#2*0.3; [C#1*0.8 C#2*0.4]; u1-u8*0.2;
%C#1$
[u1$1-u8$1*1];
%C#2$
[u1$1-u8$1*-1];
%C#3$
[u1$1-u4$1*-1]; [u1$5-u481*-1];

[これLCAじゃん？ どこに測定変動があるの？ と呆然と眺めていたんだけど、途中で気が付いた。BETWEEN = C;とは書いていないから、潜在クラスはwithinレベル変数である(学校じゃなくて生徒を分類している)。なのに、Uの閾値のクラス別指定をbetweenレベルで書いているというところがポイントなのだ。んんん？ そうだよね？]

制約のない2レベル混合モデル / Unrestricted two-level mixture model
　$$ Y_{ij} = Y_{b,j} + Y_{w,ij}$$ $$ [Y_{w,ij}|C_{ij} = k] \sim N(\mu_k, \Sigma_k) $$ $$ P(C_{ij}=k) = \frac{\exp(\alpha_j + \alpha_{jk}) }{ \sum_{s=1}^K \exp(\alpha_j + \alpha_{jk})} $$ [恥ずかしながら、ここでこの話をしている理由がよくわからない。3本目はさっきと同じ話だけど、2本目は違う話だよね？ さっきは測定変動を表すbetweenレベルパラメータ\(\tau_{pj}\)をいれた。今度はwithinレベルでクラス別に平均と分散を変えている…]

クラスタ特有の遷移確率を伴うマルチレベル潜在遷移モデル(MLTA) / Multilevel Latent Transision Analysis (MLTA) with cluster specific transition probabilities
　学校があって生徒がいて、2時点があって、生徒は各時点において2クラスに分かれている。指標はいずれの時点でも\(p\)個の二値変数である。学校による\(P(C_2|C_1)\)の変動に関心があるとする。

　まずはこういうモデルを考えよう。[元の資料ではちゃんと分母まで書いているけど、めんどくさいので分子のみメモする] $$ P(C_{1ij} = c) \propto \exp(\alpha_{1cj} + \beta_{1cj} x_{1ij})$$ $$ P(C_{2ij} = d | C_{1ij} = c) \propto \exp(\alpha_{2dj} + \gamma_{dcj} + \beta_{2dj} x_{2ij}) $$ このモデルは最尤推定できるんだけど、\(\gamma_dcj\)は実際にはクラスタ間で変動しない。実は\(\gamma_{dc}\)となる。また\(\alpha_{2j}\)と\(\alpha_{1j}\)との間で回帰したとしてもやはり2つのランダム効果があるわけで、\(P(C_2|C_1)\)は完全には変動しない。[なにいってんだかさっぱりわからんけど、要は、遷移確率をクラスタ間で変えられないってことであろう]

　今度はこういうモデルを考える。$$ P(C_{1,ij} = k_1) \propto \exp(\alpha_{jk_1})$$ $$ P(C_{2,ij} = k_2|C_{1,ij} = k_1) \propto \exp(\alpha_{j k_1 k_2}) $$ これをベイズ推定すればよい。識別のために\(\alpha_{j k_1 K} = 0\)とする。\(\alpha_{j k_1 k_2}\)を予測子に回帰してもよい。

　2つの二値潜在クラス変数を持つモデルは3つのランダム効果を持つ。[\(C_{1,ij}, C_{2,ij}\)を\(C_1, C_2\)と略記している。上のモデルでは\(\alpha_{j1} \propto \log P(C_1=1), \alpha_{j11} \propto P(C_2=1|C_1=1)\)だけど、これを書き替えている。不親切だなあ] $$ \alpha_{j1} = \log \frac{P(C_1=1)}{P(C_1=2)} $$ $$ \alpha_{j11} = \log \frac{P(C_2=1|C_1=1)}{P(C_2=2|C_1=1)}$$ $$ \alpha_{j21} = \log \frac{P(C_2=1|C_1=2)}{P(C_2=2|C_1=2)}$$

　シミュレーションするとわかるけれど、クラスタサイズは50以上でないといけない。小さいとクロス表にゼロのセルが生じて\(Var(\alpha_{j k_1 k_2})\)が過大に推定される。
　Mplusのコードはこんな感じ。指標はu11-u14, u21-u24である。クラスはc1(2), c(2)である。ANALYSIS: TYPE = mixture twolevel random;としておいて、

MODEL:
%WITHIN%
%overall%
u11-u14*1; u21-u24*1;
s1 | C2#1 on C1#1;
s2 | C2#1 on C1#2;

%BETWEEN%
%overall%
u11-u14*1; u21-u24*1;
[C1#1*0.5]; C1#1*0.05;
[s1*-0.5 s2*1]; s1-s2*0.05;

[ああ、そうか。遷移確率行列は2×2だけど、パラメータ数は3, しかし周辺分布はC1, C2の分布として表現済だから、パラメータは2つあればよくて、遷移確率行列はクラスタによって異なるから、それらはランダム傾きs1, s2として定義できるのだ。s1, s2はbetweenレベルの潜在変数で、クラス間で共通の分布を持つから(なにしろクラス間の遷移確率ですから)、%overall%に書く。コードをじっと眺めていれば理解はできるけど、こんなん自力じゃ書けないなあ]

確率パラメータ化を伴う単一レベル潜在遷移分析 / Single Level LTA with Probability Parameterization
　レベル数1のLTAの場合、遷移確率のロジット[\(s_1, s_2\)ね]はもはやランダム効果ではなく非ランダムなパラメータである。MLの場合、パラメータ化の方法は3つある(ロジット、ログリニア、確率)。ベイズ推定の場合は確率としてパラメータ化する。パラメータは(ラグ2として) \(P(C_1), P(C_2|C_1), P(C_3|C_1, C_2)\)となる。MCMCの実装は容易である。事前分布は共役でディリクレとすればよい。
　Mplusのコードはこんな感じ…[コードにMODEL POPULATIONブロックのみあってMODELブロックがない。モンテカルロ・シミュレーションの時ってMODELブロックを省略できるのだろうか？ Mplus 8ユーザーズガイドの12章のコード例を全部みたけどMODELブロックがない例なんてない。きっとMODELブロックは必要で、c2#1 on c1#1;というようなのを書くのだろう。C2については1,2, C1については1,2,3として6本かな]

(N=1の縦断データのモデル)

　話を縦断データに戻す。まずはN=1で考える。ここからは3つのモデルについて考える。

隠れマルコフモデル(HMM) / Single Level Hidden Markov Models
　\(p=P(S_t)\)はモデルパラメータでないことに注意。[…]

　[というわけでMplusのコード例が示されているのだけれど、これが全然わかんないのである。MONTECARLOブロックでgenclasses = c(2 &1); classes = c(2 &1); と書き、MODEL POPULATIONブロックではc#1 on c&1#1*0.90;という風に書いている。
　もともとMplusでは、隠れマルコフモデルは行を測定個体とし時点を横にとったwide formatのデータで時点の数だけ潜在クラスを入れる方法で実装していた(Mplus 8ユーザーズガイドのExample 8.12)。さらにMplus 8では、行を測定個体x時点とするlong formatのデータを与えて測定個体をクラスタとし、VARIABLEブロックでLAGGED = y (1);とすればMODELブロックではy&1が\(y_{t-1}\)を表すという文法も実装された(Example 9.32)。しかし、潜在クラスについてもこういう書き方ができるとしている資料は、この講義スライドのほかには見当たらない。コード例を一字一句書き写して走らせてみたけど、Mplus 8.11でではやはりエラーとなった。こんな文法は、2016年5月当時に想定されていた将来の機能であって、製品版には現在に至るまで実装されていないのではないだろうか？
　この後の事例もこの文法で書かれているので、以下Mplusコードのメモは省略する。がっくり。
　なお、試しにChatGPT君に聞いてみたところ、「Mplusにおいてclasses = c(2 &1);とは2水準の潜在クラスでクラス1を参照クラスにしますという意味です」とのことであった。いや、それはないでしょう… 平然とそれっぽいことを云ってくるなあ]

2変量のマルコフ・スイッチング自己回帰モデル / Bivariate MSAR
　こんどは混合自己回帰モデルとHMMを合体させ、マルコフ・スイッチング自己回帰モデルとする。[最後の攪乱項は\(\epsilon_{1it}, \epsilon_{2it}\)となっていたけれど、ここだけ\(i\)がついているのは変なので誤植と考えて勝手に直した。あるいは、\(Y_{1t}, Y_{2t}\)のほうが\(Y_{1it}, Y_{2it}\)の誤植なのかも？]$$ Y_{1t} = \alpha_{1,S_t} + \beta_{1,S_t} Y_{1,t-1} + \beta_{2,S_t} Y_{2,t-1} + \epsilon_{1t}$$ $$ Y_{2t} = \alpha_{2,S_t} + \beta_{3,S_t} Y_{1,t-1} + \beta_{4,S_t} Y_{2,t-1} + \epsilon_{2t}$$

マルコフ・スイッチング・カルマン・フィルタ (MSKF) / Markov Switching Kalman Filter(MSKF)
　こんどは混合カルマンフィルタ[なにそれ]とHMMを合体させ、マルコフ・スイッチング・カルマン・フィルタモデルとする。
　\(\eta_t\)が指標\(Y_{jt}, j=1,2,3\)を持っているとする。\(S_t\)を2状態カテゴリカル潜在変数とし隠れマルコフモデルを組む。で、因子についてはマルコフ・スイッチング自己回帰モデル MSAR(2) を組む。$$ Y_{jt} = \nu_j + \lambda_j \eta_t + \epsilon_{jt} $$ $$ \eta_t = \alpha_{S_t} + \beta_{1,S_t} \eta_{t-1} + \beta_{2,S_t} \eta_{t-2} + \xi_t $$ 識別のため\(\alpha_1=0, \lambda_1=1\)とする。
　[コード例と結果。上記の理由により省略…]
　モンテカルロ・シミュレーションの結果小さなバイアスがみられた。[その考察。パス]

　[なるほど、要するに動的因子モデルで因子がAR(2)になっていて指標にはラグがない奴で、構造方程式だけレジーム・スイッチするわけね。観察方程式は時変しない。それにしても、これをカルマン・フィルタというのはなぜ？ カルマン・フィルタってモデルの名前じゃなくて推定方法の名前でしょうに]

(強縦断データのモデル)

動的潜在クラス分析 / Dynamic Latent Class Analysis
200人、時点数100, クラスの指標4つ、2クラスとして、$$ P(U_{pit} = 1|S_{it} = k) = \Phi(\tau_{kp} + \tau_{ip}) $$ $$ \tau_{ip} \sim N(0, \sigma_{ip})$$ $$ P(S_{it} = 1 | S_{i,t-1} = 1) = frac{ \exp(\alpha_{i1} }{ 1+\exp(\alpha_{i1}) } $$ $$ P(S_{it} = 2 | S_{i,t-1} = 2) = frac{ \exp(\alpha_{21} }{ 1+\exp(\alpha_{i2}) } $$ $$ \alpha_{ij} \sim N(\alpha_j, \sigma_j) $$ とする。[Mplusコード例。省略するけど、100時点もあったらwide formatでは到底無理だな… どうすればいいんだろうか]

　[こうなるとだいたい想像がついてくるので、以下はメモ省略]

マルチレベル・マルコフ・スイッチング自己回帰モデル (MMSAR) / Multilevel Markov Switching Autoregressive Models (MMSAR)
マルチレベル・マルコフ・スイッチング・カルマン・フィルタ・モデル (MMSKF) / Multilevel Markov Switching Kalman Filter Models (MMSKF)

課題

  • クラス数の決定。時系列を無視して普通の方法でやる。[そうしろってこと？ 現状そうするしかないってこと？]
  • 初期値の指定。もっともベイズではたいしたことではない。
  • モデルの比較。DICとか検定とかいろいろ準備中。
  • (これは時系列に限ったことじゃないけど) MLではたくさんの解が生じることがある、ベイズではどうか。
  • ラベル・スイッチング。

————–
　やれやれ、疲れた…
　それにしても、実際には動作しないコード例が出てきたのには参った。Mplusは8.1以降の仕様変更の解説がいろんなPDFに分散していてわかりにくい(いちおうここにまとめてはみたけれど)。そろそろMplus 9をリリースしてくださらないだろうか。
　潜在クラス変数についてラグ表記ができないとなると、強縦断データのマルコフ・スイッチングは到底書けそうにない。ひょっとして、MODELブロックでC&1という表記は可能で、VARIABLEブロックでの指定の仕方がこの資料とは異なるのだろうか？ もしやLAGGED = C (1);って書いていいの？ やっぱり、もっと新しい講義資料をみたほうがよかったな…