パス解析

パス解析#

モデルのすべての変数が観測変数であり、複数の説明変数から複数の目的変数を説明するようなモデル。

Sewall Wright （シューアル・ライト）という統計遺伝学者が1920年ごろに考案した。

最初のパス解析：Wright, S. (1920). The relative importance of heredity and environment in determining the piebald pattern of guinea-pigs

逐次モデル#

逐次モデルは内生変数間の双方向のパスが無いモデル

../../_images/18cff059f20579b608ecb1b42914f24dbe4a00ea167d8005e4051df207a5d0cc.svg

この図のモデルは

\begin{array}{r} 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{array}

という同時方程式（構造方程式）を解くことになる。

行列形式で

v = A v + e

のような形で書くと次のようになる。

\begin{array}{r} [\begin{array}{c} v_{1} \\ v_{2} \\ v_{3} \\ v_{4} \end{array}] = [\begin{array}{c} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ a_{31} & a_{32} & 0 & 0 \\ 0 & a_{42} & a_{43} & 0 \end{array}] [\begin{array}{c} v_{1} \\ v_{2} \\ v_{3} \\ v_{4} \end{array}] + [\begin{array}{c} v_{1} \\ v_{2} \\ e_{3} \\ e_{4} \end{array}] \end{array}

構造方程式には次のような性質がある

構造方程式の性質

$A$ の対角成分は常に0
$v_{i}$ が外生変数であれば、 $A$ の $i$ 行はゼロベクトル $o^{T}$ である
$v_{i}$ が外生変数であれば、残差ベクトルの要素 $e_{i}$ は観測値 $v_{i}$ になる

共分散の構造化#

観測変数の共分散 $Σ$ は

\begin{array}{r} Σ = E [v v^{T}] = [\begin{array}{c} E [v_{1}^{2}] \\ 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{array}] \end{array}

で、これらの各要素は

\begin{array}{r} \begin{aligned} 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} σ_{1}^{2} + a_{32} σ_{12} \end{aligned} \end{array}

\begin{array}{r} \begin{aligned} E [v_{3} v_{2}] & = E [(a_{31} v_{1} + a_{32} v_{2} + e_{3}) v_{2}] \\ = a_{31} σ_{12} + a_{32} σ_{2}^{2} \end{aligned} \end{array}

\begin{array}{r} \begin{aligned} 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} σ_{2}^{2} + 2 a_{31} a_{32} σ_{12} + a_{32}^{2} σ_{2}^{2} + σ_{e 3}^{2} \end{aligned} \end{array}

\begin{array}{r} \begin{aligned} E [v_{4} v_{1}] & = E [(a_{42} v_{2} + a_{43} v_{3} + e_{4}) v_{1}] \\ = a_{42} σ_{21} + a_{43} σ_{31} \end{aligned} \end{array}

\begin{array}{r} \begin{aligned} E [v_{4} v_{2}] & = E [(a_{42} v_{2} + a_{43} v_{3} + e_{4}) v_{2}] \\ = a_{42} σ_{2}^{2} + a_{43} σ_{32} \end{aligned} \end{array}

\begin{array}{r} \begin{aligned} 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} σ_{21} + a_{42} a_{32} σ_{2}^{2} + a_{43} a_{31} σ_{31} + a_{43} a_{32} σ_{32} \end{aligned} \end{array}

のような形で、母数によって説明される形になっている。

行列を用いた共分散の構造化#

これを行列を用いて表記すると、

v = A v + e

の式を

\begin{array}{r} \begin{aligned} I v & = A v + e \\ (I - A) v & = e \end{aligned} \end{array}

と変形し、 $(I - A)$ に逆行列が存在することを仮定して

T = (I - A)^{- 1}

とおけば、構造方程式は

v = T e

と簡潔に表記できる。

観測変数の共分散構造は

\begin{array}{r} \begin{aligned} Σ & = E [v v^{T}] = E [T e (T e)^{T}] = E [T e e^{T} T^{T}] = T E [e e^{T}] T^{T} \\ = T Σ_{e} T^{T} \end{aligned} \end{array}

となる。

残差ベクトル $e$ の共分散行列 $Σ_{e}$ は

\begin{array}{r} Σ_{e} = [\begin{array}{c} E [v_{1}^{2}] \\ 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{array}] = [\begin{array}{c} σ_{1}^{2} \\ σ_{21} & σ_{2}^{2} \\ 0 & 0 & σ_{e 3}^{2} \\ 0 & 0 & 0 & σ_{e 4}^{2} \end{array}] \end{array}

となり、共分散構造は $A$ と $Σ_{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]])

目的関数の設定

$θ = (a_{11}, a_{12}, \dots, a_{p p}, σ_{e 1}^{2}, σ_{e 1 e 2}, \dots, σ_{e p}^{2})$ として構造化した共分散を関数とする

