構造方程式モデリング (SEM) をPythonでやりたいと思った時に、使用するパッケージとして semopy が候補に上がるかと思いますが、公式チュートリアルが何をやっているのか少し分かりにくく感じたので、動かしつつ解説をつけてみました。
- 公式チュートリアル : Tutorial - semopy
構造方程式モデリング (SEM) とは
構造方程式モデリング (SEM : Structural Equation Modeling) とは、変数間の関係性に関する仮説をモデル化し、その妥当性を検証するための統計的手法です。「共分散構造分析」などとも呼ばれます。
■ SEMの特徴
SEMの特徴として以下のようなものが挙げられます。
- 仮定した変数間の関係性を統計的に検証できる。
- 潜在変数 (構成概念) を導入した分析ができる。
- パス図を用いたビジュアル表現ができる。
- 回帰分析・因子分析など、下位モデルに様々な統計モデルを含む。
- 関係性の「探索」ではなく、「検証」に主眼を置いている。
■ 分析の流れ
実際の分析では、以下のようなPDCAサイクルを回すことになると思います。
■ SEMにおける重要な概念・用語
観測変数と潜在変数
-
「観測変数」 :
- 直接観測される変数
- 例: テストの点数、身長 etc.
-
「潜在変数」 :
- 直接は観測されない共通の要因に関する変数
- 因子分析でいうところの「因子」に相当。SEMでは「構成概念」と呼ばれる。
- 例: 数学の能力、ブランド力 etc.
パス図では観測変数は四角形で、潜在変数は円形で表現されます。
構造方程式と測定方程式
因子分析では説明変数と因子の関係性しか記述しませんが、SEMでは因子（構成概念）同士の関係性も記述することができます。これらの関係性の記述は、事前の研究結果などから導かれた仮説に基づきます。
- 「構造方程式」 ... 仮説に基づいた因果関係を記述する方程式
- 「測定方程式」 ... 潜在変数とそれを構成する観測変数の関係性を記述する方程式
■ パス図による表現
semopyで構造方程式モデリング
pip install semopy で semopy パッケージをインストールできますが、Python3.7以外では正常に動作しない可能性があります。筆者は pyenv (python3.7.6) + jupyter-lab 環境で実行しています。
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 | 一人当たり |
| x3 | 産業における労働力の割合(%) (1960) |
表から以下のようなことが読みとれます。
-
y1~y4... 1960年の民主主義に関する観測変数 -->dem60 -
y5~y8... 1965年の民主主義に関する観測変数 -->dem65 -
x1~x3... 1960年の各種産業指数 -->ind60
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型で記述します。サンプルデータセットに対応するモデル記述が既に用意されているので、今回はそちらを利用します。パラメータ推定の結果、モデルに対するデータの当てはまりが良くない場合には、このモデル記述を編集することになります。
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の定義する構文に基づいています。詳細は公式ドキュメントをご参照ください。
-
~: 回帰演算子 --> 測定方程式を記述する -
=~: 測定演算子 --> 構造方程式を記述する -
~~: 分散演算子
【Do】 モデルを定義し、パラメータを推定する
semopy.Model に引数としてモデル記述を指定し、モデルを定義します。fit() メソッドを呼び出してパラメータを推定します。
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() メソッドで、パラメータ推定結果を確認できます。
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パッケージが必要になります)
In
semopy.semplot(model, "graph.png")
【Check & Action】 データの当てはまりの良さを確認する
calc_stats メソッドにより、各種適合度指標を確認できます。指標をもとに当てはまりが良くないと判断した場合、モデルの修正を検討します。
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
以上です。