自己報告多変量時系列データ の因果沼に潜ってみた 福島県立医科大学医学部 健康リスクコミュニケーション学講座 助教 竹林由武 (公認心理師/臨床心理士) 19/3/21 社会心理学会 春の方法論セミナー2019 @明治学院大学 https://ytake2.github.io/Rsite/_site/index.html

Topics • 自己報告の時系列データ • VARの基礎 – ARモデル – グランジャー因果・インパルス応答関数 – ネットワークモデルとの親和性 • VARの拡張 – 定常性を越えて: 時変VARモデル – 定間隔を越えて：連続時間モデル 2

時系列データの構造: 経験サンプリングの場合 3

Cattel’s data box 4 ESM/EMA データ (多変数、多人数、多時点) N of 1 データ (多変数、1人、1時点) 縦断調査データ (１変数、多人数、 短期間) 横断調査データ (多変数、多人数、1時点) https://www.statmodel.com/download/HamakerDSEMforPSMG.pdf

経験サンプリング 5 スマートフォン等で、 1日の数回ランダムに通知 通知が来たら質問に回答 例) 7:30〜22:30 90間隔で1回ランダムに通知 1日で10回の回答 1週間で70回の回答

ESMで得られるデータ (N of 1) • 1人の参加者に対して複数の変数を反復測定 • 1人の参加者のみ各種統計手法が適用できる 6 Subj Time X Y 1 1 4 7 1 2 2 8 1 3 3 0 1 4 7 1 1 5 8 3 1 : : : 1 : : : 1 84 6 10 相関係数= .62 X Y

ESMで得られるデータ (N of >1) • マルチレベルモデル 7 Subj Time X Y 1 1 4 7 1 : : : 1 84 3 0 2 1 8 3 2 : : : 3 84 7 1 4 1 2 6 4 : : : 4 84 6 10 集団レベル (相関の集団平均) 個人レベル (個人ごとに異なる相関) 個人レベルの情報を集団レベルの推定に活用 個人レベルの情報を有効に使いたい

ESMが生む情報 (N of 1) • 単変量 – 平均 – 変化 (傾き) – 変動性/安定性 • 多変量 – 変数間の因果関係 平均や傾き: 通常の(反復の少ない)縦断調査でも可 変動性や因果関係:測定頻度の高いESMならでは 8 平均の違い 変化の違い 変動の違い T-1 X T-1 Y T X T Y 因果関係

単変量の時系列：ARモデル 9

AR (1)モデル Yt-1 Yt Ε Yt = 各時点 (t)のYの得点, Yt-1 = Yt の1時点前のYの値, ε = 誤差 (ホワイトノイズ) Yt = c+φYt-1 +ε 現在のデータと1時点前のデータに線形の関係がある Φ= 0.1 Φ= 0.7 De Haan-Rietdijk, S., Gottman, J. M., Bergeman, C. S., & Hamaker, E. L. (2016). Get over it! A multilevel threshold autoregressive model for state-dependent affect regulation. psychometrika, 81(1), 217-241.

