ネットワークにおける変化のモデリング

ネットワークは変化する？

• ベクトル自己回帰(VAR)モデルは基本的にパラメータの定常性（時間によって変化しないこと）を仮定している。

• 心理変数間の関係は時間に対して一定？心理的介入の前後では関係が変わるのでは？

時変ネットワークモデル

時間的変化は連続？離散？

• ネットワークのパラメータの変化には連続的変化(左)と離散的変化(右，状態間でジャンプ有り)がある

• 連続か離散かによって解析が変わるので，理論的・実証的（先行研究）にどちらが適切か検討する

時変ネットワークの推定

連続的な変化を仮定する場合

• 移動窓アプローチ（後述），一般化加法モデリング(パラメータを時間の関数として明示的にモデル化)，融合LASSOなどがある。

離散的な変化を仮定する場合

• パラメータの変化に関係する変数（治療の前後を表す変数）が調整変数か検討，隠れマルコフモデル（データから隠れ状態の遷移を推定するモデル），レジーム・スイッチングモデル（ある閾値を超えるかで状態が変化するモデル）

連続的な変化（移動窓アプローチ）

• 時系列上で窓を移動をさせていき，その窓ごとに局所モデルを多数推定する方法

• 固定窓と異なり，全てのデータポイントを利用するが特定の時点に重みをつけて推定（カーネル重み付け法，詳細は Haslbeck et al. (2020) ）

帯域幅の検討

• 移動窓アプローチでは，カーネル関数の幅（帯域幅）が重要になる

• 帯域幅が広いと多くのデータポイントを使うので，推定は安定するが時間変化の感度が下がる（オレンジ）。狭いと逆(緑)。

使用データ

• Journal of Open Psychology Dataに掲載された Kossakowski et al. (2017) のデータを用いる。

• Kossakowski et al. (2017) は，１名のうつ病患者に239日間さまざまな症状を収集（1日10回で合計1478回答）

