Pythonで構造方程式モデリング (SEM) 〜semopyのチュートリアルを試してみる〜

構造方程式モデリング (SEM) をPythonでやりたいと思った時に、使用するパッケージとして `semopy` が候補に上がるかと思いますが、公式チュートリアルが何をやっているのか少し分かりにくく感じたので、動かしつつ解説をつけてみました。

* 公式チュートリアル : [Tutorial - semopy](https://semopy.com/tutorial.html)

# 構造方程式モデリング (SEM) とは

**構造方程式モデリング** (SEM : Structural Equation Modeling) とは、変数間の関係性に関する仮説をモデル化し、その妥当性を検証するための統計的手法です。「共分散構造分析」などとも呼ばれます。

## ■ SEMの特徴

SEMの特徴として以下のようなものが挙げられます。

* 仮定した変数間の関係性を統計的に検証できる。
* **潜在変数** (構成概念) を導入した分析ができる。
* **パス図**を用いたビジュアル表現ができる。
* 回帰分析・因子分析など、下位モデルに様々な統計モデルを含む。
* 関係性の「探索」ではなく、「検証」に主眼を置いている。

## ■ 分析の流れ

実際の分析では、以下のようなPDCAサイクルを回すことになると思います。

![sem_pdca.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/d207d085-035f-b963-d395-2e573dd2143c.png)


## ■ SEMにおける重要な概念・用語

### 観測変数と潜在変数

* **「観測変数」** :
    - 直接観測される変数 
    - 例: テストの点数、身長 etc.
* **「潜在変数」** : 
    - 直接は観測されない共通の要因に関する変数
    - 因子分析でいうところの「因子」に相当。SEMでは「構成概念」と呼ばれる。
    - 例: 数学の能力、ブランド力 etc.

パス図では観測変数は四角形で、潜在変数は円形で表現されます。

![variates.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/42ef05c5-66b7-c5a2-b813-94c67d99e642.png)


### 構造方程式と測定方程式

因子分析では説明変数と因子の関係性しか記述しませんが、SEMでは因子（構成概念）同士の関係性も記述することができます。これらの関係性の記述は、事前の研究結果などから導かれた仮説に基づきます。

* **「構造方程式」** ... 仮説に基づいた因果関係を記述する方程式
* **「測定方程式」** ... 潜在変数とそれを構成する観測変数の関係性を記述する方程式

## ■ パス図による表現

![スクリーンショット 2020-12-04 13.21.36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/af5a1139-b66e-065b-4214-7db5ecbe207e.png)


# semopyで構造方程式モデリング

`pip install semopy` で `semopy` パッケージをインストールできますが、Python3.7以外では正常に動作しない可能性があります。筆者は pyenv (python3.7.6) + jupyter-lab 環境で実行しています。

```python:In
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from see import see

import semopy
from semopy import Model
```

## ■ 使用するデータセットについて

semopyがサンプルデータセットとして用意している "Political Democracy" を使用します。発展途上国における産業化と政治的民主主義に関するデータセットです。各変数の説明は以下の通りです。

|変数名|説明|
|:--|:--|
|**y1**|報道の自由の専門家評価 (1960)|
|**y2**|野党の自由 (1960)|
|**y3**|選挙の公平性 (1960)|
|**y4**|選出された議会の有効性 (1960)|
|**y5**|報道の自由の専門家評価 (1965)|
|**y6**|野党の自由 (1965)|
|**y7**|選挙の公平性 (1965)|
|**y8**|選出された議会の有効性 (1965)|
|**x1**|一人当たりGNP (1960)|
|**x2**|一人当たり|エネルギー消費 (1960)
|**x3**|産業における労働力の割合(%) (1960)|

表から以下のようなことが読みとれます。
* `y1`~`y4` ... 1960年の民主主義に関する観測変数 --> `dem60`
* `y5`~`y8` ... 1965年の民主主義に関する観測変数 --> `dem65`
* `x1`~`x3` ... 1960年の各種産業指数 --> `ind60`

```python:In
data = semopy.examples.political_democracy.get_data()
print(data.head())
```

```:Out
      y1        y2        y3        y4        y5        y6        y7        y8        x1        x2        x3
1   2.50  0.000000  3.333333  0.000000  1.250000  0.000000  3.726360  3.333333  4.442651  3.637586  2.557615
2   1.25  0.000000  3.333333  0.000000  6.250000  1.100000  6.666666  0.736999  5.384495  5.062595  3.568079
3   7.50  8.800000  9.999998  9.199991  8.750000  8.094061  9.999998  8.211809  5.961005  6.255750  5.224433
4   8.90  8.800000  9.999998  9.199991  8.907948  8.127979  9.999998  4.615086  6.285998  7.567863  6.267495
5  10.00  3.333333  9.999998  6.666666  7.500000  3.333333  9.999998  6.666666  5.863631  6.818924  4.573679
```

## 【Plan】 仮説を立て、モデルを記述する

semopyではモデルをstr型で記述します。サンプルデータセットに対応するモデル記述が既に用意されているので、今回はそちらを利用します。パラメータ推定の結果、モデルに対するデータの当てはまりが良くない場合には、このモデル記述を編集することになります。

```python:In
desc = semopy.examples.political_democracy.get_model()
print(desc)
```

```:Out
# measurement model
ind60 =~ x1 + x2 + x3
dem60 =~ y1 + y2 + y3 + y4
dem65 =~ y5 + y6 + y7 + y8
# regressions
dem60 ~ ind60
dem65 ~ ind60 + dem60
# residual correlations
y1 ~~ y5
y2 ~~ y4 + y6
y3 ~~ y7
y4 ~~ y8
y6 ~~ y8
```

モデルの記述方法はsemopyの定義する構文に基づいています。詳細は[公式ドキュメント](https://semopy.com/syntax.html)をご参照ください。

* `~` : **回帰演算子** --> 測定方程式を記述する
* `=~` : **測定演算子** --> 構造方程式を記述する
* `~~` : **分散演算子**

## 【Do】 モデルを定義し、パラメータを推定する

`semopy.Model` に引数としてモデル記述を指定し、モデルを定義します。`fit()` メソッドを呼び出してパラメータを推定します。

```python:In
model = Model(desc)
result = model.fit(data)
print(result)
```

```:Out
Name of objective: MLW
Optimization method: SLSQP
Optimization successful.
Optimization terminated successfully
Objective value: 0.508
Number of iterations: 52
Params: 2.180 1.819 1.257 1.058 1.265 1.186 1.280 1.266 1.482 0.572 0.838 0.624 1.893 1.320 2.156 7.385 0.793 5.067 0.347 3.148 1.357 4.954 0.467 0.172 0.082 0.120 3.951 2.352 3.256 3.430 0.448
```

`inspect()` メソッドで、パラメータ推定結果を確認できます。

```python:In
inspect = model.inspect()
print(inspect)
```

```:Out
     lval  op   rval  Estimate   Std. Err   z-value      p-value
0   dem60   ~  ind60  1.482379   0.399024   3.71502   0.00020319
1   dem65   ~  ind60  0.571912   0.221383   2.58336   0.00978421
2   dem65   ~  dem60  0.837574  0.0984456   8.50799            0
3      x1   ~  ind60  1.000000          -         -            -
4      x2   ~  ind60  2.180494   0.138565   15.7363            0
5      x3   ~  ind60  1.818546   0.151993   11.9646            0
6      y1   ~  dem60  1.000000          -         -            -
7      y2   ~  dem60  1.256819   0.182687   6.87965  6.00009e-12
8      y3   ~  dem60  1.058174   0.151521    6.9837  2.87503e-12
9      y4   ~  dem60  1.265186   0.145151   8.71634            0
10     y5   ~  dem65  1.000000          -         -            -
11     y6   ~  dem65  1.185743   0.168908   7.02003  2.21823e-12
12     y7   ~  dem65  1.279717   0.159996   7.99841  1.33227e-15
13     y8   ~  dem65  1.266084   0.158238   8.00114  1.33227e-15
14  dem65  ~~  dem65  0.172210   0.214861  0.801494     0.422846
15  dem60  ~~  dem60  3.950849   0.920451    4.2923  1.76835e-05
16  ind60  ~~  ind60  0.448321  0.0866766   5.17234  2.31175e-07
17     y1  ~~     y5  0.624423   0.358435   1.74208    0.0814939
18     y1  ~~     y1  1.892743    0.44456   4.25756  2.06666e-05
19     y2  ~~     y4  1.319589    0.70268   1.87794    0.0603898
20     y2  ~~     y6  2.156164   0.734155   2.93693   0.00331475
21     y2  ~~     y2  7.385292    1.37567    5.3685  7.93938e-08
22     y3  ~~     y7  0.793329   0.607642   1.30558     0.191694
23     y3  ~~     y3  5.066628   0.951722   5.32365  1.01708e-07
24     y4  ~~     y8  0.347222   0.442234  0.785154     0.432363
25     y4  ~~     y4  3.147911   0.738841    4.2606  2.03874e-05
26     y6  ~~     y8  1.357037     0.5685   2.38705    0.0169843
27     y6  ~~     y6  4.954364   0.914284   5.41884   5.9986e-08
28     x3  ~~     x3  0.466732  0.0901676   5.17628  2.26359e-07
29     x1  ~~     x1  0.081573  0.0194949   4.18432  2.86025e-05
30     x2  ~~     x2  0.119894  0.0697474   1.71897    0.0856192
31     y5  ~~     y5  2.351910   0.480369   4.89604  9.77851e-07
32     y8  ~~     y8  3.256389    0.69504   4.68518  2.79711e-06
33     y7  ~~     y7  3.430032   0.712732   4.81251  1.49045e-06
```

適合後のモデルに対して、パス図を描画・保存することができます。(`graphviz`パッケージが必要になります)

```python:In
semopy.semplot(model, "graph.png")
```

![graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/5460ace2-ce87-4ca0-e475-40b9ac5aa9cb.png)


## 【Check & Action】 データの当てはまりの良さを確認する

`calc_stats` メソッドにより、各種適合度指標を確認できます。指標をもとに当てはまりが良くないと判断した場合、モデルの修正を検討します。

```python:In
stats = semopy.calc_stats(model)
print(stats.T)
```

```:Out
                    Value
DoF             35.000000
DoF Baseline    55.000000
chi2            38.125446
chi2 p-value     0.329171
chi2 Baseline  730.654577
CFI              0.995374
GFI              0.947820
AGFI             0.918003
NFI              0.947820
TLI              0.992731
RMSEA            0.034738
AIC             60.983321
BIC            132.825453
LogLik           0.508339
```

以上です。

# 参考

* [semopy: A Python package for Structural Equation Modeling](https://arxiv.org/pdf/1905.09376.pdf)
* [PoliticalDemocracy: Industrialization And Political Democracy Dataset - rdrr.io](https://rdrr.io/cran/lavaan/man/PoliticalDemocracy.html)
* [SEM(構造方程式モデリング)とは - (株)日科技研](https://www.i-juse.co.jp/statistics/product/func/sem/sem.html)
* [構造方程式モデリングの話（初級編） - slideshare](https://www.slideshare.net/DaikiNakamura/ss-152713829)
* [100614 構造方程式モデリング基本の「き」 - slideshare](https://www.slideshare.net/MasahiroShinohara/100614-4502873)
* [構造方程式モデリングの基礎 - 京都大学大学院 教育学研究科](http://cogpsy.educ.kyoto-u.ac.jp/personal/Kusumi/datasem09/090708.pdf)
* [構造方程式モデリング ―モデル構築䛾再検討―](https://www.mizumot.com/method/2012-05_Konno.pdf)