【Pythonで統計モデル】SEM2：母数に関する制約
はじめに
株式会社GA technologiesの福中です。今、僕は、Advanced Innovation Strategy Center（以下AISC）という部署でChief Data Scientistとして働いています（詳しい自己紹介はAISC WEBサイトのプロフィールをご覧ください）。前回の記事で、semopyを使った構造方程式モデリング（SEM）の分析の流れを一通り説明しました。これにより（簡単なモデルであれば）PythonでもSEMが問題なく実行できることがわかりました。
さて、今回も前回に続きSEMの基礎的な話題を取り上げたいと思います。SEMは柔軟な手法であるとよく言われますが、その柔軟性を支えているものの１つに母数の制約があります。母数を推定する際に、SEMでは等値制約や不等式制約を比較的簡単に導入することが可能です。今回の記事では、これらの母数制約をsemopyで実行する方法について解説しようと思います。
本記事は前回の記事の補足的内容です。前回の記事を未読の方は先にそちらに目を通すことをお勧めします。
注意点
本記事は統計学およびPythonの初心者向け解説記事ではありません。普段からR等を使って統計モデルによる分析や研究を行っている研究者の方が、改めてPythonで構造方程式モデルを行う必要が出た際のリファレンスとしての利用を想定しています。したがって、SEMに不慣れな方は、先に以下で紹介するテキストをはじめとした初心者用の教科書で勉強するとよいと思います。
テキスト
本記事では東京図書から出版されている共分散構造分析［Ｒ編］―構造方程式モデリング―で取り上げられていた事例を、Pythonを使って再現していきます。分析用データは東京図書のホームページからダウンロードできますので、興味のある方はぜひ試してみてください。
環境
本記事の分析は、OS: Windows 11、Python 3.11.2の環境下で行っています。また、必要なパッケージは以下となります。インストール済みの方は割愛してください。
# packageのインストール
pip install numpy pandas
pip install graphviz
pip install semopyなおsemopyは構造方程式モデリング（SEM ）を行うためのパッケージで、ドキュメントは以下にあります。
本記事ではここから重要な部分をいくつかピックアップしてご紹介します。網羅的な解説はしませんので、適宜、必要に応じてご参照ください。
基本モデル
母数の制約に関する説明をするにあたり、まず基本となるモデルを以下で紹介します。今回も前回の記事で利用した「セミナーデータ」を利用し、以下のような2因子間に共分散を仮定した確認的因子分析モデルを基本モデルにしたいと思います。
パッケージとデータの読み込み
# packageの読み込み
import numpy as np
import pandas as pd
import semopy# データの読み込と列名の変更
X = pd.read_csv("dat/seminar.csv", encoding='shift_jis', index_col=0)
X.columns = ['x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7']分析の実行
# モデル表現
mod = """
f1 =~ x1 + x2 + x3 + x4
f2 =~ x5 + x6 + x7
"""# 分析の実行
model = semopy.Model(mod)
result = model.fit(data=X, obj='MLW')
inspect = semopy.Model.inspect(model)
print(inspect) lval op rval Estimate Std. Err z-value p-value
0 x1 ~ f1 1.000000 - - -
1 x2 ~ f1 2.251321 0.839931 2.680364 0.007354
2 x3 ~ f1 1.549742 0.615549 2.517659 0.011814
3 x4 ~ f1 1.858816 0.682278 2.724425 0.006441
4 x5 ~ f2 1.000000 - - -
5 x6 ~ f2 0.569675 0.249676 2.281657 0.02251
6 x7 ~ f2 0.385190 0.170429 2.260116 0.023814
7 f2 ~~ f2 0.866040 0.412014 2.101966 0.035556
8 f2 ~~ f1 0.137753 0.068096 2.022937 0.04308
9 f1 ~~ f1 0.106834 0.071119 1.502184 0.133049
10 x4 ~~ x4 0.904725 0.149414 6.055177 0.0
11 x3 ~~ x3 1.160349 0.167377 6.93256 0.0
12 x6 ~~ x6 1.462827 0.22986 6.363992 0.0
13 x2 ~~ x2 0.274503 0.131052 2.094616 0.036205
14 x5 ~~ x5 1.127231 0.381428 2.955291 0.003124
15 x1 ~~ x1 0.942659 0.128498 7.335975 0.0
16 x7 ~~ x7 0.715701 0.110717 6.46421 0.0母数制約の種類と方法
固定母数
推定すべき母数を何らかの値で固定したい場合、基本的に「定数*」という形で指定します。例えば、因子間共分散を0に固定したい場合、以下のように記述します。
# 固定母数：因子間共分散を０に固定する
mod = """
f1 =~ x1 + x2 + x3 + x4
f2 =~ x5 + x6 + x7
f2 ~~ 0*f1
"""本来、因子間の共分散（f2 ~~ f1）は自動で仮定されるので、シンタックスとしての記述は不要ですが、母数を何らかの定数で固定したい場合のみ、上記のように明示的に記述します。それではこのモデルで分析を実行してみましょう。
# 分析の実行
model = semopy.Model(mod)
result = model.fit(data=X, obj='MLW')
inspect = semopy.Model.inspect(model)
print(inspect)
lval op rval Estimate Std. Err z-value p-value
0 x1 ~ f1 1.000000 - - -
1 x2 ~ f1 2.085328 0.792655 2.630814 0.008518
2 x3 ~ f1 1.492594 0.58687 2.543312 0.010981
3 x4 ~ f1 1.829162 0.659734 2.772574 0.005561
4 x5 ~ f2 1.000000 - - -
5 x6 ~ f2 0.415361 0.297886 1.394361 0.163209
6 x7 ~ f2 0.331696 0.236956 1.399819 0.161567
7 f2 ~~ f1 0.000000 - - -
8 f2 ~~ f2 1.116831 0.818 1.365318 0.172153
9 f1 ~~ f1 0.117060 0.076279 1.534634 0.124874
10 x4 ~~ x4 0.881798 0.158748 5.554695 0.0
11 x3 ~~ x3 1.156754 0.170345 6.790675 0.0
12 x6 ~~ x6 1.551578 0.243113 6.38212 0.0
13 x2 ~~ x2 0.307246 0.143773 2.137017 0.032597
14 x5 ~~ x5 0.875276 0.792358 1.104647 0.269312
15 x1 ~~ x1 0.932532 0.128938 7.232402 0.0
16 x7 ~~ x7 0.721098 0.127499 5.65573 0.0結果を見るとわかるように7行目で因子間共分散が0に固定された上でその他の母数が推定されていることがわかります。
7 f2 ~~ f1 0.000000 - - -同様にして、因子の分散を1に固定してみましょう。モデル表現のスクリプトは以下のようになります。
# 固定母数：因子の分散を1に固定する
mod = """
f1 =~ x1 + x2 + x3 + x4
f2 =~ x5 + x6 + x7
f1 ~~ 1*f1
f2 ~~ 1*f2
"""# 分析の実行
model = semopy.Model(mod)
result = model.fit(data=X, obj='MLW')
inspect = semopy.Model.inspect(model)
print(inspect) lval op rval Estimate Std. Err z-value p-value
0 x1 ~ f1 1.000000 - - -
1 x2 ~ f1 0.700510 0.091402 7.66402 0.0
2 x3 ~ f1 0.663329 0.120383 5.510168 0.0
3 x4 ~ f1 0.708048 0.113378 6.245039 0.0
4 x5 ~ f2 1.000000 - - -
5 x6 ~ f2 0.520331 0.15839 3.285121 0.001019
6 x7 ~ f2 0.353281 0.110019 3.211083 0.001322
7 f1 ~~ f1 1.000000 - - -
8 f2 ~~ f2 1.000000 - - -
9 f2 ~~ f1 0.478001 0.125394 3.811986 0.000138
10 x4 ~~ x4 0.919521 0.144069 6.382493 0.0
11 x3 ~~ x3 1.106764 0.164143 6.74269 0.0
12 x6 ~~ x6 1.489731 0.221346 6.730324 0.0
13 x2 ~~ x2 0.469468 0.091095 5.153587 0.0
14 x5 ~~ x5 1.038109 0.237349 4.373772 0.000012
15 x1 ~~ x1 0.839592 0.159827 5.253135 0.0
16 x7 ~~ x7 0.726772 0.107117 6.784847 0.07行目と8行目を見るとわかるように、因子の分散がきちんと1に固定されていることがわかります。
固定母数の制約でもっともよく利用されるシチュエーションの1つが、測定方程式におけるパス係数を1に固定する際の変数の変更です。SEMでは母数を推定するために、潜在変数から観測変数へのパス係数のうちどれか1つを必ず1に固定しなければなりません。semopyでは1つ目の観測変数へのパス係数が自動的に1に固定されるわけですが、1番目以外のパス係数を1にしたい場合などがあります。その場合、1番目のパス係数を「NA*」とし、1に固定したいパス係数を「1*」にします。
# 固定母数：2番目の因子負荷を1に固定する
mod = """
f1 =~ NA*x1 + 1*x2 + x3 + x4
f2 =~ x5 + x6 + x7
"""# 分析の実行
model = semopy.Model(mod)
result = model.fit(data=X, obj='MLW')
inspect = semopy.Model.inspect(model)
print(inspect) lval op rval Estimate Std. Err z-value p-value
0 x1 ~ f1 0.445319 0.165791 2.686021 0.007231
1 x2 ~ f1 1.000000 - - -
2 x3 ~ f1 0.690279 0.208961 3.303385 0.000955
3 x4 ~ f1 0.825986 0.223251 3.699801 0.000216
4 x5 ~ f2 1.000000 - - -
5 x6 ~ f2 0.569875 0.24958 2.283341 0.02241
6 x7 ~ f2 0.385120 0.170305 2.261357 0.023737
7 f2 ~~ f2 0.866258 0.411829 2.10344 0.035427
8 f2 ~~ f1 0.310147 0.113665 2.728615 0.00636
9 f1 ~~ f1 0.540994 0.160669 3.367138 0.00076
10 x4 ~~ x4 0.904631 0.149327 6.058062 0.0
11 x3 ~~ x3 1.159880 0.167401 6.928768 0.0
12 x6 ~~ x6 1.462325 0.229807 6.363287 0.0
13 x2 ~~ x2 0.274729 0.130748 2.101211 0.035622
14 x5 ~~ x5 1.126166 0.381197 2.954292 0.003134
15 x1 ~~ x1 0.942256 0.128472 7.334311 0.0
16 x7 ~~ x7 0.715629 0.110685 6.465425 0.00行目と1行目を見るとわかるように、固定母数の制約がきちんと作用していることが見て取れます。
母数のラベルと初期値の設定
複雑なモデルを構築すると、所定の繰り返し計算の範囲内で解が求まらないようなことがたまにあります。その際、初期値を適切に設定することにより、解が求まることがあるので、初期値を設定する方法を知っておくことは重要です。
semopyで初期値を設定するためには、最初にモデル表現で推定したいパラメータにラベル（名前）を付けます。例えば、「f2 =~ x5 + a*x6 + b*x7」のようにします。この場合、x6のパス係数をa、x7のパス係数をbとネーミングしていることになります。
次に、初期値を設定するために、「START(1.5) a b」のような内容を付け加えます。この場合、母数aとbの初期値を1.5に設定し、そこから母数の推定を開始するということを意味しています。
スクリプト全体は以下のようになります。なお、今回のような単純なモデルではどのように初期値を設定しても唯一解に行きつくので、結果の出力は省略します。
# 母数のラベルと初期値の設定
mod = """
f1 =~ x1 + x2 + x3 + x4
f2 =~ x5 + a*x6 + b*x7
START(1.5) a b
"""母数のラベルと等値制約
同じモデルの中で異なるパラメータを同じ値になるように推定した場合があります。これを等値制約と呼びます。semopyで等値制約を行いたい場合、同じラベルを振ることで実現できます。例えば以下のようなスクリプトを書いたとします。
# 母数のラベルと等値制約
mod = """
f1 =~ x1 + a*x2 + b*x3 + x4
f2 =~ x5 + a*x6 + b*x7
START(1.5) a b
"""この場合、x2とx6の変数、x3とx7の変数のパス係数が同じ値として推定されます。実行してみましょう。
# 分析の実行
model = semopy.Model(mod)
result = model.fit(data=X, obj='MLW')
inspect = semopy.Model.inspect(model)
print(inspect) lval op rval Estimate Std. Err z-value p-value
0 x1 ~ f1 1.000000 - - -
1 x2 ~ f1 1.502471 0.391266 3.840027 0.000123
2 x6 ~ f2 1.502471 0.391266 3.840027 0.000123
3 x3 ~ f1 0.918290 0.276962 3.315588 0.000915
4 x7 ~ f2 0.918290 0.276962 3.315588 0.000915
5 x4 ~ f1 1.373243 0.382919 3.586245 0.000335
6 x5 ~ f2 1.000000 - - -
7 f2 ~~ f2 0.212298 0.108285 1.96055 0.049932
8 f2 ~~ f1 0.102010 0.055464 1.839224 0.065882
9 f1 ~~ f1 0.210277 0.093535 2.248111 0.024569
10 x4 ~~ x4 0.876928 0.150148 5.84042 0.0
11 x3 ~~ x3 1.175238 0.162901 7.214439 0.0
12 x6 ~~ x6 1.305105 0.246207 5.300855 0.0
13 x2 ~~ x2 0.330829 0.112703 2.935396 0.003331
14 x5 ~~ x5 1.638628 0.235468 6.959031 0.0
15 x1 ~~ x1 0.918663 0.132792 6.918058 0.0
16 x7 ~~ x7 0.703963 0.116594 6.037752 0.01行目から4行目を見るとわかるように、きちんと等値制約が行われています。
母数のラベルと推定値の境界の設定
特定の区間内で母数推定を行いたい場合、BOUND(l, r)を利用します。lには推定したい値の下限を、rには推定したい値の上限を設定します。スクリプト自体は以下のように記述します。
# 母数のラベルと推定値の境界の設定
mod = """
f1 =~ x1 + x2 + x3 + x4
f2 =~ x5 + x6 + x7
f1 ~~ a*f2
BOUND(0.5, 1) a
"""# 分析の実行
model = semopy.Model(mod)
result = model.fit(data=X, obj='MLW')
inspect = semopy.Model.inspect(model)
print(inspect) lval op rval Estimate Std. Err z-value p-value
0 x1 ~ f1 1.000000 - - -
1 x2 ~ f1 1.144927 0.24335 4.704866 0.000003
2 x3 ~ f1 0.941991 0.246288 3.824758 0.000131
3 x4 ~ f1 1.077377 0.250948 4.293234 0.000018
4 x5 ~ f2 1.000000 - - -
5 x6 ~ f2 0.434258 0.146192 2.970456 0.002974
6 x7 ~ f2 0.286917 0.099642 2.87947 0.003983
7 f1 ~~ f2 0.500000 0.139435 3.585898 0.000336
8 f1 ~~ f1 0.388710 0.14006 2.775322 0.005515
9 f2 ~~ f2 1.513065 0.49957 3.028738 0.002456
10 x4 ~~ x4 0.900705 0.143293 6.285739 0.0
11 x3 ~~ x3 1.132169 0.16516 6.854982 0.0
12 x6 ~~ x6 1.507968 0.216459 6.96653 0.0
13 x2 ~~ x2 0.394916 0.094788 4.166298 0.000031
14 x5 ~~ x5 0.897492 0.422358 2.124953 0.033591
15 x1 ~~ x1 0.882054 0.136455 6.464058 0.0
16 x7 ~~ x7 0.741024 0.104965 7.05971 0.0制約を行わない場合のf1とf2の因子間相関は0.137753でしたが、7行目を見るとわかるように、境界制約を加えることで下限の0.5になっていることが見て取れます。
母数の等式制約
ある母数が他の母数の関数になっているような場合、その方程式をモデル表現に書き込むことでその制約を実現できます。例えば以下のようなスクリプトを書いたとします。
# 母数の等式制約
mod = """
f1 =~ x1 + a*x2 + b*x3 + c*x4
f2 =~ x5 + x6 + x7
CONSTRAINT(a = (b + c)^2)
"""この場合、パス係数のbとcの和の2乗がx2のパス係数となるように母数の推定を行います。普通に分析する範囲では使うシチュエーションは少ないかもしれませんが、SEMがいかに柔軟な統計モデルであるかの一端が見えるでしょう。
# 分析の実行
model = semopy.Model(mod)
result = model.fit(data=X, obj='MLW')
inspect = semopy.Model.inspect(model)
print(inspect) lval op rval Estimate Std. Err z-value p-value
0 x1 ~ f1 1.000000e+00 - - -
1 x2 ~ f1 2.568319e+00 1.294596 1.983877 0.04727
2 x3 ~ f1 6.104322e-01 0.337581 1.808255 0.070567
3 x4 ~ f1 9.921656e-01 0.373636 2.655437 0.007921
4 x5 ~ f2 1.000000e+00 - - -
5 x6 ~ f2 5.649799e-01 0.249581 2.26371 0.023592
6 x7 ~ f2 4.148107e-01 0.180137 2.30275 0.021293
7 f2 ~~ f2 8.354820e-01 0.400569 2.085737 0.037002
8 f2 ~~ f1 1.233020e-01 0.076978 1.601784 0.109203
9 f1 ~~ f1 1.236872e-01 0.085647 1.444155 0.148695
10 x4 ~~ x4 1.039224e+00 0.145057 7.164263 0.0
11 x3 ~~ x3 1.285376e+00 0.168509 7.627926 0.0
12 x6 ~~ x6 1.477424e+00 0.229228 6.445228 0.0
13 x2 ~~ x2 1.959251e-17 0.350537 0.0 1.0
14 x5 ~~ x5 1.157093e+00 0.372219 3.108638 0.00188
15 x1 ~~ x1 1.014130e+00 0.142322 7.125601 0.0
16 x7 ~~ x7 7.001871e-01 0.11209 6.246633 0.0分析の結果を見てみると、1行目の推定値（2.568319）は、2行目の推定値（0.6104322）と3行目の推定値（0.9921656）の和の2乗になっており、きちんと制約が行われていることがわかります。
母数の不等式制約
SEMでは母数の不等式制約も行えます。不等式制約とは2つのパラメータの一方が必ず大きくなるように推定する制約や、特定のパラメータが0などの特定の値より大きくなるように推定する制約のことです。使用するシチュエーションとしては、誤差分散がマイナスになってしまった場合（不適解）の対処として、その分散パラメータを0より大きくなるように制約を加えるなどが考えられます。
# 母数の不等式制約：x2の誤差分散を0より大きくなるように推定する
mod = """
f1 =~ x1 + x2 + x3 + x4
f2 =~ x5 + x6 + x7
x2 ~~ a*x2
CONSTRAINT(a > 0)
"""なお、このモデルではもともと誤差分散は適切に推定されているので、結果の出力は省略します。
使うシチュエーションはかなり限られますが、以下ではかなり複雑な不等式制約が行えることを示しています。
# 母数の不等式制約
mod = """
f1 =~ x1 + a*x2 + b*x3 + c*x4
f2 =~ x5 + x6 + x7
CONSTRAINT(exp(a) + log(b) = 10)
CONSTRAINT(c > cos(a)^2 + sin(b)^2)
"""# 分析の実行
model = semopy.Model(mod)
result = model.fit(data=X, obj='MLW')
inspect = semopy.Model.inspect(model)
print(inspect) lval op rval Estimate Std. Err z-value p-value
0 x1 ~ f1 1.000000 - - -
1 x2 ~ f1 2.257438 0.843349 2.676753 0.007434
2 x3 ~ f1 1.554932 0.617897 2.516492 0.011853
3 x4 ~ f1 1.860176 0.683767 2.720481 0.006519
4 x5 ~ f2 1.000000 - - -
5 x6 ~ f2 0.570112 0.249665 2.283512 0.0224
6 x7 ~ f2 0.385392 0.170406 2.261606 0.023722
7 f2 ~~ f2 0.865519 0.411544 2.103105 0.035457
8 f2 ~~ f1 0.137507 0.06801 2.021873 0.04319
9 f1 ~~ f1 0.106399 0.070944 1.499754 0.133678
10 x4 ~~ x4 0.905383 0.149346 6.06233 0.0
11 x3 ~~ x3 1.160057 0.167385 6.930491 0.0
12 x6 ~~ x6 1.462219 0.229801 6.362969 0.0
13 x2 ~~ x2 0.273809 0.131187 2.087172 0.036873
14 x5 ~~ x5 1.127131 0.380965 2.958623 0.00309
15 x1 ~~ x1 0.942630 0.128459 7.338009 0.0
16 x7 ~~ x7 0.715790 0.110723 6.464692 0.0このようにかなり複雑な制約でも実行することは可能ですが、解が求まることとその不等式制約が妥当であるか否かは別の問題であることに注意してください。
なお、semopyの等式制約・不等式制約で使用できる数式表現はsympyというパッケージに準拠しています。sympyはPythonの代数計算ライブラリで、MathematicaやMapleのような数式処理を行うためのパッケージです。興味のある方は以下のサイトをご参照ください。
この記事が気に入ったらサポートをしてみませんか？