« 読了：Asparouhov, Hamaker, & Muthen (2017) ものどもひれ伏せよ、これが動的潜在クラス分析だ | メイン | Rのdplyrパッケージでプログラミングするときの注意点 »
2017年8月19日 (土)
先日リリースされたMplus 8では、時系列モデルの機能が大幅拡充されている。ただでさえ目も眩むほどに多機能なのに、MODELセクションにラグつき回帰を表す記号なんかが追加されて、もうえらいことになっている。勘弁してください、Muthen導師...
Asparouhov, T., Hamaker, E.L., Muthen, B. (2017) Dynamic Structural Equation Models. Mplus technical paper.
新機能についてのテクニカル・ペーパー。実戦投入前の儀式として読んだ。
新機能が想定しているのは、たくさんの人をたくさんの時点で測定しました、というようなintensiveな縦断データ(ILDデータ)。ambulatory assessmentsとか[なんて訳すんだろう。移動測定？]、日記データとか、ecological memontary assesmenentデータとか[これも定訳があるのかどうかわかんないけど、患者の電子日記みたいなのだと思う]、経験サンプリングとか、そういうのである。
いわく...
社会科学における縦断分析の手法として成長モデリングが挙げられる。観察変数なり潜在変数なりを時間の関数としてモデル化し、そのパラメータを個人レベルのランダム効果とみなすわけである。でも一人当たり時点数はふつう10以下。データが大きいと計算が無理になるし、時点数が大きいと単純な関数を当てはめるのは難しくなる(スプラインをつかってもいいけど補外が難しい)。
そこで動的構造方程式モデル(DSEM)を提案しよう。個人をクラスタとする2レベルモデルとして成長モデルを組み、そこにラグつきの回帰とかをいれる。Molenaar(1985)とかの動的構造モデルを2レベル化するといってもよい。[なるほど...さすがはMuthen導師、簡潔にして要を得た説明だ。しかし導師よ、Molenaar(1985)って引用文献に載ってないんですけど。Psychometrikaの動的因子分析の論文ですかね]
わがDSEMフレームワーク、それは次の４つのモデリング技術の結合である。
- マルチレベル・モデリング。すなわち、個人ごとの効果がもたらす相関に基づくモデリング。
- 時系列モデリング。すなわち、観察の近接性がもたらす相関に基づくモデリング。
- 構造方程式モデリング。すなわち、変数間の相関に基づくモデリング。
- そして時変効果モデリング(EVEM)。すなわち、進化のステージが同じであることがもたらす相関に基づくモデリング。
[...か、かっこいい...導師...濡れちゃいます...(雨で)]
能書きはこのくらいにして、モデルの話。
まずはクロス分類モデル(個人のランダム効果と時点のランダム効果を併せたモデル)から。個人$i$、時点$t (=1,2,\ldots,T_i)$ における観察ベクトルを$Y_{it}$とする。これを分解する：
$Y_{it} = Y_{1,it} + Y_{2,i} + Y_{3,t}$
順に、個人$i$時点$t$からの偏差, 個人の寄与, 時点の寄与。
第2項, 第3項のモデルがbetweenレベルのモデルとなる。まず第2項は
$Y_{2,i} = \nu_2 + \Lambda_2 \eta_{2,i} + K_2 X_{2,i} + \epsilon_{2,i}$
$\eta_{2,i} = \alpha_2 + B_2 \eta_{2,i} + \Gamma_2 X_{2,i} + \xi_{2,i}$
上が測定方程式で下が構造方程式。$X_{2,i}$は個人レベルの時間不変な共変量、$\eta_{2,i}$は個人レベルの時間不変な潜在変数、$\epsilon_{2,i}$と$\xi_{2,i}$は平均ゼロの残差である。同様に第3項は
$Y_{3,i} = \nu_3 + \Lambda_3 \eta_{3,i} + K_3 X_{3,i} + \epsilon_{3,i}$
$\eta_{3,i} = \alpha_3 + B_2 \eta_{3,i} + \Gamma_3 X_{3,i} + \xi_{3,i}$
さて、第1項のモデル、すなわちwithinレベルのモデル。いよいよ時系列が入ります。
測定方程式は
$Y_{1,it}$
$= \nu_1$
$+ \sum_{l=0}^L \Lambda_{1,l} \eta_{1,i,t-l}$
$+ \sum_{l=0}^L R_l Y_{1,i,t-l}$
$+ \sum_{l=0}^L K_{1,l} X_{1,i,t-l}$
$+ \epsilon_{1,it}$
構造方程式は
$\eta_{1,it}$
$= \alpha_1$
$+ \sum_{l=0}^L B_{1,l} \eta_{1,i,t-l}$
$+ \sum_{l=0}^L Q_l Y_{1,i,t-l}$
$+ \sum_{l=0}^L \Lambda_{1,l} X_{1,i,t-l}$
$+ \xi_{1,it}$
派手になってまいりました。順に、切片、潜在変数、観察従属変数、共変量、残差である。例によって、$X_{1,i,t}$は個人x時点の共変量、$\eta_{1,it}$は個人x時点の潜在変数。
この定式化だと、withinレベルのモデルにおけるパラメータはランダム項ではないが、ランダム項にしてもいい。つまり、たとえば$\Lambda_{1,l}$を$ \Lambda_{1,lit}$としてもいいし、$R_l$を$R_{lit}$としてもよい。この場合、withinレベルのパラメータもまた個人と時点に分解する。つまり、いまパラメータ$s$をランダムにしたい場合は
$s = s_{2,i} + s_{3,t}$
と考えるわけである。
同様に残差分散$Var(\epsilon_{1,it}), Var( \xi_{1,it})$もランダムにしてよくて...[めんどくさいので略]
以上の定式化は、時系列的特徴は自己相関しか取り入れていないように見えるけど、潜在変数があるので意外に柔軟である。たとえばARMA(1,1)は
$Y_t = \mu + a Y_{t-1} + \eta_t + b\eta_{t-1}$
としてして表現できる。
上のモデルだと、時点0以前のデータが必要になるけど、そこはその、MCMCのburninのあいだに適当な事前分布を決めてですね... [とかなんとか。なんだかわからんが、とにかくどうにかなる由]
[他にもいくつか注意書きがあったけど、省略]
推定に当たっては、モデルをブロックにわけて... [省略]
モデルの適合度としてはDICを使う。すべての従属変数が連続だとする。モデルのパラメータを$\theta$, すべての観察変数を$Y$として、デビアンスは
$D(\theta) = -2 \log p(Y|\theta)$
MCMCを通じて得たデビアンスの平均を$\bar{D}$, MCMCを通じて得たモデルパラメータの平均を$\bar{\theta}$とする。有効パラメータ数を次のように定義する：
$p_D = \bar{D} - D(\bar{\theta})$
で、DICとは
$DIC = p_D + \bar{D}$
DSEMモデルにおけるDICを厳密に定義すると... [略]。
なお、DICにはいろいろ注意すべき点がある。
- 潜在変数をパラメータとみているかどうかで値が変わってくる。たとえば1因子のモデルで、因子をパラメータとしてみるならば、$p(Y|\theta)$はその因子の下で条件づけられた尤度だ(そこではその因子の指標はその因子の下で互いに独立である)。いっぽう、因子をパラメータとみない場合、$p(Y|\theta)$は因子に条件づけられた尤度ではなく、モデルに基づく分間共分散行列を使って求める尤度だ(そこでは指標は独立でない)。というわけで、潜在変数のあるモデルではDICの定義に気を付けないといけない。
たとえば、モデル1は1因子2指標のモデルで因子をパラメータとみている。モデル2は同じモデルだが、因子の分散を0に固定している(つまり指標が独立だというモデルだ)。モデル3は因子なし、指標が相関しているというモデル。さて、モデル1とモデル2のDICを比較するのはOK。モデル2とモデル3のDICを比較するのもOK。でも、モデル1とモデル3のDICを比較するのはNG。[←えええ...] - 潜在変数をパラメータとみていると、パラメータの数が増え、パラメータ推定は収束しているのにDICは収束していないということが起きる。シードを変えてMCMCしなおして確認すべし。
モデルの評価のためには、標本統計量とモデルによるその推定量を比べるという手もある。たとえば、2レベルDSEMで従属変数がひとつ($Y$)のとき、個人$i$における$Y$の標本平均を標本統計量を$\bar{Y}_{i*} = \sum ^{T_i} _t Y_{it} / T_{i}$、そのモデルによる推定量を$\mu_i$として、$R = Cor(\mu_i, \bar{Y}_{i*})$ とか$MSE = \sum_i^{N} (\mu_i - \bar{Y}_{i*})^2 / N$とかを調べたりする。
... 以上、前半のメモ。
後半はシミュレーションの紹介。扱われている問題は以下の通り。
- センタリング
- 個人内分散
- ARMA(1,1)とMEAR(1)
- AR(1), MEAR(1)に共変量を追加する方法
- 動的因子分析
- 時点が個人によって異なる場合、等間隔でない場合
- 時点ごとの効果
長いので、関心のある部分だけつまみ食い。
センタリングについて。
次の単変量AR(1)モデルを考える。
$Y_{it} = \mu_i + \phi_i (Y_{i,t-1} - \mu_i) + \xi_{it}$
$\xi_{it}$は平均0のホワイトノイズ、$\mu_i$と$\phi_i$は二変量正規とする。
DSEMならこれを一発推定できるけど、普通の2レベル回帰だと、標本平均で中心化して
$Y_{it} = \mu_i + \phi_i (Y_{i,t-1} - \bar{Y}_{i*}) + \xi_{it}$
を推定することになるわね。すると$\phi_i$にバイアスがかかる。これを動的パネルバイアス、ないしNickellのバイアスと呼ぶ。 Nickellさんの近似式によればバイアスは$-(1+\phi)/(T-1)$である。 [←へぇー]
そこでシミュレーションしてみると... DSEMで一発推定した場合、バイアスは小さく、時系列が100時点以上ならほぼゼロ。標本平均で中心化した場合、Nickellの近似式はほぼ当たっている。云々...[後略]
ARMA(1,1)とMEAR(1)について。
ARMA(1,1)モデルとは
$Y_t = \mu + \phi Y_{t-1} + \epsilon_t + \theta \epsilon_{t-1}$
だけど、これは次のモデルと等価である。
$Y_t = \mu + f_t + \xi_t$
$f_t = \phi f_{t-1} + e_t$
つまり、潜在変数$f_t$がAR(1)に従っているんだけど、その測定に誤差が乗っていると考えるわけである。これをmeasurement error AR(1)モデル、略してMEAR(1)と呼ぶことにする。云々... [あー、この話面白そう。時間がなくて読み飛ばしたけど、いつかきちんと勉強したい]
動的因子分析(DFA)について。
一般的なDFAモデルとして、直接自己回帰因子得点モデル(DAFS)とホワイトノイズ因子得点モデル(WNFS)がある。DAFSっていうのは
$Y_t = \nu + \Lambda \eta_t + \epsilon_t$
$\eta_t = \sum^L_{l=1} B_l \eta_{t-l} + \xi_t$
というモデルで、観察変数と潜在変数の間にはラグがなく、潜在変数がAR過程に従う。観察変数は結局ARMA(p,p)に従う(pは観察変数の数)。いっぽうWNFSってのは
$Y_t = \nu + \sum^L_{l=0} \Lambda_l \eta_{t-l} + \epsilon_t$
というモデルで、観察変数と潜在変数の間にラグがあり、潜在変数はホワイトノイズ。観察変数は結局MA(L)に従う。
これをハイブリッドにしたモデルを考えることもできる。すなわち
$Y_t = \nu + \sum^L_{l=0} \Lambda_l \eta_{t-l} + \epsilon_t$
$\eta_t = \sum^L_{l=1} B_l \eta_{t-l} + \xi_t$
1因子5指標、$L=1$、100人100時点でシミュレーションしてみると、[...中略...]、ちゃんとDICで正しいモデルを選択できました。云々。[途中で疲れてきて読み飛ばしたが、この節、必要になったらきちんと読もう]
...そんなこんなで、後半はほとんど読めてないんだけど、まあいいや。このtech. paperのmplusのコードは公開されているはずなので、いずれ余力ができたら勉強する、ということで。もっとも、どうせ余力なんて永遠に手に入らないのだが。
論文：データ解析(2015-) - 読了：Asparouhov, Hamaker, & Muthen (2017) ものども跪け、これがMplus8の新機能「動的SEM」だ