パス解析 - データサイエンス関連+αのメモ

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

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

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

逐次モデル¶

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

# import graphviz

# graphviz.Source("""digraph G {
#     newrank = true;
#     rankdir = RL;

#     node [shape = rectangle];
#     y1; y2; x1; x2;
#     node [shape = circle];
#     e1; e2;

#     subgraph cluster_0 {
#         color = white;
#         x1 -> y1
#         e1 -> y1
#     }
#     subgraph cluster_1 {
#         color = white;
#         x2 -> y2
#         e2 -> y2
#         { rank = max; e1 }
#     }

#     x2 -> y1
#     x1 -> y2

#     x1 -> x2 [tailport = w; headport = w; label=" "]
#     x2 -> x1 [tailport = w; headport = w; label=" "]
    
#     { rank = same; y1; y2; }
#     { rank = same; x1; x2; }
#     { rank = max; e1; e2; }
# }""")

import graphviz

graphviz.Source("""
digraph g {
    graph [rankdir = LR];
    node [shape = rectangle];
    v1 -> v3
    v2 -> v3 -> v4
    v2 -> v4

    v1 -> v2 [dir=both style=dashed label=" "]
    { rank = same; v1; v2 }
    { rank = same; v3; v4 }

    node [shape = circle]
    e3 -> v3
    e4 -> v4
    { rank = max; e3; e4 }
}
""")

Loading...

この図のモデルは

v_3 = a_{31} v_1 + a_{32} v_2 + e_3\\ v_4 = a_{42} v_2 + a_{43} v_3 + e_4

(1)

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

行列形式で

\newcommand{\b}[1]{\boldsymbol{#1}} \b{v} = \b{A v} + \b{e}

(2)

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

\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}

(3)

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

共分散の構造化¶

観測変数の共分散\b{\Sigma}は

Undefined control sequence: \b at position 1: \̲b̲{\Sigma}
= E[\b…

\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}

で、これらの各要素は

\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}

(5)

\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}

(6)

\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}

(7)

\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}

(8)

\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}

(9)

\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}

(10)

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

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

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

Undefined control sequence: \b at position 1: \̲b̲{v} = \b{Av} + …

\b{v} = \b{Av} + \b{e}

の式を

Undefined control sequence: \b at position 15: \begin{align}
\̲b̲{I v} &= \b{A v…

\begin{align}
\b{I v} &= \b{A v} + \b{e}\\
(\b{I} - \b{A})\b{v} &= \b{e}
\end{align}

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

Undefined control sequence: \b at position 1: \̲b̲{T} = (\b{I} - …

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

とおけば、構造方程式は

Undefined control sequence: \b at position 1: \̲b̲{v} = \b{T e}

\b{v} = \b{T e}

と簡潔に表記できる。

観測変数の共分散構造は

Undefined control sequence: \b at position 15: \begin{align}
\̲b̲{\Sigma}
&= E[\…

\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}

となる。

残差ベクトル\b{e}の共分散行列\b{\Sigma}_{\b{e}}は

Undefined control sequence: \b at position 1: \̲b̲{\Sigma_e}
=
\b…

\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}

となり、共分散構造は\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})として構造化した共分散を関数とする

Undefined control sequence: \b at position 1: \̲b̲{\Sigma}(\b{\th…

\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

負の対数尤度

Undefined control sequence: \b at position 20: …L} = \text{tr}(\̲b̲{\Sigma}(\b{\th…

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)

-1354.3911742902071

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[162], 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)

Loading...

# パス図
g = semopy.semplot(model, filename="/tmp/path_diagram.png", plot_covs=True)

with g.subgraph() as s:
    s.attr(rank = "same")
    s.node("v3")
    s.node("v4")

with g.subgraph() as s:
    s.attr(rank = "same")
    s.node("v1")
    s.node("v2")
    
g.attr(rankdir = "RL")
g.attr(pad = "0.2") # padding
g

Loading...

import graphviz

graphviz.Source("""
digraph g {
    graph [rankdir = RL];
    node [shape = rectangle];
        X1 [label=<<I>X<SUB>1</SUB></I>>]
        X2 [label=<<I>X<SUB>2</SUB></I>>]
        Y1 [label=<<I>Y<SUB>1</SUB></I>>]
        Y2 [label=<<I>Y<SUB>2</SUB></I>>]
    edge [];
        {X1, X2} -> Y1
        {X1, X2} -> Y2
        X1 -> X2 [dir=both style=dashed label=" "]
        {rank = same; X1; X2;}
}
""")

Loading...

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)

Loading...

# パス図
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

Loading...

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

import graphviz

graphviz.Source("""
digraph g {
    graph [rankdir = LR];
    node [shape = rectangle];
    v1 -> v2 -> v3
    v3 -> v2

    { rank = same; v2; v3 }

    node [shape = circle]
    e2 -> v2
    e3 -> v3
    { rank = max; e2; e3 }
}
""")

Loading...

\begin{align} v_2 &= a_{21} v_1 + a_{23} v_3 + e_2\\ v_3 &= a_{32} v_2 + e_3 \end{align}

(19)

\newcommand{\b}[1]{\boldsymbol{#1}} \b{v} = \b{A v} + \b{e}

(20)

\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}

(21)

Undefined control sequence: \b at position 1: \̲b̲{\Sigma}
= E[\b…

\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}

\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}

(23)

ゆえに

Undefined control sequence: \b at position 1: \̲b̲{\Sigma}
=
\beg…

\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}

\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}

(25)

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

IV推定量¶

同時方程式モデル

Y_1 = \beta_1 Y_2 + \gamma_1 Z_1 + \varepsilon_1\\ Y_2 = \beta_2 Y_1 + \gamma_2 Z_2 + \varepsilon_2

(26)

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

行列形式で $\boldsymbol{\beta}_1 = (\beta_1, \gamma_1)^T$ のようにまとめて

\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

(27)

と書くことにする。外生変数からなる行列を\b{Z} = [\b{Z}_1|\b{Z}_2]と表すと

$$

MIMICモデル¶

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

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

参考文献¶

パス解析

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

非逐次モデル

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