パス解析と構造方程式モデリング入門

ある企業の従業員の「満足度」は何によって決まるのでしょうか。「給与」が高ければ満足度も高い、という直接的な関係はありそうです。しかし、「給与が高い → 仕事のストレスが増える → 満足度が下がる」という間接的な経路もあるかもしれません。このような複数の経路を通じた因果的な影響を定量化するにはどうすればよいでしょうか。

あるいは教育の場面を考えてみましょう。学生の「最終成績」は、「授業への出席率」「自習時間」「モチベーション」といった要因が複雑に絡み合って決まります。出席率が高いと理解度が上がり、理解度が上がるとモチベーションも高まり、モチベーションが高まるとさらに自習時間が増える — こうした変数の連鎖的な影響を、重回帰分析だけで解きほぐすのは困難です。

パス解析(Path Analysis)は、変数間の因果関係を方向付きの経路(パス)として図示し、各経路の強さ(パス係数)を統計的に推定する手法です。1920年代に遺伝学者シューアル・ライトが提案したこの方法は、重回帰分析の自然な拡張であり、変数間の直接効果と間接効果を分離できます。

パス解析をさらに発展させ、直接観測できない潜在変数因子分析のモデル)を組み込んだものが構造方程式モデリング(SEM: Structural Equation Modeling)です。

パス解析とSEMを理解すると、以下のような応用が開けます。

  • 社会科学: 教育年数→職業→収入という社会的地位の達成モデル
  • 心理学: 態度→意図→行動という計画的行動理論のモデル
  • マーケティング: 品質→満足度→リピート購入の因果連鎖
  • 医学: リスク要因→中間変数→疾患の発症経路の分析
  • 教育学: 授業設計→理解度→成績の因果メカニズムの定量化

本記事では、パス解析の基本理論を定式化し、パス係数の推定方法を解説します。ライトのパス定理を厳密に導出し、効果の分解を数値例で確認します。SEMへの拡張の概要も示し、Pythonで実装して動作を確認します。

本記事の内容

  • パス解析の基本概念(パス図、直接効果と間接効果)
  • パス係数と標準化回帰係数の関係
  • 効果の分解(直接効果・間接効果・総合効果)
  • ライトのパス定理の導出と検証
  • 構造方程式モデリング(SEM)への拡張
  • モデルの識別と適合度評価
  • Pythonでの実装と可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

パス解析の基本概念

パス図の構成要素

パス解析の出発点は、研究者が仮説する因果関係をパス図(path diagram)として描くことです。パス図は有向非巡回グラフ(DAG: Directed Acyclic Graph)として表現され、以下の要素で構成されます。

  • 矩形(四角): 観測変数(直接測定できる変数)
  • 楕円: 潜在変数(SEMで使用、直接観測できない変数)
  • 片方向の矢印(→): 因果関係の方向。原因変数から結果変数へ向かう
  • 双方向の矢印(↔): 相関関係(因果の方向を仮定しない共変動)
  • パス係数: 矢印に付記される数値で、影響の強さと方向を表す

たとえば、$x_1 \to x_2 \to y$ というパス図は「$x_1$ が $x_2$ を通じて間接的に $y$ に影響する」というモデルを表します。$x_1 \to y$ という矢印も加えれば「$x_1$ は $y$ に直接的にも影響する」というモデルになります。

パス図を描くことの利点は、複雑な因果仮説を視覚的に明確にできることです。回帰分析だけでは「どの変数がどの変数を通じて影響しているか」という構造的な情報が失われますが、パス図はその構造を明示的に保持します。

外生変数と内生変数

パス図の変数は2種類に分けられます。

  • 外生変数(exogenous variables): モデル内の他の変数から矢印を受けない変数。因果連鎖の出発点であり、モデルの外部で決まる原因です
  • 内生変数(endogenous variables): モデル内の他の変数から矢印を受ける変数。少なくとも1つの原因変数を持ちます

各内生変数には攪乱変数(disturbance) $\epsilon$ が付随します。攪乱変数は、モデルに含めなかった変数の効果やランダムな変動を表し、内生変数の変動のうちモデルで説明できない部分に対応します。攪乱変数同士が相関する場合もあり、これは「モデルに含めていない共通の原因がある」ことを意味します。

たとえば、「教育年数($x_1$) → 職業的地位($x_2$) → 収入($y$)」というモデルでは、教育年数が外生変数、職業的地位と収入が内生変数です。職業的地位には「教育年数以外の要因」を表す攪乱変数 $\epsilon_2$ が付き、収入にも攪乱変数 $\epsilon_y$ が付きます。

変数の種類と構造が定まったところで、パス係数の正確な意味を見ていきましょう。

パス係数の定義と解釈

パス係数は、変数を標準化した(平均0、標準偏差1に変換した)上での標準化偏回帰係数です。標準化によって、元の変数のスケール(単位)に依存しない相対的な影響の強さが得られます。