Σ (θ) = T Σ_{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_{M L} = tr (Σ (θ)^{- 1} S) - \log | Σ (θ)^{- 1} 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_15269/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
0	v3	~	v1	2.93	0.60	0.11	27.65
1	v3	~	v2	1.87	0.40	0.10	18.18
2	v4	~	v2	5.08	0.23	0.11	44.84
3	v4	~	v3	4.00	0.87	0.02	166.72
4	v3	~~	v3	0.95	0.48	0.04	22.36
5	v4	~~	v4	0.96	0.02	0.04	22.36

../../_images/873b47d1f60215f776f3a0f5fa166a7fddb930c6f70b6cd432b493c2932d6c57.svg

../../_images/0291abaabc44267eb38be04c561a87444339482c1e3db971970c5881cb59e139.svg

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

../../_images/6d59b64d203102ff3dc01e95b27240091c04d88112903d0dc91a4238cd0bd961.svg

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）が必要である。

../../_images/ab68ab8f022f7a95e9f08f38594196c2ec73ab64e294d66fdac88865fab53936.svg

\begin{array}{r} \begin{aligned} v_{2} & = a_{21} v_{1} + a_{23} v_{3} + e_{2} \\ v_{3} & = a_{32} v_{2} + e_{3} \end{aligned} \end{array}

v = A v + e

\begin{array}{r} [\begin{array}{c} v_{1} \\ v_{2} \\ v_{3} \end{array}] = [\begin{array}{c} 0 & 0 & 0 \\ a_{21} & 0 & a_{23} \\ 0 & a_{32} & 0 \end{array}] [\begin{array}{c} v_{1} \\ v_{2} \\ v_{3} \end{array}] + [\begin{array}{c} v_{1} \\ e_{2} \\ e_{3} \end{array}] \end{array}

\begin{array}{r} Σ = E [v v^{T}] = [\begin{array}{c} E [v_{1}^{2}] \\ E [v_{2} v_{1}] & E [v_{2}^{2}] \\ E [v_{3} v_{1}] & E [v_{3} v_{2}] & E [v_{3}^{2}] \end{array}] \end{array}

\begin{array}{r} \begin{aligned} E [v_{1}^{2}] & = σ_{1}^{2} \\ E [v_{2} v_{1}] & = E [(a_{21} v_{1} + a_{23} v_{3} + e_{2}) v_{1}] \\ = a_{21} σ_{1}^{2} + a_{23} σ_{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} σ_{1}^{2} + 2 a_{21} a_{23} σ_{13} + a_{23}^{2} σ_{3}^{2} + σ_{e 2}^{2} \\ E [v_{3} v_{1}] & = E [(a_{32} v_{2} + e_{3}) v_{1}] \\ = a_{32} σ_{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} σ_{21} + a_{32} a_{23} σ_{23} \\ E [v_{3}^{2}] & = E [(a_{32} v_{2} + e_{3})^{2}] \\ = a_{32}^{2} σ_{2}^{2} + σ_{e 3}^{2} \end{aligned} \end{array}

ゆえに

\begin{array}{r} Σ = [\begin{array}{c} σ_{1}^{2} \\ a_{21} σ_{1}^{2} + a_{23} σ_{31} & a_{21}^{2} σ_{1}^{2} + 2 a_{21} a_{23} σ_{13} + a_{23}^{2} σ_{3}^{2} + σ_{e 2}^{2} \\ a_{32} σ_{21} & a_{32} a_{21} σ_{21} + a_{32} a_{23} σ_{23} & a_{32}^{2} σ_{2}^{2} + σ_{e 3}^{2} \end{array}] \end{array}

\begin{array}{r} \begin{aligned} σ_{1}^{2} & = σ_{1}^{2} \\ σ_{21} & = a_{21} σ_{1}^{2} + a_{23} σ_{31} \\ σ_{2}^{2} & = a_{21}^{2} σ_{1}^{2} + 2 a_{21} a_{23} σ_{13} + a_{23}^{2} σ_{3}^{2} + σ_{e 2}^{2} \\ σ_{31} & = σ_{32} σ_{21} \\ σ_{32} & = a_{32} a_{21} σ_{21} + a_{32} a_{23} σ_{23} \\ σ_{3}^{2} & = a_{32}^{2} σ_{2}^{2} + σ_{e 3}^{2} \end{aligned} \end{array}

考察：2SLS推定量とのつながり#

IV推定量#

同時方程式モデル

\begin{array}{r} Y_{1} = β_{1} Y_{2} + γ_{1} Z_{1} + ε_{1} \\ Y_{2} = β_{2} Y_{1} + γ_{2} Z_{2} + ε_{2} \end{array}

があるとする。外生変数が $Z$ である。

行列形式で $β_{1} = (β_{1}, γ_{1})^{T}$ のようにまとめて

\begin{array}{r} Y_{1} = [Y_{2} | Z_{1}] β_{1} + ε_{1} = X_{1} β_{1} + ε_{1} \\ Y_{2} = [Y_{1} | Z_{2}] β_{2} + ε_{2} = X_{2} β_{2} + ε_{2} \end{array}

と書くことにする。外生変数からなる行列を $Z = [Z_{1} | Z_{2}]$ と表すと

MIMICモデル#

60年代にGoldbergerなど経済学者や社会科学者がパス解析を再発見した（『因果推論の科学』, p.159）

Goldberger & Karl Jöreskog (1975) はMIMIC modelを提案

参考文献#

パス解析

豊田秀樹（1998）『共分散構造分析入門編』
豊田秀樹（2014）『共分散構造分析 R編』

非逐次モデル

Williams, R (2015) Nonrecursive Models – Highlights
- 2SLSで解いている

パス解析

Contents

パス解析#

逐次モデル#

共分散の構造化#

行列を用いた共分散の構造化#

データを使って推定#

考察#

非逐次モデル#

考察：2SLS推定量とのつながり#

IV推定量#

MIMICモデル#

参考文献#