AR (1)モデルの適用例 • 感情慣性 (Emotional inertia): 感情指標の自己相関 11 Φ= 0.1 Φ= 0.7 １度嫌な気分になると引きずる（自己相関が高い） 1時的に嫌な気分になることがあるが、引きずらない (自己相関が低い） Suls, J., Green, P., & Hillis, S. (1998). Emotional reactivity to everyday problems, affective inertia, and neuroticism. Personality and Social Psychology Bulletin, 24(2), 127-136.

AR (1)モデルの適用例 Kuppens, P., Allen, N. B., & Sheeber, L. B. (2010). Emotional inertia and psychological maladjustment. Psychological science, 21(7), 984-991. Koval, P., Sütterlin, S., & Kuppens, P. (2016). Emotional inertia is associated with lower well- being when controlling for differences in emotional context. Frontiers in psychology, 6, 2016. Suls, J., Green, P., & Hillis, S. (1998). Emotional reactivity to everyday problems, affective inertia, and neuroticism. Personality and Social Psychology Bulletin, 24(2), 127-136. 健康な地域住民 (25-48才)48名にDRM. 神経症傾向の高い人は、ネガティブ気分の自己相関高く、 神経症傾向の低い人は、ネガティブ気分の自己相関低い 大学生80名にESM. 自尊心の低い人は、抑うつがある人は、 ネガティブ気分でもポジティブ気分でも自己相関高い 大学生約100名が気分誘導映像視聴後に気分を反復評定 (ポジティブ気分の自己相関よりも)ネガティブ気分の自己相関の 高さが抑うつ症状と関連(人生の満足度よりも)

推定コード： Rの場合 • ARモデルの推定 lm(y ~ y_lag1, data=data) • マルチレベルARモデルの推定 lmer(y ~ y_lag1+(1+y_lag1|id), data=data) 13

多変量の時系列：VARモデル 14

VARモデル 15 D T-1 D T Act T-1 Act T T T-1 例）一時点前のD得点が一時点後のACT得点を説明 個々の時系列の自己相関に加えて、 時系列間の因果(的な)関係を検討したい 自己回帰 自己回帰 Depression Activation

個々の時系列のAR（1）モデル Y1t-1 Y1t Y2t-1 Y2t ε1 ε2 変数Y1 のAR（1）プロセス 変数Y2 のAR（1）プロセス Y1t = c+ φ11 Y1t-1 +ε1 Y2t = c+ φ21 Y2t-1 +ε2

2変量のVARモデル Y1t-1 Y1t Y2t-1 Y2t ε1 ε2 Y1t = c+ φ11 Y1t-1 + φ12 Y2t-1 +ε1 Y2t =c+ φ22 Y1t-1 +φ21 Y2t-1 +ε2 残差間に相関がありうる

VARモデル Y1t = c+ φ11 Y1t-1 + φ12 Y2t-1 +ε1 Y2t = c+ φ22 Y1t-1 + φ21 Y2t-1 +ε2 Yt = C+φYt-1 +εt n*nの係数行列 n個の変数ベクトル n*1の定数ベクトル Yt = c+φYt-1 +ε ARモデルが多変量に拡張しただけ VARモデル ARモデル

Granger因果 Y1t-1 Y1t Y2t-1 ε1 ① 変数Y1 のAR1 ② 変数Y1 のAR(1)の説明変数にY2t1 を導入 Y1t-1 Y1t ε1 ①でY1を予測した場合と②でY1を予測した場合に、①よりも➁で残差（予測 誤差）が小さくなれば（予測精度が向上すれば)、因果的な関係あり

Granger因果のイメージ 20 Granger因果による時系列データの因果推定（因果フェス2015）: 尾崎 隆 https://www.slideshare.net/takashijozaki1/granger2015 あくまで Grangerの意味で 因果関係あり (因果推論において、 十分条件を満たすわけ ではない） この結果のみで、因果関係を 述べるのは避けるべき

Granger因果性の検定 • 2Fは漸近的にχ2(2)に従う • χ2(2)の95%点を比較して2Fが大きい • Y2→Y1のgranger因果性が存在しないという帰無仮説が 棄却 • Y2→Y1のgranger因果性ありという対立仮説支持 21 F = (SSR 0 - SSR 1 ) / 2 SSR 1 / (T -5) ①の残差平方和: SSR1 ②の残差平方和：SSR0 サンプルサイズ： T

Granger因果の注意点 22 ・F検定に依拠するGranger因果性検定は、 単位根過程を含むVARモデルでは適用できない ・時系列A、Bが共和分関係にある場合、 1階差VARモデル自体が定義できないので、Granger 因果性の検定がそもそも成立しない ・Granger因果はあくまで2変数間での関係のみを言及 第三の変数の影響は考慮されていない （偏相関の発想で、共変量を統制した偏Granger因果という検定法もある）