変数 $y$ が変数 $x_1, x_2, \dots, x_k$ から矢印を受けているとき、$y$ の構造方程式は次のようになります。

$$ \begin{equation} y^* = p_{y \cdot x_1}x_1^* + p_{y \cdot x_2}x_2^* + \dots + p_{y \cdot x_k}x_k^* + p_{y \cdot \epsilon}\epsilon^* \end{equation} $$

ここで $*$ は標準化を表し、$p_{y \cdot x_j}$ がパス係数です。パス係数の解釈は「他の原因変数を一定に保ったとき、$x_j$ が1標準偏差増えると $y$ が $p_{y \cdot x_j}$ 標準偏差変化する」です。

パス係数と通常の(非標準化)回帰係数 $\beta_j$ の関係は次の通りです。

$$ \begin{equation} p_{y \cdot x_j} = \beta_j \cdot \frac{s_{x_j}}{s_y} \end{equation} $$

ここで $s_{x_j}$ と $s_y$ はそれぞれ $x_j$ と $y$ の標準偏差です。標準化回帰係数は変数のスケールを統一するため、異なるパスの影響力を直接比較できます。

パス係数と回帰分析の関係がわかったところで、具体的な推定方法に進みましょう。

パス係数の推定

回帰分析による推定

パス係数の推定は、各内生変数について、その変数を目的変数、直接の原因変数(矢印の発射元)を説明変数とする重回帰分析を行うことに帰着します。

パス図が非巡回(循環する因果経路がない)であれば、モデル全体のパス係数は各内生変数について独立に回帰分析を行うだけで推定できます。これを逐次的最小二乗法と呼びます。

具体例で考えてみましょう。次のモデルを仮定します。

$$ x_2 = p_{x_2 \cdot x_1} x_1 + \epsilon_2 $$

$$ y = p_{y \cdot x_1} x_1 + p_{y \cdot x_2} x_2 + \epsilon_y $$

推定は2段階で行います。

ステップ1: $x_2$ を $x_1$ で回帰し、パス係数 $\hat{p}_{x_2 \cdot x_1}$ を得る

ステップ2: $y$ を $x_1$ と $x_2$ で重回帰し、パス係数 $\hat{p}_{y \cdot x_1}$ と $\hat{p}_{y \cdot x_2}$ を得る

この手順は直感的にわかりやすいだけでなく、DAGの仮定のもとで一致推定量を与えることが知られています。

攪乱変数のパス係数

内生変数の攪乱変数のパス係数 $p_{y \cdot \epsilon}$ は、回帰のあてはまりから求められます。回帰の決定係数 $R^2$ を使うと次のようになります。

$$ \begin{equation} p_{y \cdot \epsilon} = \sqrt{1 – R^2} \end{equation} $$

$R^2$ が大きいほど攪乱変数のパス係数は小さくなり、モデルが変数 $y$ の変動をよく説明していることを意味します。$R^2 = 1$ であればモデルは完全に変動を説明し、攪乱変数は不要になります。

効果の分解という核心的なトピックに進みましょう。

効果の分解

直接効果・間接効果・総合効果

パス解析の最大の貢献は、変数間の総合効果直接効果間接効果に分解できることです。

$$ \begin{equation} \text{総合効果} = \text{直接効果} + \text{間接効果} \end{equation} $$

各効果の定義は次の通りです。

  • 直接効果(direct effect): $x$ から $y$ への直接のパス係数。他の変数を介さない影響
  • 間接効果(indirect effect): $x$ から1つ以上の中間変数を経由して $y$ に至るパスの係数の積の和。中間変数を介した影響
  • 総合効果(total effect): 直接効果と間接効果の和。$x$ が $y$ に及ぼす全体的な影響

間接効果の計算について、具体的に見ていきましょう。$x$ から $y$ に至る間接経路が複数ある場合、各経路の効果はその経路上のパス係数の積で計算されます。複数の間接経路がある場合は、それらの積をすべて合計します。

たとえば、$x \to m \to y$ と $x \to y$ の両方のパスがあるモデル(媒介モデル)では次のようになります。

$$ \text{直接効果} = p_{yx}, \quad \text{間接効果} = p_{mx} \times p_{ym}, \quad \text{総合効果} = p_{yx} + p_{mx} \times p_{ym} $$

複数の媒介変数がある場合

媒介変数が複数ある場合、間接効果はそれぞれの経路ごとに計算して合算します。たとえば、$x \to m_1 \to y$ と $x \to m_2 \to y$ の2つの間接経路がある場合、間接効果は次のようになります。

$$ \text{間接効果} = p_{m_1 x} \times p_{y m_1} + p_{m_2 x} \times p_{y m_2} $$

さらに、$x \to m_1 \to m_2 \to y$ のような逐次的な媒介経路がある場合、その効果は3つのパス係数の積 $p_{m_1 x} \times p_{m_2 m_1} \times p_{y m_2}$ となります。