• データはOSF(https://osf.io/j4fg8/)で配布されているので，ESMdata.zipをダウンロードして展開する。

• 同じデータがmgm(Haslbeck & Waldorp, 2020)パッケージに入っているが(symptom_data)，実際の解析に近づけるためにcsvファイルから読み込む。

使用パッケージ

• データの前処理にtidyverse(Wickham et al., 2019)を用い，日付や時間の情報の処理にlubridate(Grolemund & Wickham, 2011)を用いる

• 時変自己相関モデルにmgm(Haslbeck & Waldorp, 2020)を用い，可視化にqgraph(Epskamp et al., 2012)を用いる

• mgmの使用法は，Haslbeckのブログ(Haslbeck, 2020)が参考になる

library("tidyverse")
library("lubridate") 
library("mgm")
library("qgraph")

データの読み込みと日付データの処理

• csvファイルを読み込んで，dmy_hms()で日付データを処理する
• time_normは最初の測定が0，最後の測定が1になるように標準化

data <- read_csv("ESMdata/ESMdata.csv") |> 
  mutate(datetime = dmy_hms(paste(date, resptime_e))) 

data$time_norm <- as.numeric(data$datetime-data$datetime[1])/as.numeric((data$datetime-data$datetime[1])[nrow(data)])
head(data)

# A tibble: 6 × 88
   ...1 date     phase concentrat dayno beepno beeptime resptime_s resptime_e
  <dbl> <chr>    <dbl>      <dbl> <dbl>  <dbl> <time>   <time>     <time>    
1     1 13/08/12     1        150   226      1 08:58    08:58:56   09:00:15  
2     2 14/08/12     1        150   227      5 14:32    14:32:09   14:33:25  
3     3 14/08/12     1        150   227      6 16:17    16:17:13   16:23:16  
4     4 14/08/12     1        150   227      8 18:04    18:04:10   18:06:29  
5     5 14/08/12     1        150   227      9 20:57    20:58:23   21:00:18  
6     6 14/08/12     1        150   227     10 21:54    21:54:15   21:56:05  
# ℹ 79 more variables: resp_abort <dbl>, mood_relaxed <dbl>, mood_down <dbl>,
#   mood_irritat <dbl>, mood_satisfi <dbl>, mood_lonely <dbl>,
#   mood_anxious <dbl>, mood_enthus <dbl>, mood_suspic <dbl>,
#   mood_cheerf <dbl>, mood_guilty <dbl>, mood_doubt <dbl>, mood_strong <dbl>,
#   pat_restl <dbl>, pat_agitate <dbl>, pat_worry <dbl>, pat_concent <dbl>,
#   se_selflike <dbl>, se_ashamed <dbl>, se_selfdoub <dbl>, se_handle <dbl>,
#   soc_who1 <dbl>, soc_enjoy_alone <dbl>, soc_prefcomp <dbl>, …

変数を選択し，欠測を除外する

• select()で使う変数を選択（気分12変数，時間関連5変数）

• 欠測があるので（５つの測定），filter()とcomplete.cases()で完全ケースに絞る

data_select <- data |> 
  select(mood_relaxed,mood_down,mood_irritat,
         mood_satisfi,mood_lonely,mood_anxious,
         mood_enthus,mood_suspic,mood_cheerf,
         mood_guilty,mood_doubt,mood_strong,
         datetime,time_norm,beepno,dayno,phase)

data_complete <- data_select |> 
  filter(complete.cases(data_select))

データの可視化(対象者の状態)

• 対象者の状態（1＝ベースライン，２=減薬前:盲検，3=減薬中:盲検，4=減薬後，5=研究後）をggplot2でプロットする(aes()でx軸やy軸を決め，geom_line()で折れ線グラフを描く)。

data_complete |> 
  ggplot(aes(x=datetime,y=phase)) + 
  geom_line() 

データの可視化(1)

• データを前処理し（selectで-を使って使わない変数を除外，pivot_longerでロング型に，mutateで変数を因子型に），ggplot2でプロット。情報多すぎ…

data_complete |> 
  select(-time_norm,-beepno,-dayno,-phase) |>   
  pivot_longer(cols=-datetime,
               names_to = "variables",
               values_to = "score") |> 
  mutate(variables = as.factor(variables)) |> 
  ggplot(aes(x=datetime,y=score, color=variables)) + 
  geom_line() +
  geom_point()

データの可視化(2)

• 全データだと見にくかったので，filterで最初の５日間にしぼってプロットしてみる。

data_complete |> 
  filter(datetime < "2012-08-18") |> 
  select(-time_norm,-beepno,-dayno,-phase) |>   
  pivot_longer(cols=-datetime,
               names_to = "variables",
               values_to = "score") |> 
  mutate(variables = as.factor(variables)) |> 
  ggplot(aes(x=datetime,y=score, color=variables)) + 
  geom_line() +
  geom_point()

tvmvar関数

• mgmパッケージのtvmvar()関数を使うと，カーネル重み付き法による時変混合ベクトル自己回帰モデルを用いることができる。
• mgmパッケージは，Mixed graphical models(MGMs) のパッケージであり，2値・順序・連続変数が混ざったデータに適用できる。
• mgmパッケージの詳細は，CRANのページを参照

帯域幅の決定

• カーネル重み付き法では帯域幅の設定が重要（広いと推定は安定するが時間変化の感度が下がる）

• bwSelect()でクロスバリデーションによって最も予測のエラーが小さい帯域幅を探索できる

帯域幅の決定(bwSelect)

• dataには回答データのみを入れる，typeは変数のタイプ（ガウス,ポワソン,カテゴリカル），levelは変数のカテゴリ数（連続は1にする），bwSeqに検討する帯域幅のベクトルをいれる（今回は0.01から1まで10等分），bwFoldsとbwFoldsizeはクロスバリデーションの設定， modeltypeはtvmvarの場合は”mvar”とする，lagsは１とした，timepointsに標準化したタイムポイント，beepvarやdayvarは時間の情報

set.seed(123)
bwSeq <- seq(0.01, 1, length = 10)
bw_object <- bwSelect(data = data_complete[,1:12],
                      type = rep("g", 12),
                      level = rep(1, 12),
                      bwSeq = bwSeq,
                      bwFolds = 1,
                      bwFoldsize = 10,
                      modeltype = "mvar",
                      lags = 1,
                      timepoints = data_complete$time_norm,
                      beepvar = data_complete$beepno,
                      dayvar = data_complete$dayno,
                      pbar = FALSE)
bandwidth <- bwSeq[which.min(bw_object$meanError)]
bandwidth

[1] 0.45

tvmvarによる推定

• 帯域幅が決まったので(0.45)，tvmvar()で時変混合ベクトル自己回帰モデルで推定を行う

• bwSelect()と同じ引数は省略。lambdaSelは正則化（今回はクロスバリデーション”CV”），estpointsは時間の分割（今回は20分割した），bandwidthにbwSelect()で決めた帯域幅を指定する。

result_tvmvar <- tvmvar(data = data_complete[,1:12],
                    type = rep("g", 12),
                    level = rep(1, 12), 
                    lambdaSel = "CV",
                    timepoints = data_complete$time_norm, 
                    estpoints = seq(0, 1, length = 20), 
                    bandwidth = bandwidth,
                    lags = 1,
                    beepvar = data_complete$beepno,
                    dayvar = data_complete$dayno,
                    pbar = TRUE)

予測可能性を検討

• predict()関数で予測のエラーを計算する。
• tvmvar()で推定した結果をobjectに指定し，errorConでRMSEとR二乗を選ぶ，tvMethod は”weighted”を選ぶ。

pred_tvmvar <- predict(object = result_tvmvar, 
                    data = data_complete[,1:12],
                    errorCon = c("R2", "RMSE"),
                    tvMethod = "weighted", 
                    beepvar = data_complete$beepno,
                    dayvar = data_complete$dayno)

推定されたモデルを可視化

• ノードのラベルを用意して，1・10・20のタイムポイントの結果を可視化する（推定したオブジェクト$wadjに重みの情報，$edgecolorに色の情報がある）
• ネットワークのプロットはqgraphを用いる

node_labels <- c("Relaxed","Down","Irritat","Satisfi","Lonely","Anxious","Enthus","Suspic","Cheer","Guilty","Doubt","Strong") 

par(mfrow=c(1,3))
for(tp in c(1,10,20)) {
  qgraph(t(result_tvmvar$wadj[, , 1, tp]), 
         layout = "circle",
         edge.color = t(result_tvmvar$edgecolor[, , 1, tp]), 
         labels = node_labels, 
         mar = rep(5, 4), 
         vsize=14, esize=15, asize=13,maximum = .5, 
         pie = pred_tvmvar$tverrors[[tp]][, 3],
         title = paste0("Estimation point = ", tp), title.cex=1.2)
}

推定されたモデルを可視化

パラメータの時間変動を可視化(1)

• あまりパラメータの時間変化のないもの(Down-Lonely:落ち込み-孤独)，時間に応じて増加したもの（Strong-Satisfi：力強い-満足），時間に応じて減少したもの(Doubt-Satisfi:迷う-満足)をプロット

tibble(time = 1:20,
       Doubt_Satisfi = result_tvmvar$wadj[4,11,1,],
       Down_Lonely = result_tvmvar$wadj[5,2,1,],
       Strong_Satisfi = result_tvmvar$wadj[4,12,1,]) |>
  pivot_longer(cols=-time,
               names_to = "variables",
               values_to = "score") |> 
  mutate(variables = as.factor(variables)) |> 
  ggplot(aes(x=time,y=score, color=variables)) + 
  geom_line() +
  geom_point()+
  labs(x = "時間", y = "偏相関係数",color="変数") +
  theme(text = element_text(family = "HiraKakuProN-W3")) 

パラメータの時間変動を可視化(2)

• phaseを同時にプロットしてみる
• 力強い-満足（Strong-Satisfi)は，減薬後から研究後に増加している？