ESMデータへのVARモデルの適用例 23 Tschacher, W., Zorn, P., & Ramseyer, F. (2012). Change mechanisms of schema-centered group psychotherapy with personality disorder patients. PloS one, 7(6), e39687. パーソナリティ障害患者69名を対象にスキーマ中心感情行動療法 (Schema centered emotion-behavior therapy: 初耳)を、 集団形成(100分のセッションを週二回）で実施。合計50回程度実施。 毎セッションで記入してもらった自己評定尺度の得点にVARを適用。 Granger因果検定で 有意な関連のみを図示 患者1名の時系列データの例 “洞察”の増加が、拒絶過敏性や感情的覚醒の改善に影響する 個人にVARを適用し得られた回帰パラメータについてt検定をし集団の推定結果

VAR+ネットワークモデル • VARは、変数間の関係を総当たりで検討しているので、 推定結果がネットワークモデルで可視化できる – 観測変数→ノード – クロス回帰係数→エッジ – 自己回帰係数→ セルフフィードバック 24 AさんのVAR BさんのVAR ネットワークモデルの指標に 基づいた解釈が可能 中心性 (centrality) 密度 (density) Bringmann, L. F., Vissers, N., Wichers, M., Geschwind, N., Kuppens, P., Peeters, F., ... & Tuerlinckx, F. (2013). A network approach to psychopathology: new insights into clinical longitudinal data. PloS one, 8(4), e60188.

VAR+ネットワークモデル 25 Groen, R. N., Snippe, E., Bringmann, L. F., Simons, C. J., Hartmann, J. A., Bos, E. H., & Wichers, M. (2019). Capturing the risk of persisting depressive symptoms: A dynamic network investigation of patients' daily symptom experiences. Psychiatry Research, 271, 640-648. ネットワーク中心性 (マルチレベル)VAR うつ病患者69名にESMを使った介入を実施。 介入前から介入後6ヶ月の抑うつ尺度の評定値から、症状維持群の改善群に分類。 介入期間中のESMデータについて(マルチレベル)VAR+ネットワークモデルを適用 ネットワーク中心性の高い指標（全てが億劫に感じる）が、他の指標よリも予後を予測

VAR+ネットワークモデル 26 うつ病患者(53名）と健常成人(53名） にESMを実施し、マルチレベルVAR+ ネットワークモデルを適用。 うつ病患者群は健常群よりも、 高密度な気分ネットワーク。 特にネガティブな感情のネットワーク が、高密度。 ネットワークの群間差も検定が可能 （mlVAR::mlVARcompare） Pe, M. L., Kircanski, K., Thompson, R. J., Bringmann, L. F., Tuerlinckx, F., Mestdagh, M., ... & Kuppens, P. (2015). Emotion-network density in major depressive disorder. Clinical Psychological Science, 3(2), 292-300.

ネットワーク推定の正確性/安定性 • ネットワークの安定性評価 （bootnet package） 27 エッジのブートスラップ信頼区間 中心性のブートスラップ信頼区間 サンプルからランダムに標本抽出を繰り返して パラメータの信頼区間を構成 Epskamp, S., Borsboom, D., & Fried, E. I. (2018). Estimating psychological networks and their accuracy: a tutorial paper. Behavior Research Methods, 50(1), 195-212.

インパルス応答 • VARである変数の時系列の式の誤差項にショック (通常1標準偏差分)を与えることで、一方の変数 に動きがどうなるかをみるもの • 一方の変数に変化があった場合に、他の変数に どの程度の影響があるか、影響がどれくらい続く か、を検討する。 • 誤差項に相関がある場合には（大抵ある）、コレス キー分解で直交化 →直交化インパルス応答関数

インパルス応答分析 • 当時点である変数にある変化があった際に、数時 点後の他の変数がどの程度変化するか y2にある時点で生じた変化は、いずれ の時点後においてもy1に影響しない y2のある時点変化は、時点5後までy3 に影響を与える 直交化インパルス応答関数 ブルーの範囲はブートストラップ95%信頼区間

ESMデータへのVARモデルの適用例 • レンポンス応答関数 – インパルス応答関数で信頼区間幅が±0.01以下になる 時点を回復時点 (recovery time)とする – 感情の変化が元に戻る時間の早さ=レジリエンス 30 Yang, X., Ram, N., Gest, S. D., Lydon-Staley, D. M., Conroy, D. E., Pincus, A. L., & Molenaar, P. (2018). Socioemotional Dynamics of Emotion Regulation and Depressive Symptoms: A Person-Specific Network Approach. Complexity, 2018.

ESMデータへのVARモデルの適用例 • レンポンス応答関数 – 個々人の回復時点とうつ尺度得点に有意な関連 - ネガティブライフイベント経験が多い場合、 回復時点と抑うつ症状の関連が強い - 回復時点が大きいほど抑うつが強い 31 Yang, X., Ram, N., Gest, S. D., Lydon-Staley, D. M., Conroy, D. E., Pincus, A. L., & Molenaar, P. (2018). Socioemotional Dynamics of Emotion Regulation and Depressive Symptoms: A Person-Specific Network Approach. Complexity, 2018.

• vars – N of 1のVAR, granger因果検定, インパル応答関数 • graphicalVAR: – N of 1データのVAR+ネットワークモデル, VARのパラメータ推定に generalized lassoを使用 • pompom: – N of 1データのVAR+ネットワークモデル+インパルス応答関数, VARをuSEM で推定 • mlVAR – graphicalVARのマルチレベル版 VARモデル推定: Rの場合 32 Sacha Epskamp (2016). graphicalVAR: Graphical VAR for Experience Sampling Data. R package version 0.1.4. https://CRAN.R-project.org/package=graphicalVAR Xiao Yang, Nilam Ram and Peter Molenaar (2018). pompom: Person-Oriented Method and Perturbation on the Model. R package version 0.2.0. https://CRAN.R-project.org/package=pompom Sacha Epskamp, Marie K. Deserno and Laura F.Bringmann (2019). mlVAR: Multi-Level Vector Autoregression. R package version 0.4.2. https://CRAN.R project.org/package=mlVAR

VARモデルの拡張 33

VARのモデルの限界 34 • 時系列のどの時点でもパラメータは一定 → 時点によって変化するVARモデル time-varying VAR model • 等間隔時間が前提 → 連続時間 VAR continuous-time VAR • 誤差に正規性を仮定 → 一般化VARモデルなど？ Bringmann, L. F., Ferrer, E., Hamaker, E. L., Borsboom, D., & Tuerlinckx, F. (2018). Modeling nonstationary emotion dynamics in dyads using a time-varying vector-autoregressive model. Multivariate behavioral research, 53(3), 293-314. Driver, C. C., & Voelkle, M. C. (2017). Introduction to Hierarchical Continuous Time Dynamic Modelling With ctsem. R package Vignette. Available online at: https://cran. r-project. org/web/packages/ctsem/index. html.

時変VAR 35 • 時間経過で、VARのパラメータが変わる • 一般化加法モデル (GAM)で推定 – GLMで説明変数の最適な次数をデータドリブンに決定 男性の幸せが女性の幸せに90日間ずっと影響しない 60日以降、女性の幸せは男性の幸せを下げる 90日間の1組のカップルのポジティブ感情の時変VAR → このカップルはうまくいってない（確信） Yt = Ct +φt Yt-1 +εt

連続時間 (V)AR 36 • 確率微分方程式 α 切片ベクトル B ドリフト行列 （状態のクロス回帰と自己回帰のパラメータ） G ディフージョン行列 （状態の誤差相関と自己相関） W(t) 連続時間誤差 (ウェーナープロセス、ブラウニア運動) AR(1)の連続時間版はOrnstein-Uhlenbeck process (OUP) OUPの2変量版が連続時間VAR、階層モデルも可能 (ベイズで) Oravecz, Z., Tuerlinckx, F., & Vandekerckhove, J. (2011). A hierarchical latent stochastic differential equation model for affective dynamics. Psychological methods, 16(4), 468.

連続時間AR(1) 37 150名成人(18-89才)を対象に、加齢と感情調整のダイナミクスとの関 係をESMで検討 Attractor point (μ): 得点の均衡状態 平均的な得点 BPS-related reactivity (γ)：外的要因への反応性 得点の分散：情動システムへの外的な要因による影響 Attractor strength(β): attractorへの引き戻り速度 自己相関：情動システムへの内的な制御力 Wood, J., Oravecz, Z., Vogel, N., Benson, L., Chow, S. M., Cole, P., ... & Ram, N. (2017). Modeling intraindividual dynamics using stochastic differential equations: Age differences in affect regulation. The Journals of Gerontology: Series B, 73(1), 171- 184.

• 時変VAR – mgcv: 一般化加法モデルのパッケージ • 連続時間AR or VARモデル – ctsem （階層ベイズ） • 連続時間モデルをsemやmcmcで解く – Matlabのツールボックス（階層ベイズ） いずれもRで簡単に実施可能 38 Driver, C. C., Oud, J. H., & Voelkle, M. C. (2017). Continuous time structural equation modeling with R package ctsem. https://www.researchgate.net/profile/Joachim_Vandekerckhove/publication/265060395_BHOUM_A_ MATLAB_toolbox_for_Bayesian_Hierarchial_Ornstein- Uhlenbeck_modeling/links/54d11b6d0cf25ba0f040b4ee.pdf Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36

Take home message • 自己報告時系列データに、AR/VARモデルを適用 することで、自己回帰、クロス回帰が推定可能 • VARで、granger因果やインパルス応答を検討す ることで、時系列間の因果的な関係を検討可能 • VARはネットワークモデルを組み合わせると解釈 に深みが出る • 時変VARや連続時間モデルでさらに深みが出る 複雑なモデルも無料（R）で簡単に推定できる素敵な時 代に僕らは生きている 39

Take home message （ポエム） • VARモデルや連続時間モデルなどで時系列デー タを扱えるが、その解釈を支えるのは、測定指標 の背景に想定している心理学的な理論。 • 心理学的な構成概念のダイナミックな変化を数理 的にモデル化することで、直接的にモデルを評価 できるし、パラメータの解釈がより有意味なものと なる。 • 時系列解析を活用して、自己報告データであって も活発にモデリングを。それが心理学の理論的発 展に大きく貢献するだろう。 40