TWS-007「心理ネットワークアプローチ入門」

class: center, middle, inverse, title-slide

.title[
# <font size='5'>TWS-007「心理ネットワークアプローチ入門」</font>
]
.subtitle[
## 横断的ネットワーク解析：基礎編
]
.author[
### 国里愛彦（専修大学）<img style="float:right" src="fig/yk.png" width="250">
]
.date[
### 日本心理学会(2022年)
]

---

# 配布資料

![](fig/qr.png)
---
## 心理ネットワークとは？

- 項目（心理変数）間のネットワーク構造を検討する方法.med[(Epskamp, Borsboom & Fried, 2018)]

- ネットワークは，観察可能な変数を表すノード（円）と統計的関係を表すエッジ（線）から構成される

![](fig/net01.png)

---

## 心理ネットワーク分析の手順

1. 統計的手法によってノード間のエッジを推定する

2. 推定されたエッジをもとにネットワークを図示する

3. ネットワーク特性の指標を検討する

4. ネットワークの正確度・安定性を検討する

</br>

- エッジはデータから推定される（重み付きネットワーク）

→(1)サンプルサイズに影響されるので正確度・安定性を検討

→(2)グローバルなネットワーク特性.med[(small worldnessなど)]は使えないので，ローカルなネットワーク特性.med[(中心性指標)]を用いる
---

## 有向・無向ネットワーク

