パス解析#
モデルのすべての変数が観測変数であり、複数の説明変数から複数の目的変数を説明するようなモデル。
Sewall Wright （シューアル・ライト）という統計遺伝学者が1920年ごろに考案した。
最初のパス解析：Wright, S. (1920). The relative importance of heredity and environment in determining the piebald pattern of guinea-pigs
逐次モデル#
逐次モデルは内生変数間の双方向のパスが無いモデル
この図のモデルは
\[\begin{split}
v_3 = a_{31} v_1 + a_{32} v_2 + e_3\\
v_4 = a_{42} v_2 + a_{43} v_3 + e_4
\end{split}\]
という同時方程式（構造方程式）を解くことになる。
行列形式で
\[
\newcommand{\b}[1]{\boldsymbol{#1}}
\b{v} = \b{A v} + \b{e}
\]
のような形で書くと次のようになる。
\[\begin{split}
\begin{bmatrix}
v_1 \\ v_2 \\ v_3 \\ v_4
\end{bmatrix}
=
\begin{bmatrix}
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
a_{31} & a_{32} & 0 & 0\\
0 & a_{42} & a_{43} & 0\\
\end{bmatrix}
\begin{bmatrix}
v_1 \\ v_2 \\ v_3 \\ v_4
\end{bmatrix}
+
\begin{bmatrix}
v_1 \\ v_2 \\ e_3 \\ e_4
\end{bmatrix}
\end{split}\]
構造方程式には次のような性質がある
構造方程式の性質
\(\b{A}\)の対角成分は常に0
\(v_i\)が外生変数であれば、\(\b{A}\)の\(i\)行はゼロベクトル\(\b{o}^T\)である
\(v_i\)が外生変数であれば、残差ベクトルの要素\(e_i\)は観測値\(v_i\)になる
共分散の構造化#
観測変数の共分散\(\b{\Sigma}\)は
\[\begin{split}
\b{\Sigma}
= E[\b{v}\b{v}^T]
=
\begin{bmatrix}
E[v^2_1]\\
E[v_2 v_1] & E[v_2^2]\\
E[v_3 v_1] & E[v_3 v_2] & E[v_3^2]\\
E[v_4 v_1] & E[v_4 v_2] & E[v_4 v_3] & E[v_4^2]
\end{bmatrix}
\end{split}\]
で、これらの各要素は
\[\begin{split}
\begin{align}
E[v_3 v_1] &= E[(a_{31} v_1 + a_{32} v_2 + e_3) v_1]\\
&= a_{31} E[v_1^2] + a_{32} E[v_1 v_2] + E[e_3 v_1]\\
&= a_{31} \sigma^2_1 + a_{32} \sigma_{12}
\end{align}
\end{split}\]
\[\begin{split}
\begin{align}
E[v_3 v_2] &= E[(a_{31} v_1 + a_{32} v_2 + e_3) v_2]\\
&= a_{31} \sigma_{12} + a_{32} \sigma^2_2
\end{align}
\end{split}\]
\[\begin{split}
\begin{align}
E[v_3 v_3]
&= E[(a_{31} v_1 + a_{32} v_2 + e_3)^2]\\
&= E[(a_{31} v_1 + a_{32} v_2 + e_3)(a_{31} v_1 + a_{32} v_2 + e_3)]\\
&= E[
(a_{31} v_1)^2
+ 2 (a_{31} v_1) (a_{32} v_2)
+ (a_{31} v_1) e_3
+ (a_{32} v_2)^2
+ (a_{32} v_2) e_3
+ (a_{31} v_1) e_3
+ (a_{32} v_2) e_3
+ e_3^2
]\\
&= a_{31}^2 E[v_1^2]
+ 2 a_{31} a_{32} E[v_1 v_2]
+ a_{32}^2 E[v_2^2]
+ E[e_3^2]
\\
&= a_{31}^2 \sigma^2_2
+ 2 a_{31} a_{32} \sigma_{12}
+ a_{32}^2 \sigma^2_2
+ \sigma^2_{e3}
\end{align}
\end{split}\]
\[\begin{split}
\begin{align}
E[v_4 v_1] &= E[(a_{42} v_2 + a_{43} v_3 + e_4) v_1]\\
&= a_{42} \sigma_{21} + a_{43} \sigma_{31}
\end{align}
\end{split}\]
\[\begin{split}
\begin{align}
E[v_4 v_2] &= E[(a_{42} v_2 + a_{43} v_3 + e_4) v_2]\\
&= a_{42} \sigma^2_2 + a_{43} \sigma_{32}
\end{align}
\end{split}\]
\[\begin{split}
\begin{align}
E[v_4 v_3]
&= E[(a_{42} v_2 + a_{43} v_3 + e_4)
(a_{31} v_1 + a_{32} v_2 + e_3)]\\
&= a_{42} a_{31} \sigma_{21}
+ a_{42} a_{32} \sigma^2_2
+ a_{43} a_{31} \sigma_{31}
+ a_{43} a_{32} \sigma_{32}
\end{align}
\end{split}\]
のような形で、母数によって説明される形になっている。
行列を用いた共分散の構造化#
これを行列を用いて表記すると、
\[
\b{v} = \b{Av} + \b{e}
\]
の式を
\[\begin{split}
\begin{align}
\b{I v} &= \b{A v} + \b{e}\\
(\b{I} - \b{A})\b{v} &= \b{e}
\end{align}
\end{split}\]
と変形し、\((\b{I} - \b{A})\)に逆行列が存在することを仮定して
\[
\b{T} = (\b{I} - \b{A})^{-1}
\]
とおけば、構造方程式は
\[
\b{v} = \b{T e}
\]
と簡潔に表記できる。
観測変数の共分散構造は
\[\begin{split}
\begin{align}
\b{\Sigma}
&= E[\b{v v}^T]
= E[\b{Te}(\b{Te})^T]
= E[\b{Te}\b{e}^T \b{T}^T]
= \b{T} E[\b{ee}^T] \b{T}^T\\
&= \b{T \Sigma_e T}^T
\end{align}
\end{split}\]
となる。
残差ベクトル\(\b{e}\)の共分散行列\(\b{\Sigma}_{\b{e}}\)は
\[\begin{split}
\b{\Sigma_e}
=
\begin{bmatrix}
E[v^2_1]\\
E[v_2 v_1] & E[v_2^2]\\
E[e_3 v_1] & E[e_4 v_2] & E[e_3^2]\\
E[e_4 v_1] & E[e_4 v_2] & E[e_4 e_3] & E[e_4^2]
\end{bmatrix}
=
\begin{bmatrix}
\sigma^2_1\\
\sigma_{21} & \sigma^2_2\\
0 & 0 & \sigma^2_{e3}\\
0 & 0 & 0 & \sigma^2_{e4}
\end{bmatrix}
\end{split}\]
となり、共分散構造は\(\b{A}\)と\(\b{\Sigma_e}\)が特定されれば一意に定まることがわかる。
データを使って推定#
import numpy as np
import pandas as pd
# データ生成
n = 1000
np.random.seed(0)
# 真のパラメータ
A = np.array([
[0, 0, 0, 0],
[0, 0, 0, 0],
[3, 2, 0, 0],
[0, 5, 4, 0],
])
# NOTE: eは実際には観測できない。推定前に得らるのはvのみ
v1 = np.random.uniform(size=n)
v2 = np.random.uniform(size=n)
e3 = np.random.normal(size=n)
e4 = np.random.normal(size=n)
v3 = A[2, 0] * v1 + A[2, 1] * v2 + e3
v4 = A[3, 1] * v2 + A[3, 2] * v3 + e4
data = pd.DataFrame(dict(v1=v1, v2=v2, v3=v3, v4=v4))
# 標本共分散
V = np.array([v1, v2, v3, v4]).T
S = np.cov(V, rowvar=False)
S.round(2)
array([[ 0.08, 0. , 0.25, 0.99],
[ 0. , 0.09, 0.17, 1.13],
[ 0.25, 0.17, 1.99, 8.82],
[ 0.99, 1.13, 8.82, 42.01]])
目的関数の設定
\(\b{\theta} = (a_{11}, a_{12}, \dots, a_{pp}, \sigma_{e1}^2, \sigma_{e1 e2}, \dots, \sigma^2_{ep})\)として構造化した共分散を関数とする
\[
\b{\Sigma}(\b{\theta}) = \b{T \Sigma_e T}^T
\]
def calc_sigma(A, Sigma_e):
p = A.shape[1]
I = np.eye(p, p)
T = np.linalg.inv(I - A)
return T @ Sigma_e @ T.T
負の対数尤度
\[
f_{ML} = \text{tr}(\b{\Sigma}(\b{\theta})^{-1} \b{S})
- \log | \b{\Sigma}(\b{\theta})^{-1} \b{S}| - n
\]
def log_likelihood(Sigma, S, n):
SigmaS = np.linalg.inv(Sigma) @ S
return np.trace( SigmaS ) - np.log( np.linalg.det(SigmaS) ) - n
p = S.shape[0]
A_hat = np.random.uniform(size=(p, p)) # 初期値
Sigma_e_hat = np.random.uniform(size=(p, p)) # 初期値
Sigma_hat = calc_sigma(A_hat, Sigma_e_hat)
log_likelihood(Sigma_hat, S, n=n)
/tmp/ipykernel_7124/594764428.py:3: RuntimeWarning: invalid value encountered in log
return np.trace( SigmaS ) - np.log( np.linalg.det(SigmaS) ) - n
nan
import scipy.linalg.lapack as lapack
x = V.copy()
c, info = lapack.dpotrf(x)
lapack.dpotri(c, overwrite_c=1)
c + c.T - np.diag(c.diagonal())
---------------------------------------------------------------------------
error Traceback (most recent call last)
Cell In[8], line 3
1 import scipy.linalg.lapack as lapack
2 x = V.copy()
----> 3 c, info = lapack.dpotrf(x)
4 lapack.dpotri(c, overwrite_c=1)
5 c + c.T - np.diag(c.diagonal())
error: (shape(a,0)==shape(a,1)) failed for 1st argument a
import semopy
desc = """
v3 ~ v1 + v2
v4 ~ v2 + v3
"""
model = semopy.Model(desc)
model.fit(data)
model.inspect(std_est=True).round(2)
| lval | op | rval | Estimate | Est. Std | Std. Err | z-value | p-value | |
|---|---|---|---|---|---|---|---|---|
| 0 | v3 | ~ | v1 | 2.93 | 0.60 | 0.11 | 27.65 | 0.0 |
| 1 | v3 | ~ | v2 | 1.87 | 0.40 | 0.10 | 18.18 | 0.0 |
| 2 | v4 | ~ | v2 | 5.08 | 0.23 | 0.11 | 44.84 | 0.0 |
| 3 | v4 | ~ | v3 | 4.00 | 0.87 | 0.02 | 166.72 | 0.0 |
| 4 | v3 | ~~ | v3 | 0.95 | 0.48 | 0.04 | 22.36 | 0.0 |
| 5 | v4 | ~~ | v4 | 0.96 | 0.02 | 0.04 | 22.36 | 0.0 |
semopyはRの{lavaan}とは異なり、外生変数間の共変関係や分散を推定してくれない（そのため自由度が上がり適合度系の指標がおかしくなったりする）
なので
x1 ~~ x2
x1 ~~ x1
x2 ~~ x2
も指定してやる必要がある
import semopy
import numpy as np
import pandas as pd
# データ生成
n = 1000
np.random.seed(0)
x1 = np.random.uniform(size=n)
x2 = np.random.uniform(size=n)
e1 = np.random.normal(size=n)
e2 = np.random.normal(size=n)
y1 = 10 + 3 * x1 + 5 * x2 + e1
y2 = 5 + 10 * x1 + 15 * x2 + e2
data = pd.DataFrame(dict(y1=y1, y2=y2, x1=x1, x2=x2))
# モデルを構築
desc = """
y1 ~ x1 + x2
y2 ~ x1 + x2
x1 ~~ x2
x1 ~~ x1
x2 ~~ x2
"""
model = semopy.Model(desc)
model.fit(data)
model.inspect(std_est=True)
| lval | op | rval | Estimate | Est. Std | Std. Err | z-value | p-value | |
|---|---|---|---|---|---|---|---|---|
| 0 | y1 | ~ | x1 | 2.925450 | 0.436100 | 0.105792 | 27.652921 | 0.000000 |
| 1 | y1 | ~ | x2 | 4.869091 | 0.746496 | 0.102865 | 47.334923 | 0.000000 |
| 2 | y2 | ~ | x1 | 9.907971 | 0.528034 | 0.106670 | 92.884242 | 0.000000 |
| 3 | y2 | ~ | x2 | 15.085703 | 0.826853 | 0.103719 | 145.448189 | 0.000000 |
| 4 | x1 | ~~ | x2 | 0.000515 | 0.005930 | 0.002748 | 0.187510 | 0.851261 |
| 5 | x1 | ~~ | x1 | 0.084481 | 1.000000 | 0.003778 | 22.360680 | 0.000000 |
| 6 | x2 | ~~ | x2 | 0.089357 | 1.000000 | 0.003996 | 22.360680 | 0.000000 |
| 7 | y1 | ~~ | y1 | 0.945464 | 0.248700 | 0.042282 | 22.360680 | 0.000000 |
| 8 | y2 | ~~ | y2 | 0.961229 | 0.032317 | 0.042987 | 22.360680 | 0.000000 |
# パス図
g = semopy.semplot(model, filename="/tmp/path_diagram.png", plot_covs=True)
# rank = same
with g.subgraph() as s:
s.attr(rank = "same")
s.node("x1")
s.node("x2")
g.attr(rankdir = "RL")
g.attr(pad = "0.2") # padding
g
print(g.source)
digraph G {
overlap=scale splines=true
edge [fontsize=12]
node [fillcolor="#cae6df" shape=circle style=filled]
node [shape=box style=""]
y1 [label=y1]
y2 [label=y2]
x1 [label=x1]
x2 [label=x2]
x1 -> y1 [label="2.925\np-val: 0.00"]
x2 -> y1 [label="4.869\np-val: 0.00"]
x1 -> y2 [label="9.908\np-val: 0.00"]
x2 -> y2 [label="15.086\np-val: 0.00"]
x2 -> x1 [label="0.001\np-val: 0.85" dir=both style=dashed]
{
rank=same
x1
x2
}
rankdir=RL
pad=0.2
}
考察#
上記x1,x2は共変関係は考えるが一方向への関係（例えばx2→x1）は考えていない。そのため、もしデータ側にそういう関係があった場合は欠落変数バイアスが発生する様子。
非逐次モデル#
非逐次モデル（non-recursive model）は単方向の矢印だけをたどって起点となった元の変数に戻ることができるモデル、つまり、内生変数間の双方向のパスがあるモデルである。
非逐次モデルの推定には、双方向のパスを仮定する2つの変数のうち、一方の変数にだけ直接影響を与え、残りの変数には影響を与えない外生変数である道具的変数（instrumental variable）が必要である。
\[\begin{split}
\begin{align}
v_2 &= a_{21} v_1 + a_{23} v_3 + e_2\\
v_3 &= a_{32} v_2 + e_3
\end{align}
\end{split}\]
\[
\newcommand{\b}[1]{\boldsymbol{#1}}
\b{v} = \b{A v} + \b{e}
\]
\[\begin{split}
\begin{bmatrix}
v_1 \\ v_2 \\ v_3
\end{bmatrix}
=
\begin{bmatrix}
0 & 0 & 0 \\
a_{21} & 0 & a_{23} \\
0 & a_{32} & 0
\end{bmatrix}
\begin{bmatrix}
v_1 \\ v_2 \\ v_3
\end{bmatrix}
+
\begin{bmatrix}
v_1 \\ e_2 \\ e_3
\end{bmatrix}
\end{split}\]
\[\begin{split}
\b{\Sigma}
= E[\b{v}\b{v}^T]
=
\begin{bmatrix}
E[v^2_1]\\
E[v_2 v_1] & E[v_2^2]\\
E[v_3 v_1] & E[v_3 v_2] & E[v_3^2]
\end{bmatrix}
\end{split}\]
\[\begin{split}
\begin{align}
E[v_1^2] &= \sigma^2_1\\
E[v_2 v_1] &= E[(a_{21} v_1 + a_{23} v_3 + e_2) v_1]\\
&= a_{21} \sigma^2_1 + a_{23} \sigma_{31}\\
E[v_2^2] &= E[(a_{21} v_1 + a_{23} v_3 + e_2)
(a_{21} v_1 + a_{23} v_3 + e_2)]\\
&= a_{21}^2 \sigma^2_1
+ 2 a_{21} a_{23} \sigma_{13}
+ a_{23}^2 \sigma^2_3
+ \sigma_{e2}^2
\\
E[v_3 v_1]
&= E[(a_{32} v_2 + e_3) v_1]\\
&= a_{32} \sigma_{21}\\
E[v_3 v_2]
&= E[(a_{32} v_2 + e_3) (a_{21} v_1 + a_{23} v_3 + e_2)]\\
&= a_{32} a_{21} \sigma_{21}
+ a_{32} a_{23} \sigma_{23}\\
E[v_3^2]
&= E[(a_{32} v_2 + e_3)^2]\\
&= a_{32}^2 \sigma^2_2
+ \sigma^2_{e3}
\end{align}
\end{split}\]
ゆえに
\[\begin{split}
\b{\Sigma}
=
\begin{bmatrix}
\sigma^2_1\\
a_{21} \sigma^2_1 + a_{23} \sigma_{31}
& a_{21}^2 \sigma^2_1
+ 2 a_{21} a_{23} \sigma_{13}
+ a_{23}^2 \sigma^2_3
+ \sigma_{e2}^2
\\
a_{32} \sigma_{21}
& a_{32} a_{21} \sigma_{21} + a_{32} a_{23} \sigma_{23}
& a_{32}^2 \sigma^2_2 + \sigma^2_{e3}
\end{bmatrix}
\end{split}\]
\[\begin{split}
\begin{align}
\sigma^2_1 &= \sigma_1^2\\
\sigma_{21} &= a_{21} \sigma^2_1 + a_{23} \sigma_{31}\\
\sigma^2_2 &= a_{21}^2 \sigma^2_1
+ 2 a_{21} a_{23} \sigma_{13}
+ a_{23}^2 \sigma^2_3
+ \sigma_{e2}^2\\
\sigma_{31} &= \sigma_{32} \sigma_{21}\\
\sigma_{32} &= a_{32} a_{21} \sigma_{21} + a_{32} a_{23} \sigma_{23}\\
\sigma_3^2 &= a_{32}^2 \sigma^2_2 + \sigma^2_{e3}
\end{align}
\end{split}\]
考察：2SLS推定量とのつながり#
IV推定量#
同時方程式モデル
\[\begin{split}
Y_1 = \beta_1 Y_2 + \gamma_1 Z_1 + \varepsilon_1\\
Y_2 = \beta_2 Y_1 + \gamma_2 Z_2 + \varepsilon_2
\end{split}\]
があるとする。外生変数が\(Z\)である。
行列形式で\(\boldsymbol{\beta}_1 = (\beta_1, \gamma_1)^T\)のようにまとめて
\[\begin{split}
\newcommand{\b}[1]{\boldsymbol{#1}}
\b{Y}_1 = [\b{Y}_2|\b{Z}_1] \b{\beta}_1 + \b{\varepsilon}_1 = \b{X}_1 \b{\beta}_1 + \b{\varepsilon}_1\\
\b{Y}_2 = [\b{Y}_1|\b{Z}_2] \b{\beta}_2 + \b{\varepsilon}_2 = \b{X}_2 \b{\beta}_2 + \b{\varepsilon}_2
\end{split}\]
と書くことにする。外生変数からなる行列を\(\b{Z} = [\b{Z}_1|\b{Z}_2]\)と表すと
\[\]
MIMICモデル#
60年代にGoldbergerなど経済学者や社会科学者がパス解析を再発見した（『因果推論の科学』, p.159）
Goldberger & Karl Jöreskog (1975) はMIMIC modelを提案
参考文献#
パス解析
豊田秀樹（1998）『共分散構造分析 入門編』
豊田秀樹（2014）『共分散構造分析 R編』
非逐次モデル