« 読了: Im, et al.(2016) cross-classifiedマルチレベルデータをただの階層データとみなしたときの弊害 (多群CFAで測定不変性を検証する編) | メイン | 覚え書き：cross-classifiedなマルチレベルモデル »
2019年1月14日 (月)
Chow, S., Ho, M.R., Hamaker, E., Dolan, C.V. (2010) Equivalence and Differences Between Structural Equation Modeling and State-Space Modeling Techniques. Structural Equation Modeling, 17, 303-332.
仕事の都合で泣く泣く読んだ奴。いま気がついたのだが、著者に動的因子分析のHamakerさん, 状態空間モデルのDolanさんがはいっている。先に著者名をちゃんとみておこうよ、俺...
SEM(構造方程式モデリング)と状態空間モデリングのどこが同じでどこがちがうか説明しましょう、という論文。
壮大な問題設定なのでびびってしまうが、あとでみるように、この論文がSEMと呼んでいるのはいわゆるSEMよりちょっと狭くて、測定変数間にパスがなく、因子に外生的な説明変数がないようなやつのことである(つまり、CFAの因子間にパス引いた奴)。SEM関連の本は全部オフィスの本棚に並べておるのでいまわからんが、普通そういうモデルはなんて呼ぶんですかね。多重指標モデル？
イントロは省略して...
まずSEMについて説明しよう。[原文の添字を簡略化してメモする]
$y_i = \tau + \Lambda \eta_i + \epsilon_i$
$\eta_i = \alpha + B \eta_i + \zeta_i$
ただし$y_i$, $\tau$, $\epsilon_i$は長さ$p$のベクトル, $\eta_i$, $\alpha$, $\zeta_i$は長さ$w$のベクトル。$\epsilon_i$の共分散行列を$\Theta$, $\zeta_i$の共分散行列を$\Psi$とする。
母共分散行列と母平均は
$\Sigma = \Lambda(I-B)^{-1} \Psi (I-B)^{-1'} \Lambda_{'} + \Theta$
$\mu = \tau + \Lambda (I-B)^{-1} \alpha$
SEMの場合、複数時点で測定しているデータについては、ふつうは$y_i$を二次元にするのではなくて長さを伸ばす。$T$時点の縦断データだったら、$T$時点がすべて$y_i$, $\eta_i$に叩き込まれるわけである。
パラメータ推定にあたっては、ふつはMVNを仮定して尤度を求める。この場合は標本共分散行列$S$と標本平均ベクトル$m$があればよい。でもローデータから尤度を求めることもできる(以下RMLと呼ぶ。FIMLともいう)。時点数が個人数を超えている反復測定データなんかの場合はRMLを使わざるを得ない。
RMLでは、個々の個人のデータベクトルの対数尤度を足し上げた
$LL(\theta) = (1/2)\sum_i [-p \log(2\pi) - \log|\Sigma| - (y_i - \mu)^{'} \Sigma^{-1} (y_i - \mu)]$
を最大化する。MVNを仮定すると、この式はこう書き換えられる。
$F(\theta) = (1/2)[\log|\Sigma|+tr(S\Sigma^{-1}) - \log|S| - p] + (1/2)[(m-\mu)^{'} \Sigma^{-1} (m-\mu)]$
第1項はモデルの対数尤度、第1項は完全に飽和したモデル(つまり$\Sigma = S, \mu=m$であるモデル)の対数尤度である。これに$-2 (N-1)$を掛けるとカイ二乗統計量になる。[途中から写経モードに... まあ信じますよ先生]
今度は状態空間モデルについて。
$y_{it} = \tau + \Lambda \eta_{it} + \epsilon_{it}$
$\eta_{it} = \alpha + B \eta_{i,t-1} + \zeta_{it}$
標準的な状態空間表記ではすべてのパラメータを時変とみるのだが、ここでは不変とみる。また外生変数もないことにする。
パラメータ推定はふつうカルマン・フィルタ(KF)でやる。
まず$t=0$における状態ベクトル$\eta_{i,0|0}$とその共分散行列$P_{0|0}$を決める。たとえば前者は0のベクトル、後者は$\kappa I$($\kappa$は十分に大きな数)とか。さらに、$\epsilon_{it} \sim N(0, \Theta)$, $\zeta_{it} \sim N(0, \Psi)$と仮定する。
次に、各時点における$i$の状態ベクトルとその共分散行列を、それまでの観察によって予測する。
$\eta_{i,t|t-1} = \alpha + B \eta_{i,t-1|t-1}$
$P_{t|t-1} = B P_{t-1|t-1}B^{'} + \Psi$
以下では
$e_{i,t|t} = y_{it}-(\tau + \Lambda \eta_{i,t|t-1})$
をイノベーションと呼ぶ(one-step-ahead予測誤差ともいう)。期待値は0である。
さて、この2つの推定値を、その時点の観察を使って更新する。
$\eta_{i,t|t} = \eta_{i,t|t-1} + K_{t|t} e_{i,t|t}$
$P_{t|t} = [I-K_{t|t}\Lambda] P_{t|t-1}$
ただし
$K_{t|t} = P_{t|t-1} \Lambda^{'} [\Lambda P_{t|t-1} \Lambda^{'} + \Theta]^{-1}$
これをカルマン・ゲインという。
さて、全員について$T$時点の推定を行った暁には、イノベーション$e_{i,t|t-1}$とその共分散行列
$G_t = E[e_{i,t|t-1}e_{i,t|t-1}^{'}]$
が手に入る。後者はすべての$i$について$\Lambda P_{t|t-1}\Lambda^{'}+\Theta$である。
ここから対数尤度関数が手に入る。
$LL(\theta) = (1/2)\sum_i \sum_t [-p \log(2\pi) - \log|G_t| -(e_{i,t|t-1}^{'})G_t^{-1} (e_{i,t|t-1}) ]$
これを予測誤差分解(PED)関数という。これを最大化すればいい。
要約しよう。まずパラメータに初期値を与える。つぎにカルマンフィルタで個々人の各時点での状態を推定する(このときパラメータは定数とみなす)。次に状態推定を既知と見なしてパラメータをPED関数で推定する。このパラメータを定数とみなして新しい状態を推定する。これを繰り返す。
いよいよ本題。SEMと状態空間モデリングは、どこが同じでどこがちがうか。
ひとことでいうと次の通りである。どちらも「プールされた個人内ダイナミクス」を表現している。SEMは潜在変数間の同時的な構造的関係と、その関係の個人間の差を捉えるのに適している。状態空間モデルは、もっと複雑な個人内ダイナミクスを表現するのに適している。
もっと細かくいうと次のとおり。以下、SEMのパラメータに添字$S$, 状態空間モデルのパラメータに添字$K$をつける。
状態推定について。$B_K=0$とすれば、カルマンフィルタは回帰法による因子得点推定と同じである。また$B_K=0$で$\Psi_K$が無限大に近いとき(つまり事前情報がないとき)、カルマンフィルタはバートレット法による因子得点推定に近づく。[...中略...]
パラメータ推定について。$T=1$の場合、$B_S = B_K = 0$となり、2つのモデルは等しい。$T > 1$の場合、$p_S = T \times p_K$ならば、対数尤度関数が同じになる。さらに$B_S = B_K = 0$なら、パラメータ推定値は同じになる。
問題は$B_S, B_K$が0でない場合。特別な制約を付けない限りパラメータは変わってくる。その制約がどんなものかは、実例をみたほうが早かろう。
実例1, 横断共通因子モデル。ここでみなさんに、「状態空間モデルはある個人の反復測定データを扱う手法である」というよくある思い込みから脱して頂きます。
SEMの場合。$B$が0になるので、
$y_i = \tau + \Lambda \eta_i + \epsilon_i$
$\eta_i = \alpha + \zeta_i$
共分散行列を$\Theta = E(\epsilon_i \epsilon_i^{'}), \ \Psi = E(\zeta_i^2)$として、
$\Sigma = \Lambda \Psi \Lambda^{'} + \Theta$
$\mu = \tau + \Lambda \alpha$
$LL = (1/2) \sum_i [-p \log(2\pi) - \log |\Sigma| - (y_i - \mu)^{'} \Sigma^{-1} (y_i - \mu)]$
対数尤度関数は状態空間モデルと同一になる。
[ほんとなの？ 試してみよう。
$y_{it} = \tau + \Lambda \eta_{it} + \epsilon_{it}$
$\eta_{it} = \alpha + \zeta_{it}$
イノベーションは本来は$e_{i,t|t-1} = y_{it}-(\tau + \Lambda \eta_{i,t|t-1})$だけど、$\eta_{i,t|t-1} = \alpha$だから
$e_{i,t|t-1} = y_{it}-(\tau + \Lambda \alpha)$
$G_t = \Lambda \Psi \Lambda^{'}+\Theta$
$LL = (1/2) \sum_i \sum_t [-p \log(2\pi) - \log|G_t| -(e_{i,t|t-1}^{'})G_t^{-1} (e_{i,t|t-1}) ]$
なるほど、同じだわ]
シミュレーション。$N=200, T=1, p=4$のデータについて$w=1$の共通因子モデルを組みます。パラメータの真値を適当に決めて架空データをつくり、LISRELとmkfm6(という状態空間モデルのプログラムがある由)で分析してみる。結果はほぼおなじ。
実例その2, 単変量LDSモデル(Nが大きくてTが小さい場合)。
McArdle & Hamagamiのlatent difference scoreモデルについて考える。(ああ、読んでおいて良かった... なにがいつ役に立つかわからんねぇ...)
LDSモデルでは、観察の後ろに潜在変数$\eta_{it}$を考え、さらに「潜在変数の前期からの差」という潜在変数$\Delta \eta_{it}$を考える。
$y_{it} = \eta_{it} + e_{it}$
$\eta_{i, t+1} = \eta_{it} + \Delta \eta_{i,t+1}$
でもって、
$\Delta \eta_{it} = \eta_{si} + \beta \eta_{i,t-1}$
と考える。$\eta_{si}$は$i$さんの定数項で変化の強さを表す。$\beta$は全員共通の定数項で、自己フィードバックを表す。以上をまとめると
$\eta_{i,t+1} = \eta_{si} + (1+\beta) \eta_{it}$
$i$さんの時点0における切片項と変化の強さを、平均と偏差に分解して
$\eta_{0i} = \mu_{\eta 0} + \zeta_{\eta 0 i}$
$\eta_{si} = \mu_{\eta s} + \zeta_{\eta s i}$
$E(\zeta_{\eta 0 i}^2) = \sigma^2_{\eta 0}$, $E(\zeta_{\eta s i}^2) = \sigma^2_{\eta s}$とする。
LDSモデルは要するに成長曲線モデルだと考えてよい。ちがいは、潜在傾きをさらにモデル化しているという点。
SEMで表現してみましょう。
長さ$T$のベクトル$y_i$を考える。測定モデル
$y_i = \Lambda \eta_i + \epsilon_i$
で、潜在変数$\eta_i$を、長さ$T+1$のベクトル$[\eta_{i1}, \ldots, \eta_{iT}, \eta_{si}]'$とする。$\Lambda$は$T \times T$の単位行列の右に、すべて0の列をくっつけた奴にする。要するに、$\eta_{si}$は測定モデルには出てこないんだけど、話の都合上、潜在変数ベクトルにいれておくわけである。トリックですね。
構造モデル
$eta_i = \alpha + B \eta_i + \zeta_i$
では、$\alpha$は先頭が$\mu_{\eta 0}$、末尾が$\mu_{\eta s}$、ほかはすべて0のベクトルにする。$\zeta_i$も同様で、先頭が$\zeta_{0i}$、末尾が$\zeta_{si}$、ほかはすべて0である。
$B$はどうなるか。時点1については(つまり$B$の1行目)、すべて0でよい。時点$2,\ldots,T$については(つまり$B$の2行目から$T$行目)、対角に$1-\beta$をいれた$(T-1) \times (T-1)$の対角行列の右に、すべて0の列($\eta_{iT}$にかかる係数)、すべて1の列($\eta_{si}$にかかる列)をくっつければ良い。最後の行は$\eta_{si}$をつくるための係数だから、すべて0。ね？ 実は$\Delta \eta_{it}$なんていらないのだ。[←なるほど...]
状態空間で表現してみましょう。[原文とはちがうが、話がややこしくなるので、状態空間モデルでの潜在変数はアスタリスクをつけて$\eta^{*}$とする]
観察方程式
$y_{it} = \Lambda \eta_{it}^{*} + \epsilon_{it}$
で、潜在変数$\eta_{it}^{*}$をベクトル$[\eta_{it}, \eta_{si}]'$にする。$\eta_{si}$は時間不変だし、そもそも観察変数には出てこない奴なので、変な感じだが、これもトリックですね。$\Lambda = [1,0]'$である。
遷移方程式
$\eta_{it}^{*} = B \eta_{i,t-1}^{*} + \zeta_{it}$
はどうなるか[ここ、原文に誤植があると思うので勝手に直した]。$\zeta_{it}$はすべて0。$B$は$2 \times 2$行列で、1行目($\eta_{it}$をつくる係数)は$[1-\beta, 1]'$。2行目は、時間不変だから$[0,1]'$。
最後に、初期値$\eta^{*}_{i1|1}$を$\mu_{\eta 0}, \mu_{\eta 1}$, その共分散行列$P_{i1|1}$の対角要素を$\sigma^2_{\eta 0}, \sigma^2{\eta s}$とする。
シミュレーション。$N=400, T=10$のデータについてLDSモデルを組む。パラメータの真値を適当に決めて架空データをつくり、LISRELとOx(という状態空間モデルのプログラムがある由)で分析してみる。結果はほぼおなじ。
というわけで、LDSモデルはSEMとしても状態空間モデルとしても推定できるんだけど、後者のほうが便利。なぜなら、SEMだと(1)$T>N$のときに困るし、(2)$T$が大きいだけでえらいことになる。
実例その３，自己回帰動的因子モデル($T$が大きくて$N=1$)。
長さ$p$の観察ベクトル$y_{it}$について
$y_{it} = \Lambda \eta_{it} + u_{it}$
$\eta_{it} = A_1 \eta_{i,t-1} + A_2 \eta_{i,t-2} + \ldots + \zeta_{it}$
というのがこのモデル[原文に誤植があると思うので勝手に直した]。以下、2因子、6指標、ラグ1までのモデルを考えます。
状態空間表現の場合。観察方程式
$y_{it} = \Lambda \eta_{it} + \epsilon_{it}$
の$\eta$は長さ2のベクトル。$\Lambda$は$6 \times 2$の行列で、1列目を$[1, \lambda_{21}, \lambda_{31}, 0, 0, 0]'$, 2列目を$[0,0,0,1, \lambda_{52}, \lambda_{62}]'$とする(識別のためにひとつ1をいれている)。
遷移方程式
$\eta_{it}^{*} = B \eta_{i,t-1} + \zeta_{it}$
の$B$は$2 \times 2$の行列で、ラグ1の自己回帰係数。$\zeta_{it}$はこのVAR(1)プロセスで説明できない因子得点である。
SEMの場合。データ行列は$N \times T p$となるわけで、共分散行列をつかった最尤推定はできない。ひとつの路線はRMLを使うことだが、計算が大変である。もうひとつはブロック・トープリッツ行列を使うという路線。こっちで考えます。
[あれ、ほかにも方法があるのではなかろうか。まあいいけどさ]
$N=10, T=200$でモンテカルロシミュレーションしてみよう。[...メモ中略...] とこのように、ブロック・トープリッツ行列路線はパラメータ推定値がけっこう歪む。それに個人レベルの因子得点も推定できない。お勧めできない。
考察。
心理学のみなさんは状態空間モデルをあまり使わない。もったいない。SEMにも状態空間モデルにもそれぞれ美点があるので、使い分けましょう。
云々、云々。
論文：データ解析(2018-) - 読了：Chow, Ho, Hamaker, Dolan (2010) SEMと状態空間モデルのどこが同じでどこが違うか