- 有向ネットワークはエッジに方向性があるネットワーク.med[(左図）]，無向ネットワークは方向性がないネットワーク.med[（右図）]

- エッジの矢印は因果関係を示唆しているが，横断データではそれは難しい(→無向ネットワークの利用)

![](fig/net06.png)
---

## エッジの見方

- **エッジの色** 符号を意味しており，青・緑の場合は正，赤の場合は負を表す

- **エッジの太さ** 全体のノードからの影響を考慮した時の２つのノード間の関係の強さ

- **エッジの長さ** エッジの強さの逆数。関係が強ければ強いほど短くなる

![](fig/net01.png)
---
## ペアワイズマルコフ確率場

- 無向ネットワークモデルの推定では，ペアワイズ・マルコフ確率場.med[(pairwise Markov random field)]がよく使われる.med[(Epskamp, Haslbeck, Isvoranu,& Van Borkulo, 2022)]

- 心理ネットワーク分析の主なペアワイズ・マルコフ確率場のモデルには以下の３つがある

1. **Gaussian graphical model (GGM)**: 連続変数のデータ.med[(エッジの値は-1から1の範囲)]

2. **Ising model**: 2値変数のデータ.med[(エッジの値は-∞から∞の範囲)]

3. **mixed graphical models (MGMs)**: 2値・順序・連続変数の混ざったデータ.med[(エッジの値は-1から1の範囲)]

---
## Gaussian graphical model(GGM)

- GGMでは，多変量正規分布を用いて，連続変数間のネットワークを推定する
- `$\boldsymbol{Y}$`.med[(ランダムに選んだ参加者からのn個の反応のベクトル)]は，平均 `$\boldsymbol{\mu}$` と分散共分散行列 `$\boldsymbol{\Sigma}$` の多変量正規分布に従う

$$
\boldsymbol{Y} \sim N\left(\boldsymbol{\mu}, \boldsymbol{\Sigma}\right)
$$

- GGMでは，分散共分散行列 `$\boldsymbol{\Sigma}$`からその逆行列である精度行列 `$\boldsymbol{K}$` を計算し( `$\boldsymbol{K}=\boldsymbol{\Sigma}^{-1}$`)，精度行列 `$\boldsymbol{K}$` から偏相関行列を計算する
- 偏相関係数をエッジにしてネットワークをプロットする

※具体的な計算例は[配布ドキュメント](https://ykunisato.github.io/ccp-lab-slide/tws-psynet-2022/doc-psynet-basic.html)参照

---
## Ising モデル

- 元々は磁性体に関する統計力学的モデルだが，２値データのネットワークの記述に使える

- ２値反応のペアの同時確率をモデル化.med[(Epskamp et al., 2022) ]

`$$\operatorname{Pr}(\boldsymbol{Y}=\boldsymbol{y})=\frac{1}{Z} \exp \left(\sum_{i} y_{i} \tau_{i}+\sum_{<i, j>} y_{i} y_{j} \omega_{i j}\right)%$$`  
- `$\tau_{i}$` は閾値.med[(どっちの反応が出やすいか）]， `$\omega_{i j}$` はノード間の反応の類似度を表す。Zは確率の合計が１になるように機能する

- Zの計算は， `$\boldsymbol{Y}$` の全ての可能性について考慮するので20変数以上の多変量推定は難しい.med[(その場合は，単変量推定する)]

---
## Mixed graphical models(MGMs)

- 2値・順序・連続変数が混ざったデータに適用できる.med[(Epskamp et al., 2022)]

- データのタイプに合わせてノードごとに確率分布が適用.med[(指数分布族: 正規分布，ポワソン分布，指数分布，多項分布）]

- ノードごとに一般化線形モデルによって推定し，推定されたパラメータが最終モデルで統合される

- MGMsは，ペアワイズである必要もなくて，高次の交互作用を含めたmoderated network modelも検討できる

---
## モデル選択

- ネットワークは，エッジが１つもないものから全てのノード間にエッジがある飽和モデルの間に，可能性のあるネットワークが大量にある.med[(Blanken, Isvoranu & Epskamp, 2022)]

- オッカムの剃刀に従い，パフォーマンスを下げないようにして，シンプルなモデルを選ぶ

- 心理ネットワーク分析のモデル選択には以下の４つがある

1. しきい値(Thresholding)

2. 刈り込み(Pruning)

3. モデル探索(model search)

4. 正則化(regularization)

---
**しきい値(Thresholding)**

- しきい値では，ある基準.med[（p値, false discovery rate, credibility interval, Bayes factor）]に基づいて，エッジを除去する.med[（エッジを0に設定）]

- 基準以下のエッジを見せてないだけで影響は残る

**刈り込み(Pruning)**

- 刈り込みでは，ある基準.med[（p値, false discovery rate, credibility interval, Bayes factor）]に基づいて，除去した上で.med[（エッジを0に設定）]，**再推定**する

- 再推定するので，基準以下のエッジを0としたときのネットワークが推定される

---
**モデル探索(model search)**

- 可能なネットワークを反復的に探索して最もデータフィットしたモデルを選ぶ.med[（AIC, BIC, EBICなどの情報量規準を用いる])

**正則化(regularization)**

- モデルの複雑さに罰則をかけて，ネットワークを疎.med[(sparse)]にする。影響力が弱いエッジはゼロに近づき単純化できる

- LASSO .med[(least absokute shrinkage and selection operator)]の一種のGraphical LASSOを使う。罰則の強さの調整は，cross-validation(CV)やEBICを使う
---
### GGMでのモデル選択の推奨.med[(Isvoranu & Epskamp, 2021)]

(1)サンプルサイズが300の場合，ネットワーク構造の発見には，正則化が好ましい

(2)サンプルサイズが1000の場合，ネットワーク構造や強いエッジに興味があるなら正則化，特定のエッジに関心があるならモデル探索が好ましい

(3)サンプルサイズが5000の場合，真のネットワークの探索には，モデル探索が好ましい

### IsingモデルとMGMsでのモデル選択の推奨.med[(Blanken et al, 2022)]

Isingモデルの場合，小中サンプルサイズでは正則化手法が見やすく，大サンプルサイズでは非正則化手法も使える。現状，MGMsで使えるのは正則化手法のみ

---
** 順序・非正規変数での注意点.med[(Blanken, Isvoranu & Epskamp, 2022)]**

- 順序変数ではポリコリック相関を使うが，小サンプルサイズでは不安定になる.med[（psychonetricsパッケージの推定法やBGGMパッケージを使う）]
- 正規分布していないデータの場合は，スピアマンの順位相関係数を用いてGGMする

**サンプルサイズ**

- ノードは３０より少ないのが望ましく，できるだけサンプルサイズが大きいほうが望ましい.med[(Isvoranu & Epskamp, 2021)]
- ネットワークの構造を調べるより，特定のエッジに焦点をしぼった方が大きなサンプルサイズが必要

**欠測処理**

- bootnetでは欠測値除去が必要だが，psychonetricsでは完全情報最尤推定法が使える.med[(estimator = "FIML")]

---
### エッジの解釈法

- PRMFのエッジは，２変数以外の影響も考慮した２変数の関係を表す.med[（ネットワーク全体の影響を考慮する）]

- エッジの解釈としては以下の２つがある

(1)**予測可能性の示唆**: 他の変数の影響を考慮しても，ある変数から別の変数が予測できる可能性がある

(2)**因果関係の可能性の示唆**: 因果の方向性は不明だが，そのエッジにはなんらかの因果的な効果がある可能性はある.med[(逆にいうと，強い理論的予測による仮説検証の結果のような解釈はしない)]

---
## ノードの中心性指標

心理ネットワーク分析で使えるネットワークの指標として，以下の中心性指標がある

- **Strength**:あるノードがつながっている全てのエッジの強さの合計。あるノードの全体に対する影響力を示している

- **Closeness**: あるノードと他の全てのノード間の最短経路長の合計の逆数。あるノードがどのくらい間接的な影響を含めて他のノードとつながっているのかを示している

- **Betweeness**: 2つのノード間の最短経路上に，あるノードが何回あるのかを示している。あるノードが2つのノードの経路上でどのくらい重要かを示している

---
## エッジの正確度と中心性指標の安定性

- 推定されたエッジはサンプルやサンプルサイズの影響をうけるので，正確度(accuracy)や安定性(Stability)を検討する.med[(Epskamp, Borsboom,& Fried, 2018)]

1. ブートストラップ法による信頼区間をプロットすることでエッジの重みの正確度を検討する .med[(有意性の検討には使わない)]

2. データのサブセットを用いた場合の中心性指標の安定性を検討する .med[(サンプルサイズを小さくしても安定するか検討，CS係数が0.25未満は不適切であり0.5を超えている必要がある)]

3. エッジの重みと中心性指標間でブートストラップ差異検定を行って，それらが有意に異なるのかを検討する

.med[※1は必須，2は中心性指標を使う場合，3は目的に合わせて]
---
## Rパッケージ

- アムステルダム大学の研究者達が開発した以下のRパッケージで，横断データの心理ネットワーク分析は実施できる

- [qgraph](https://cran.r-project.org/web/packages/qgraph/index.html):ネットワーク推定(GGM)とプロット

- [bootnet](https://cran.r-project.org/web/packages/bootnet/index.html): ネットワーク推定(GGM & Isingモデル)，ネットワークの正確度と安定性の検討

- [psychonetrics](https://cran.r-project.org/web/packages/psychonetrics/index.html)：より洗練化した書き方でネットワーク推定(GGM & Isingモデル)

- [mgm](https://cran.r-project.org/web/packages/mgm/index.html):ネットワーク推定(MGMs)

アムステルダム大学の統計ソフト[JASP](https://jasp-stats.org/)には，心理ネットワーク分析が入っており，ちょっと試すのに便利

---
### 心理ネットワーク分析を試してみよう！
.med[
[Jordan et al.(2017)](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0182162)が，3404名からGAD-7(不安)，PHQ-9(うつ),PHQ-15(身体症状)を収集したデータで試す
]
### 使用するRパッケージの読み込み

```r
library(tidyverse)
library(foreign)
library(bootnet)
library(qgraph)
library(psychonetrics)
```

---
### データの読み込み
.med[
Jordan et al.(2017)のデータをダウンロードして，読み込む
]

```r
download.file("https://doi.org/10.1371/journal.pone.0182162.s004","pone.0182162.s004.sav")
data <- read.spss("pone.0182162.s004.sav", to.data.frame=TRUE)
```

.med[
GAD-7（不安に関する７項目）のデータを使う。renameで変数名を整理する
]

```r
# データの整理
data_gad <- data %>% 
  rename(gad7a = S_GAD7_a, gad7b = S_GAD7_b, gad7c = S_GAD7_c, gad7d = S_GAD7_d,gad7e = S_GAD7_e, gad7f = S_GAD7_f, gad7g = S_GAD7_g) %>% 
  select(gad7a, gad7b, gad7c, gad7d, gad7e, gad7f, gad7g)
```

---
### モデル選択： しきい値
.med[
bootnetのestimateNetwork()関数で，thresholdとalphaを指定する

GAD-7は4件法で順序変数なので，「corMethod = "cor_auto"」と指定して，自動的にポリコリック相関で推定
]

```r
results_gad <- estimateNetwork(data_gad,default = "pcor", threshold = "sig", alpha = 0.05, corMethod = "cor_auto")
plot(results_gad, theme = "colorblind", cut = 0,layout = "spring",labels = TRUE)
```

![](slide-psynet-basic_files/figure-html/unnamed-chunk-5-1.png)

---
### モデル選択： 刈り込み
.med[
psychonetricsパッケージのprune()関数を用いる。psychonetricsはパイプ演算子を使って可読性の高いコードが書きやすい
]

```r
ggm(data_gad) %>%
  runmodel %>%
  prune(alpha = 0.05) %>%
  getmatrix("omega") %>%
  qgraph(theme = "colorblind", cut = 0,layout = "spring",labels = TRUE)
```

![](slide-psynet-basic_files/figure-html/unnamed-chunk-6-1.png)

---
### モデル選択：モデル探索

.med[モデル探索はmodelsearch()で行う]

```r
ggm(data_gad) %>% runmodel %>%
  modelsearch(criterion = "bic") %>%
  getmatrix("omega") %>% qgraph(theme = "colorblind", cut = 0,layout = "spring",labels = TRUE)
```

![](slide-psynet-basic_files/figure-html/unnamed-chunk-7-1.png)

---
### モデル選択：正則化
.med[bootnetのestimateNetwork()関数を使って，EBICでパラメータ調整をするGLASSOを使う(「default="EBICglasso"」を指定する）]

```r
results_gad <- estimateNetwork(data_gad,default = "EBICglasso", corMethod = "cor_auto")
plot(results_gad, theme = "colorblind", cut = 0,layout = "spring",labels = TRUE)
```

![](slide-psynet-basic_files/figure-html/unnamed-chunk-8-1.png)

---
### エッジの重みの正確度
.med[bootnetのbootnet()関数を使う。summaryでも推定結果を確認できるが，プロットのほうが確認がしやすい]

```r
accuracy_edge <- bootnet(results_gad, nBoots = 2500, nCores =4, statistics =   c("edge", "strength", "closeness", "betweenness"))
plot(accuracy_edge, labels = FALSE, order = "sample")
```

---
## 中心性の指標のプロット

.med[qgraphパッケージのcontralityPlot関数でプロットできる。デフォルトではStrengthのみが出力されるので，includeで"Strength", "Betweenness", "Closeness"の３つを指定する]

```r
centralityPlot(results_gad, include = c("Strength", "Betweenness", "Closeness"))
```

![](slide-psynet-basic_files/figure-html/unnamed-chunk-10-1.png)

---
## 中心性指標の安定性(Stability)
.med[
bootnet()関数を使う。「type = ”case”」を指定するとケースドロップ時の安定性を推定できる。statisticsで"strength", "closeness", "betweenness"を指定する
]

```r
stability_centrality <- bootnet(results_gad, nBoots = 2500, type = "case", nCores =4, statistics =  c("strength", "closeness", "betweenness"))
```

---
.med[
CS係数を確認する。StrengthとClosenessはCS係数が0.5を超えているが，Betweenessは低いので，解釈を控える必要がある
]

```r
corStability(stability_centrality)
```

![](fig/plot3.png)
---
## 有意差検定
.med[
bootnetのdifferenceTest()関数を使う

Strengthがgad7bとgad7fで差があるかを検討すると，95%CIが0をまたいでおらず，gad7bはgad7fよりも有意にstrengthが大きい
]

```r
differenceTest(accuracy_edge, 2, 6, "strength")  
```

![](fig/plot4.png)
---
エッジ間の差の検定結果をプロットする.med[(エッジ間で有意な差がある場合は黒く塗りつぶされ，有意ではない場合は灰色に塗りつぶされる)]

```r
plot(accuracy_edge, "edge", plot = "difference", onlyNonZero = TRUE, order = "sample")  
```

![](fig/plot5.png)
---
Strengthにおけるノード間の差の検定結果をプロットする.med[(ノード間で有意な差がある場合は黒色，有意ではない場合は灰色に塗りつぶす)]

```r
plot(accuracy_edge,"strength")
```

![](fig/plot6.png)
---
## まとめ

- 心理ネットワーク解析は新しい方法だが，偏相関をエッジとしたネットワークを描く方法と考えると難しくない

- 便利なRパッケージも増えており実施は簡単だが，モデル選択・エッジ・中心性指標の評価などは理論的な理解も必要

- 日々新しい方法の開発と解析法の推奨の変化があるので，アムステルダム大学の[Psych Systems](http://psychosystems.org/)のサイトを確認したり，最新の動向をおさえる

</br>

※ [配布ドキュメント](https://ykunisato.github.io/ccp-lab-slide/tws-psynet-2022/doc-psynet-basic.html)に，IsingモデルやMGMsを用いた解析法や詳細な補足説明を書いているので参照ください

---
## 参考文献
.med[
- Blanken, T. F., Isvoranu, A. M., & Epskamp, S. (2022). Chapter 7. Estimating network structures using model selection. In Isvoranu, A. M., Epskamp, S., Waldorp, L. J., & Borsboom, D. (Eds.). Network psychometrics with R: A guide for behavioral and social scientists. Routledge, Taylor & Francis Group.
- Epskamp, S., Borsboom, D., & Fried, E. I. (2018). Estimating psychological networks and their accuracy: A tutorial paper. Behavior Research Methods, 50(1), 195–212.
- Epskamp, S., Haslbeck, J. M. B., Isvoranu, A. M., & Van Borkulo, C. D. (2022). Chapter 6. Pairwise Markov random fields. In Isvoranu, A. M., Epskamp, S., Waldorp, L. J., & Borsboom, D. (Eds.). Network psychometrics with R: A guide for behavioral and social scientists. Routledge, Taylor & Francis Group.
- Isvoranu, A., & Epskamp, S. (2021). Which Estimation Method to Choose in Network Psychometrics? Deriving Guidelines for Applied Researchers. preprint
- Jordan, P., Shedden-Mora, M. C., & Löwe, B. (2017). Psychometric analysis of the Generalized Anxiety Disorder scale (GAD-7) in primary care using modern item response theory. PloS One, 12(8), e0182162.
]