このように、パス解析は複雑な因果構造を系統的に分解する枠組みを提供します。

数値例による効果の分解

具体的な数値で効果の分解を確認しましょう。以下のパスモデルを考えます。

$$ m = 0.6 \cdot x + \epsilon_m, \quad y = 0.3 \cdot x + 0.5 \cdot m + \epsilon_y $$

$x$ から $y$ への効果を分解すると次のようになります。

  • 直接効果: $p_{yx} = 0.3$
  • 間接効果: $p_{mx} \times p_{ym} = 0.6 \times 0.5 = 0.3$
  • 総合効果: $0.3 + 0.3 = 0.6$

この例では、直接効果と間接効果がちょうど同じ大きさです。つまり、$x$ が $y$ に与える影響のうち半分は直接的であり、残りの半分は媒介変数 $m$ を通じて間接的に伝わっています。このような分解は、介入の設計に直接的に役立ちます。もし $m$ の経路をブロックすれば、$x$ の $y$ への影響は半減するということがわかるからです。

効果の分解の妥当性は、ライトのパス定理によって理論的に裏付けられます。次にこの重要な定理を見ていきましょう。

ライトのパス定理

定理の内容

ライトは、変数間の相関係数をパス係数から計算する規則(パス定理)を導きました。この定理は、パス解析の理論的基盤を形成する最も重要な結果です。

ライトのパス定理: 2つの変数 $x_i$ と $x_j$ の間の相関係数は、それらを結ぶすべての合成パス(compound path)に沿ったパス係数の積の総和に等しい。

$$ \begin{equation} r_{x_i x_j} = \sum_{\text{全ての合成パス}}\prod_{\text{パスに沿った}}\text{パス係数} \end{equation} $$

合成パスの規則

合成パスとは、以下の規則に従って $x_i$ から $x_j$ へ辿ることができるパスです。

  1. パスの中で方向を変える(→から←へ、あるいは←から→へ)のは最大1回まで
  2. 同じ変数を2回以上通過しない
  3. 双方向矢印は方向変更とみなす

これらの規則から、合成パスは次の3種類に分類されます。

  • 直列パス: $x_i \to \cdots \to x_j$(すべて同じ方向)
  • 分岐パス: $x_i \leftarrow z \to x_j$(共通の原因 $z$ からの分岐。方向変更は $z$ で1回)
  • 合流パス: $x_i \to z \leftarrow x_j$(このパスは相関に寄与しない — $z$ が制御されていない場合)

定理の導出(簡略版)

3変数モデル $x \to m \to y$、$x \to y$ を例に、ライトのパス定理を導出します。構造方程式は次の通りです。

$$ m^* = p_{mx}x^* + p_{m\epsilon}\epsilon_m^*, \quad y^* = p_{yx}x^* + p_{ym}m^* + p_{y\epsilon}\epsilon_y^* $$

$x$ と $y$ の相関 $r_{xy}$ を直接計算します。$y^*$ に $m^*$ の式を代入すると次のようになります。

$$ y^* = p_{yx}x^* + p_{ym}(p_{mx}x^* + p_{m\epsilon}\epsilon_m^*) + p_{y\epsilon}\epsilon_y^* $$

$$ = (p_{yx} + p_{ym}p_{mx})x^* + p_{ym}p_{m\epsilon}\epsilon_m^* + p_{y\epsilon}\epsilon_y^* $$

$x^*$ と $\epsilon_m^*$ および $\epsilon_y^*$ は独立(相関ゼロ)なので、次のようになります。

$$ r_{xy} = E[x^* \cdot y^*] = p_{yx} + p_{ym}p_{mx} $$

これはまさに「$x$ から $y$ への直接パス $p_{yx}$ と、$x \to m \to y$ の間接パス $p_{ym}p_{mx}$ の和」であり、パス定理の結果と一致します。

パス定理によるモデルの検証

パス定理の重要な応用は、モデルの適合度の検証です。パス定理から計算される「モデルが含意する相関」と、データから計算される「標本相関」を比較することで、モデルの妥当性を評価できます。

過剰識別されたモデル(推定すべきパラメータの数より、データから得られる相関の数が多いモデル)では、モデルが含意する相関と標本相関が一致しない可能性があります。両者の乖離が大きければ、仮定した因果構造が不適切であることを示唆します。

SEMへの拡張を見ていきましょう。

構造方程式モデリング(SEM)への拡張

パス解析の限界

パス解析は観測変数のみを扱います。しかし、「知能」「満足度」「動機づけ」「社会的地位」といった概念は直接観測できず、テストスコアやアンケート回答などの観測指標を通じて間接的に測定されます。