phase <- NULL
split_20 <- seq(0, 1, length = 20)
for (i in 1:length(split_20)) {
  min_index <- which.min(abs(data_complete$time_norm - split_20[i]))
  phase[i] <- data_complete$phase[min_index]
}
phase <- (phase/5)*0.16
tibble(time = 1:20,
       Doubt_Satisfi = result_tvmvar$wadj[4,11,1,],
       Down_Lonely = result_tvmvar$wadj[5,2,1,],
       Strong_Satisfi = result_tvmvar$wadj[4,12,1,],
       Phase = phase) |>
  pivot_longer(cols=-time,
               names_to = "variables",
               values_to = "score") |> 
  mutate(variables = as.factor(variables)) |> 
  ggplot(aes(x=time,y=score, color=variables)) + 
  geom_line() +
  geom_point()+
  labs(x = "時間", y = "偏相関係数",color="変数") +
  theme(text = element_text(family = "HiraKakuProN-W3")) 

その他の推定方法

• ベクトル自己回帰ではなく，Mixed Graphical Models (MGMs)を用いてパラメータの時間変動を検討する場合は，mgmパッケージのtvmgm()が使える。
• 一般化加法モデル(パラメータを時間の関数として明示的にモデル化)を用いる場合は，tvvarGAMが使える。こちらは，CRANではなくgithub上で公開されている。

まとめ

• 時系列を通じてネットワークは定常？（時変ネットワークモデルを検討）
• 時間的変化は離散的か連続的か？
• 連続的変化は移動窓アプローチなど，離散的変化は隠れマルコフモデルなど
• 時変自己相関モデルには帯域幅の推定もできるmgmパッケージが便利