これらの直接観測できない概念(潜在変数、latent variables)をモデルに組み込むためにパス解析を拡張したのがSEMです。潜在変数を導入する利点は主に2つあります。

  1. 測定誤差の考慮: 観測指標には必ず測定誤差が含まれます。潜在変数を使うことで、真のスコアと測定誤差を分離し、測定誤差によるバイアス(減衰)を補正できます
  2. 構成概念の表現: 複数の指標で1つの概念を測定することで、より信頼性の高い測定が実現します

SEMの構成

SEMは2つのサブモデルで構成されます。

測定モデル(measurement model): 潜在変数と観測指標の関係を記述します。これは因子分析のモデルそのものです。

$$ \begin{equation} \bm{x} = \bm{\Lambda}_x\bm{\xi} + \bm{\delta} \end{equation} $$

$$ \begin{equation} \bm{y} = \bm{\Lambda}_y\bm{\eta} + \bm{\varepsilon} \end{equation} $$

ここで $\bm{\Lambda}_x$ と $\bm{\Lambda}_y$ は因子負荷行列、$\bm{\xi}$ は外生潜在変数、$\bm{\eta}$ は内生潜在変数、$\bm{\delta}$ と $\bm{\varepsilon}$ は測定誤差です。

構造モデル(structural model): 潜在変数間の因果関係を記述します。これはパス解析のモデルです。

$$ \begin{equation} \bm{\eta} = \bm{B}\bm{\eta} + \bm{\Gamma}\bm{\xi} + \bm{\zeta} \end{equation} $$

ここで $\bm{B}$ は内生変数間のパス係数行列(対角成分はゼロ)、$\bm{\Gamma}$ は外生変数から内生変数へのパス係数行列、$\bm{\zeta}$ は構造攪乱(disturbance)です。

$\bm{B}$ の非対角成分 $B_{ij}$ は「潜在変数 $\eta_j$ が $\eta_i$ に与える直接効果」を表します。

モデルが含意する共分散行列

SEMの分析において中核となるのは、モデルが含意する共分散行列 $\bm{\Sigma}(\bm{\theta})$ と、データから計算される標本共分散行列 $\bm{S}$ の比較です。

構造モデルから $\bm{\eta}$ を解くと次のようになります。

$$ \bm{\eta} = (\bm{I} – \bm{B})^{-1}(\bm{\Gamma}\bm{\xi} + \bm{\zeta}) $$

これを使って、観測変数の共分散行列をモデルパラメータ $\bm{\theta}$ の関数として表現できます。SEMの推定は、$\bm{\Sigma}(\bm{\theta})$ と $\bm{S}$ の乖離を最小化するパラメータ $\hat{\bm{\theta}}$ を求める問題です。

モデルの識別

SEMでは、モデルが識別されている(identified)かどうかを確認する必要があります。識別とは、データ(共分散行列の要素)からモデルの全パラメータが一意に決定できることを意味します。

識別の必要条件は次の通りです。

$$ \begin{equation} \text{推定すべきパラメータ数} \leq \frac{p(p+1)}{2} \end{equation} $$

ここで $p$ は観測変数の数、$p(p+1)/2$ は共分散行列の非重複要素の数です。右辺はデータから得られる独立な情報量を表します。

モデルの適合度指標

SEMではモデルの適合度を評価するために以下の指標が使われます。

  • カイ二乗統計量: $\chi^2 = (n-1)F_{\min}$。モデルが含意する共分散行列と標本共分散行列の乖離を測る。小さいほど適合がよい。ただしサンプルサイズに影響されるため、他の指標と併用することが推奨される
  • RMSEA(Root Mean Square Error of Approximation): 近似の誤差を自由度で調整した指標。0.05以下が良好、0.08以下が許容範囲とされる
  • CFI(Comparative Fit Index): 独立モデル(変数間に関係がないモデル)と比較した適合度。0.95以上が良好
  • TLI(Tucker-Lewis Index): CFIと同様だが自由度にペナルティを課す。0.95以上が良好
  • SRMR(Standardized Root Mean Square Residual): 標準化残差の平均。0.08以下が良好

複数の指標を総合的に判断することが推奨されます。1つの指標だけでモデルの適合度を判定するのは危険です。

Pythonで実装して確認しましょう。

Pythonによる実装と可視化

パス解析の実装

3つの変数のパス解析を実装し、効果の分解を行います。$x \to m \to y$ と $x \to y$ の両方のパスを持つ媒介モデルです。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
n = 500

# --- 真のパスモデル: x -> m -> y, x -> y ---
# x: 外生変数
x = np.random.randn(n)

# m: 中間変数 (x -> m)
p_mx = 0.6  # パス係数
m = p_mx * x + np.sqrt(1 - p_mx**2) * np.random.randn(n)

# y: 結果変数 (x -> y, m -> y)
p_yx = 0.3   # 直接効果
p_ym = 0.5   # m -> y のパス係数
noise_var = 1 - p_yx**2 - p_ym**2 - 2*p_yx*p_ym*p_mx
y = p_yx * x + p_ym * m + np.sqrt(max(noise_var, 0.01)) * np.random.randn(n)

# 標準化
x = (x - x.mean()) / x.std()
m = (m - m.mean()) / m.std()
y = (y - y.mean()) / y.std()

# --- パス係数の推定(回帰分析) ---
# Step 1: m を x で回帰
X1 = np.column_stack([np.ones(n), x])
beta_mx = np.linalg.lstsq(X1, m, rcond=None)[0]
p_mx_hat = beta_mx[1]

# Step 2: y を x と m で回帰
X2 = np.column_stack([np.ones(n), x, m])
beta_y = np.linalg.lstsq(X2, y, rcond=None)[0]
p_yx_hat = beta_y[1]  # 直接効果
p_ym_hat = beta_y[2]  # m -> y

# --- 効果の分解 ---
direct_effect = p_yx_hat
indirect_effect = p_mx_hat * p_ym_hat
total_effect = direct_effect + indirect_effect

# --- 攪乱変数のパス係数 ---
R2_m = 1 - np.sum((m - X1 @ beta_mx)**2) / np.sum((m - m.mean())**2)
R2_y = 1 - np.sum((y - X2 @ beta_y)**2) / np.sum((y - y.mean())**2)
p_epsilon_m = np.sqrt(1 - R2_m)
p_epsilon_y = np.sqrt(1 - R2_y)

print("=== パス解析の結果 ===")
print(f"\nパス係数:")
print(f"  x -> m: {p_mx_hat:.4f} (真値: {p_mx})")
print(f"  x -> y (直接): {p_yx_hat:.4f} (真値: {p_yx})")
print(f"  m -> y: {p_ym_hat:.4f} (真値: {p_ym})")
print(f"\n攪乱変数のパス係数:")
print(f"  ε_m -> m: {p_epsilon_m:.4f}")
print(f"  ε_y -> y: {p_epsilon_y:.4f}")
print(f"\n効果の分解:")
print(f"  直接効果 (x -> y): {direct_effect:.4f}")
print(f"  間接効果 (x -> m -> y): {indirect_effect:.4f}")
print(f"  総合効果: {total_effect:.4f}")
print(f"\n相関 r(x, y) = {np.corrcoef(x, y)[0,1]:.4f}")
print(f"総合効果と相関の差: {abs(total_effect - np.corrcoef(x, y)[0,1]):.4f}")

この結果から、パス解析の基本的な性質が確認できます。各パス係数の推定値は真の値に近く、$n = 500$ のサンプルサイズで良好な推定が得られています。攪乱変数のパス係数は各内生変数の説明されない変動の割合を表しており、$R^2$ と対応関係にあります。

注目すべきは、総合効果と $x$-$y$ 間の相関がほぼ一致している点です。これはライトのパス定理の数値的な確認です。外生変数から内生変数への総合効果は、両者の相関と等しくなるという理論が、データにより裏付けられています。

パス図と効果の分解の可視化

パス図を描画し、効果の分解を棒グラフで可視化します。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
n = 500
x = np.random.randn(n)
p_mx, p_yx, p_ym = 0.6, 0.3, 0.5
m = p_mx * x + np.sqrt(1 - p_mx**2) * np.random.randn(n)
noise_var = 1 - p_yx**2 - p_ym**2 - 2*p_yx*p_ym*p_mx
y = p_yx * x + p_ym * m + np.sqrt(max(noise_var, 0.01)) * np.random.randn(n)
x = (x - x.mean()) / x.std()
m = (m - m.mean()) / m.std()
y = (y - y.mean()) / y.std()

X1 = np.column_stack([np.ones(n), x])
beta_mx = np.linalg.lstsq(X1, m, rcond=None)[0]
p_mx_hat = beta_mx[1]
X2 = np.column_stack([np.ones(n), x, m])
beta_y = np.linalg.lstsq(X2, y, rcond=None)[0]
p_yx_hat = beta_y[1]
p_ym_hat = beta_y[2]
direct_effect = p_yx_hat
indirect_effect = p_mx_hat * p_ym_hat
total_effect = direct_effect + indirect_effect

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# (a) パス図
ax = axes[0]
ax.set_xlim(-0.5, 4.5)
ax.set_ylim(-1, 3)
for pos, label in [((0, 1), "x"), ((2, 2.2), "m"), ((4, 1), "y")]:
    ax.add_patch(plt.Rectangle((pos[0]-0.4, pos[1]-0.3), 0.8, 0.6,
                                fill=True, facecolor="lightblue",
                                edgecolor="black", linewidth=1.5))
    ax.text(pos[0], pos[1], label, fontsize=14, ha="center", va="center",
            fontweight="bold")
ax.annotate("", xy=(1.6, 2.2), xytext=(0.4, 1.3),
            arrowprops=dict(arrowstyle="->", lw=2, color="steelblue"))
ax.text(0.7, 2.0, f"{p_mx_hat:.3f}", fontsize=11, color="steelblue",
        fontweight="bold")
ax.annotate("", xy=(3.6, 1.3), xytext=(2.4, 2.0),
            arrowprops=dict(arrowstyle="->", lw=2, color="steelblue"))
ax.text(3.3, 2.0, f"{p_ym_hat:.3f}", fontsize=11, color="steelblue",
        fontweight="bold")
ax.annotate("", xy=(3.6, 1), xytext=(0.4, 1),
            arrowprops=dict(arrowstyle="->", lw=2, color="coral"))
ax.text(2, 0.6, f"{p_yx_hat:.3f}", fontsize=11, color="coral",
        fontweight="bold")
ax.set_title("(a) Path diagram", fontsize=12)
ax.axis("off")

# (b) 効果の分解
ax = axes[1]
effects = [direct_effect, indirect_effect, total_effect]
labels = ["Direct\n(x->y)", "Indirect\n(x->m->y)", "Total"]
colors = ["coral", "steelblue", "gray"]
bars = ax.bar(labels, effects, color=colors, alpha=0.7, edgecolor="black")
for bar, val in zip(bars, effects):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.01,
            f"{val:.3f}", ha="center", fontsize=11)
ax.set_ylabel("Effect size", fontsize=11)
ax.set_title("(b) Decomposition of effects", fontsize=12)
ax.grid(True, alpha=0.3, axis="y")

plt.tight_layout()
plt.savefig("path_analysis.png", dpi=150, bbox_inches="tight")
plt.show()

この可視化から、パス解析の力が確認できます。

  1. パス図(図a): 推定されたパス係数が図上に表示されています。直接パス(赤、$x \to y$)と間接パスの各要素(青、$x \to m$ と $m \to y$)の強さが視覚的に比較できます。間接経路のうち $x \to m$ が0.6程度と最も強く、$m$ を通じた影響が大きいことが一目でわかります

  2. 効果の分解(図b): $x$ から $y$ への総合効果が直接効果と間接効果に分解されています。間接効果($x \to m \to y$)が直接効果とほぼ同程度の大きさを持っており、中間変数 $m$ を通じた経路が重要であることがわかります。もし $m$ への介入(たとえば $m$ を固定する)を行うと、間接効果が遮断され、$x$ の $y$ への影響は約半分に減少します

ライトのパス定理の検証

ライトのパス定理が成り立つことを、すべての変数ペアについて検証します。

import numpy as np

np.random.seed(42)
n = 500
x = np.random.randn(n)
p_mx, p_yx, p_ym = 0.6, 0.3, 0.5
m = p_mx * x + np.sqrt(1 - p_mx**2) * np.random.randn(n)
noise_var = 1 - p_yx**2 - p_ym**2 - 2*p_yx*p_ym*p_mx
y = p_yx * x + p_ym * m + np.sqrt(max(noise_var, 0.01)) * np.random.randn(n)
x = (x - x.mean()) / x.std()
m = (m - m.mean()) / m.std()
y = (y - y.mean()) / y.std()

# 推定
X1 = np.column_stack([np.ones(n), x])
p_mx_hat = np.linalg.lstsq(X1, m, rcond=None)[0][1]
X2 = np.column_stack([np.ones(n), x, m])
betas = np.linalg.lstsq(X2, y, rcond=None)[0]
p_yx_hat, p_ym_hat = betas[1], betas[2]

# 標本相関
r_xm = np.corrcoef(x, m)[0, 1]
r_xy = np.corrcoef(x, y)[0, 1]
r_my = np.corrcoef(m, y)[0, 1]

# パス定理から計算される相関
# r(x, m) の合成パス: x -> m だけなので r = p_mx
r_xm_model = p_mx_hat

# r(x, y) の合成パス: x -> y と x -> m -> y
r_xy_model = p_yx_hat + p_mx_hat * p_ym_hat

# r(m, y) の合成パス: m -> y と m <- x -> y
r_my_model = p_ym_hat + p_mx_hat * p_yx_hat

print("=== ライトのパス定理の検証 ===")
print(f"{'ペア':>10s} {'標本相関':>10s} {'モデル含意':>10s} {'差':>10s}")
print(f"{'r(x,m)':>10s} {r_xm:10.4f} {r_xm_model:10.4f} {abs(r_xm - r_xm_model):10.4f}")
print(f"{'r(x,y)':>10s} {r_xy:10.4f} {r_xy_model:10.4f} {abs(r_xy - r_xy_model):10.4f}")
print(f"{'r(m,y)':>10s} {r_my:10.4f} {r_my_model:10.4f} {abs(r_my - r_my_model):10.4f}")

この検証結果から、ライトのパス定理の正確さが確認できます。全ての変数ペアについて、パス定理から計算されるモデル含意の相関と標本相関がよく一致しています。特に注目すべきは $r(m, y)$ の計算です。$m$ から $y$ への直接パス $p_{ym}$ だけでなく、$m \leftarrow x \to y$ という分岐パス(共通原因 $x$ を通じた間接的な関連)も含まれています。この分岐パスの効果 $p_{mx} \times p_{yx}$ を合算しないと、標本相関と一致しません。

4変数モデルの拡張

より複雑な4変数のパスモデルを実装し、複数の間接経路がある場合の効果分解を行います。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
n = 1000

# --- 4変数モデル: x -> m1 -> y, x -> m2 -> y, x -> y ---
x = np.random.randn(n)

# m1, m2: 2つの媒介変数
p_m1x = 0.5
p_m2x = 0.4
m1 = p_m1x * x + np.sqrt(1 - p_m1x**2) * np.random.randn(n)
m2 = p_m2x * x + np.sqrt(1 - p_m2x**2) * np.random.randn(n)

# y: 結果変数
p_yx = 0.2
p_ym1 = 0.3
p_ym2 = 0.25
var_noise = 1 - p_yx**2 - p_ym1**2 - p_ym2**2 - 2*p_yx*p_m1x*p_ym1 - 2*p_yx*p_m2x*p_ym2 - 2*p_ym1*p_ym2*p_m1x*p_m2x
y = p_yx*x + p_ym1*m1 + p_ym2*m2 + np.sqrt(max(var_noise, 0.01))*np.random.randn(n)

# 標準化
for arr in [x, m1, m2, y]:
    arr -= arr.mean()
    arr /= arr.std()

# 推定
ones = np.ones(n)
p_m1x_hat = np.linalg.lstsq(np.column_stack([ones, x]), m1, rcond=None)[0][1]
p_m2x_hat = np.linalg.lstsq(np.column_stack([ones, x]), m2, rcond=None)[0][1]
betas_y = np.linalg.lstsq(np.column_stack([ones, x, m1, m2]), y, rcond=None)[0]
p_yx_hat, p_ym1_hat, p_ym2_hat = betas_y[1], betas_y[2], betas_y[3]

# 効果の分解
direct = p_yx_hat
indirect_via_m1 = p_m1x_hat * p_ym1_hat
indirect_via_m2 = p_m2x_hat * p_ym2_hat
total_indirect = indirect_via_m1 + indirect_via_m2
total = direct + total_indirect

print("=== 4変数パスモデルの結果 ===")
print(f"パス係数:")
print(f"  x -> m1: {p_m1x_hat:.4f} (真値: {p_m1x})")
print(f"  x -> m2: {p_m2x_hat:.4f} (真値: {p_m2x})")
print(f"  x -> y:  {p_yx_hat:.4f} (真値: {p_yx})")
print(f"  m1 -> y: {p_ym1_hat:.4f} (真値: {p_ym1})")
print(f"  m2 -> y: {p_ym2_hat:.4f} (真値: {p_ym2})")
print(f"\n効果の分解:")
print(f"  直接効果 (x -> y):       {direct:.4f}")
print(f"  間接効果 (x -> m1 -> y): {indirect_via_m1:.4f}")
print(f"  間接効果 (x -> m2 -> y): {indirect_via_m2:.4f}")
print(f"  間接効果合計:            {total_indirect:.4f}")
print(f"  総合効果:                {total:.4f}")
print(f"  相関 r(x,y):             {np.corrcoef(x, y)[0,1]:.4f}")

# 可視化
fig, ax = plt.subplots(figsize=(8, 5))
effects_list = [direct, indirect_via_m1, indirect_via_m2, total_indirect, total]
labels_list = ["Direct\n(x->y)", "Indirect\n(x->m1->y)", "Indirect\n(x->m2->y)",
               "Total\nIndirect", "Total"]
colors_list = ["coral", "steelblue", "mediumpurple", "teal", "gray"]
bars = ax.bar(labels_list, effects_list, color=colors_list, alpha=0.7, edgecolor="black")
for bar, val in zip(bars, effects_list):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.005,
            f"{val:.3f}", ha="center", fontsize=10)
ax.set_ylabel("Effect size", fontsize=11)
ax.set_title("Decomposition of effects (4-variable model)", fontsize=12)
ax.grid(True, alpha=0.3, axis="y")
plt.tight_layout()
plt.savefig("path_analysis_4var.png", dpi=150, bbox_inches="tight")
plt.show()

この4変数モデルの結果から、以下が読み取れます。

  1. 複数の間接経路: $x$ から $y$ への間接効果が $m_1$ 経由と $m_2$ 経由の2つに分解されています。$m_1$ 経由の間接効果の方がやや大きく、これは $p_{m_1 x}$ と $p_{y m_1}$ が共に大きいことに起因します

  2. 直接効果の相対的な大きさ: 直接効果は間接効果合計よりも小さくなっています。つまり、$x$ が $y$ に与える影響の大部分は媒介変数を通じた間接的な経路によるものです

  3. 介入の設計への示唆: もし $m_1$ の経路だけを遮断すると、$x$ の $y$ への影響は約0.15減少します。$m_2$ の経路も遮断すれば、直接効果の約0.2だけが残ります。このように、効果分解の結果は具体的な介入戦略の設計に直結します

サンプルサイズと推定精度の検証

パス係数の推定精度がサンプルサイズにどう依存するかを、モンテカルロシミュレーションで検証します。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
p_mx, p_yx, p_ym = 0.6, 0.3, 0.5
sample_sizes = [50, 100, 200, 500, 1000, 2000]
n_sim = 500

results = {ss: {"direct": [], "indirect": [], "total": []} for ss in sample_sizes}

for ss in sample_sizes:
    for _ in range(n_sim):
        x = np.random.randn(ss)
        m_sim = p_mx * x + np.sqrt(1 - p_mx**2) * np.random.randn(ss)
        nv = 1 - p_yx**2 - p_ym**2 - 2*p_yx*p_ym*p_mx
        y_sim = p_yx * x + p_ym * m_sim + np.sqrt(max(nv, 0.01)) * np.random.randn(ss)
        x_s = (x - x.mean()) / x.std()
        m_s = (m_sim - m_sim.mean()) / m_sim.std()
        y_s = (y_sim - y_sim.mean()) / y_sim.std()

        ones = np.ones(ss)
        pmx = np.linalg.lstsq(np.column_stack([ones, x_s]), m_s, rcond=None)[0][1]
        b = np.linalg.lstsq(np.column_stack([ones, x_s, m_s]), y_s, rcond=None)[0]
        results[ss]["direct"].append(b[1])
        results[ss]["indirect"].append(pmx * b[2])
        results[ss]["total"].append(b[1] + pmx * b[2])

fig, axes = plt.subplots(1, 3, figsize=(16, 5))
true_vals = {"direct": p_yx, "indirect": p_mx * p_ym, "total": p_yx + p_mx * p_ym}
titles = {"direct": "Direct effect", "indirect": "Indirect effect", "total": "Total effect"}

for idx, key in enumerate(["direct", "indirect", "total"]):
    ax = axes[idx]
    data = [results[ss][key] for ss in sample_sizes]
    bp = ax.boxplot(data, labels=[str(ss) for ss in sample_sizes],
                    patch_artist=True)
    for patch in bp["boxes"]:
        patch.set_facecolor("lightblue")
    ax.axhline(true_vals[key], color="red", linestyle="--", linewidth=1.5,
               label=f"True = {true_vals[key]:.2f}")
    ax.set_xlabel("Sample size $n$", fontsize=11)
    ax.set_ylabel("Estimated effect", fontsize=11)
    ax.set_title(f"({chr(97+idx)}) {titles[key]}", fontsize=12)
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("path_analysis_mc.png", dpi=150, bbox_inches="tight")
plt.show()

モンテカルロシミュレーションの結果から、パス解析の推定特性が明らかになります。

  1. 一致性: すべてのサンプルサイズで推定値の中央値は真の値の近くにあり、サンプルサイズが大きくなるほど箱ひげ図の箱が狭くなっています。これは推定量の一致性(サンプルサイズ増加とともに真の値に収束する性質)を視覚的に確認しています

  2. 間接効果の精度: 間接効果(図b)は直接効果(図a)と比較して、同じサンプルサイズでのばらつきがやや大きい傾向があります。これは間接効果が2つのパス係数の積であるため、推定誤差が乗法的に蓄積されるためです

  3. 必要サンプルサイズ: $n = 200$ 以上であれば、効果の推定値の分布がかなり安定しています。$n = 50$ ではばらつきが大きく、偶然による影響を受けやすい状態です。パス解析では一般に $n \geq 200$ が推奨されています

まとめ

本記事では、パス解析と構造方程式モデリングの基本理論を解説しました。

  • パス解析: 変数間の因果関係をパス図で表現し、パス係数を回帰分析で推定する手法。1920年代にシューアル・ライトが提案しました
  • パス係数: 標準化偏回帰係数であり、「他の原因を一定にしたとき、原因変数が1標準偏差変わると結果変数が何標準偏差変わるか」を表します
  • 効果の分解: 総合効果 = 直接効果 + 間接効果。中間変数を通じた間接的な影響を定量化でき、介入の設計に直結します
  • ライトのパス定理: 変数間の相関はすべての合成パスに沿ったパス係数の積の総和に等しくなります。モデルの適合度検証に利用できます
  • SEMへの拡張: 潜在変数を組み込むことで、直接観測できない概念間の因果関係をモデル化できます。測定モデルと構造モデルの2つのサブモデルで構成されます
  • 注意: パス解析は仮定した因果構造のもとでパス係数を推定するものであり、因果関係そのものを証明するものではありません。因果的な解釈は、理論的な根拠や研究デザイン(実験、時間的順序など)に依拠する必要があります

次のステップとして、以下の記事も参考にしてください。