「ある薬を飲んだら病気が治った」「広告を出したら売上が上がった」——こうした経験から因果関係を結論づけてしまいがちですが、本当にその介入が原因なのでしょうか。単なる相関と真の因果を区別することは、統計学における最も重要な問題の一つです。
因果推論（causal inference） は、観察データやランダム化実験のデータから、ある介入（処置）が結果に与える因果効果を推定するための統計的枠組みです。本記事では、Donald Rubinが体系化した 潜在結果フレームワーク（potential outcomes framework） を中心に、因果効果の定義、識別条件、推定手法を数学的に導出します。
本記事の内容
- 因果推論の根本問題（反事実の観測不可能性）
- ルービンの因果モデル（潜在結果フレームワーク）
- ATE・ATT・ATUの定義と関係
- SUTVA仮定
- 無視可能性条件とランダム化実験
- 傾向スコアの定義とIPW推定量の導出
- 傾向スコアマッチング
- Pythonでの合成データを用いたATE推定
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
因果推論の根本問題
相関と因果の違い
統計学の基本原則として、相関は因果を意味しない（correlation does not imply causation） というものがあります。例えば、「アイスクリームの売上が増えると溺水事故が増える」というデータがあっても、アイスクリームが溺水を引き起こしているわけではありません。気温という交絡因子（confounding variable）が両方に影響しています。
反事実と因果の定義
因果効果を厳密に定義するには、反事実（counterfactual） という概念が必要です。ある個人が薬を飲んだとき病気が治ったとします。因果効果を知るには、「同じ個人が薬を飲まなかったとしたら、どうなっていたか」という反事実の情報が必要です。
しかし、同じ個人が同時に「薬を飲む」状態と「薬を飲まない」状態の両方を経験することはできません。これが 因果推論の根本問題（fundamental problem of causal inference） です。
ルービンの因果モデル
潜在結果フレームワーク
$n$ 人の個人（単位）$i = 1, 2, \dots, n$ を考えます。各個人に対して二値の処置変数 $W_i$ を定義します。
$$ W_i = \begin{cases} 1 & \text{（処置群：処置を受けた場合）} \\ 0 & \text{（対照群：処置を受けなかった場合）} \end{cases} $$
潜在結果
各個人 $i$ に対して、2つの 潜在結果（potential outcomes） を定義します。
$$ Y_i(1) : \text{個人 } i \text{ が処置を受けた場合の結果} $$
$$ Y_i(0) : \text{個人 } i \text{ が処置を受けなかった場合の結果} $$
$Y_i(1)$ と $Y_i(0)$ は処置の割り当てに関わらず存在する値であり、処置前から決まっていると考えます。
個人因果効果
個人 $i$ に対する 個人因果効果（Individual Causal Effect, ICE） は、
$$ \tau_i = Y_i(1) – Y_i(0) $$
と定義されます。しかし、根本問題により、観測できるのは
$$ Y_i^{\text{obs}} = W_i \cdot Y_i(1) + (1 – W_i) \cdot Y_i(0) $$
のみです。つまり、$W_i = 1$ のとき $Y_i(1)$ のみ、$W_i = 0$ のとき $Y_i(0)$ のみが観測されます。$\tau_i$ を直接計算することは不可能です。
観測データの構造
個人 $i$ について観測されるデータは三つ組 $(Y_i^{\text{obs}}, W_i, \bm{X}_i)$ です。
- $Y_i^{\text{obs}}$：観測された結果変数
- $W_i$：処置の割り当て
- $\bm{X}_i$：共変量（交絡因子を含む背景変数）
因果効果の定義
ATE（Average Treatment Effect）
個人因果効果の期待値（母集団全体の平均因果効果）を ATE と呼びます。
$$ \begin{equation} \tau_{\text{ATE}} = E[Y_i(1) – Y_i(0)] = E[Y(1)] – E[Y(0)] \end{equation} $$
ATT（Average Treatment Effect on the Treated）
処置を実際に受けた集団における平均因果効果を ATT と呼びます。
$$ \begin{equation} \tau_{\text{ATT}} = E[Y(1) – Y(0) \mid W = 1] = E[Y(1) \mid W = 1] – E[Y(0) \mid W = 1] \end{equation} $$
第2項 $E[Y(0) \mid W = 1]$ は「処置を受けた人が、もし処置を受けなかったとしたらどうなっていたか」という反事実の期待値であり、直接観測できません。
ATU（Average Treatment Effect on the Untreated）
処置を受けなかった集団における平均因果効果を ATU と呼びます。
$$ \tau_{\text{ATU}} = E[Y(1) – Y(0) \mid W = 0] $$
3つの因果効果の関係
処置群の割合を $\pi = P(W = 1)$ とすると、
$$ \begin{align} \tau_{\text{ATE}} &= \pi \cdot \tau_{\text{ATT}} + (1 – \pi) \cdot \tau_{\text{ATU}} \end{align} $$
導出：
$$ \begin{align} \tau_{\text{ATE}} &= E[Y(1) – Y(0)] \\ &= E[Y(1) – Y(0) \mid W = 1]P(W = 1) + E[Y(1) – Y(0) \mid W = 0]P(W = 0) \\ &= \pi \cdot \tau_{\text{ATT}} + (1 – \pi) \cdot \tau_{\text{ATU}} \end{align} $$
全確率の法則を用いました。
SUTVA仮定
定義
SUTVA（Stable Unit Treatment Value Assumption） は以下の2つの条件からなります。
- 干渉なし（no interference）：個人 $i$ の潜在結果は他の個人の処置割り当てに影響されない
$$ Y_i(W_1, \dots, W_n) = Y_i(W_i) $$
- 処置の一貫性（consistency）：処置のバージョンは1つのみ。同じ $W_i = 1$ でも処置の内容が人によって異なることはない
SUTVAが成り立たない例としては、感染症のワクチン接種があります。ある人がワクチンを接種すると、その人の周囲の人の感染リスクも下がるため、干渉なしの仮定が破れます。
ナイーブな比較のバイアス
観測データから素朴（ナイーブ）に処置効果を推定すると、
$$ \hat{\tau}_{\text{naive}} = E[Y^{\text{obs}} \mid W = 1] – E[Y^{\text{obs}} \mid W = 0] = E[Y(1) \mid W = 1] – E[Y(0) \mid W = 0] $$
しかし、ATEとの差（バイアス）は、
$$ \begin{align} \hat{\tau}_{\text{naive}} – \tau_{\text{ATE}} &= E[Y(1) \mid W = 1] – E[Y(0) \mid W = 0] – E[Y(1)] + E[Y(0)] \end{align} $$
$E[Y(1)] = \pi E[Y(1)|W=1] + (1-\pi)E[Y(1)|W=0]$ などを用いて整理すると、
$$ \begin{align} \hat{\tau}_{\text{naive}} – \tau_{\text{ATE}} &= (1-\pi)\{E[Y(1)|W=1] – E[Y(1)|W=0]\} \\ &\quad + \pi\{E[Y(0)|W=1] – E[Y(0)|W=0]\} \end{align} $$
右辺は 選択バイアス（selection bias） を表しています。処置群と対照群で潜在結果の分布が異なる場合（つまり処置の割り当てが潜在結果と独立でない場合）、ナイーブな比較にはバイアスが生じます。
無視可能性条件とランダム化実験
強い無視可能性条件
条件1：条件付き独立性（unconfoundedness / ignorability）
$$ \begin{equation} (Y(1), Y(0)) \perp\!\!\!\perp W \mid \bm{X} \end{equation} $$
共変量 $\bm{X}$ で条件づけると、処置の割り当て $W$ は潜在結果と独立になります。
条件2：正値性（positivity / overlap）
$$ 0 < P(W = 1 \mid \bm{X} = \bm{x}) < 1 \quad \text{for all } \bm{x} $$
どの共変量の値でも、処置を受ける確率と受けない確率がともに正であることを意味します。
これら2つの条件を合わせて 強い無視可能性（strong ignorability） と呼びます。
無視可能性の下でのATE識別
条件付き独立性の下では、
$$ E[Y(1) \mid \bm{X}] = E[Y(1) \mid W = 1, \bm{X}] = E[Y^{\text{obs}} \mid W = 1, \bm{X}] $$
第1の等号は条件付き独立性、第2の等号は一貫性の仮定によります。同様に、
$$ E[Y(0) \mid \bm{X}] = E[Y^{\text{obs}} \mid W = 0, \bm{X}] $$
したがって、条件付きATE（CATE）は、
$$ \tau(\bm{x}) = E[Y^{\text{obs}} \mid W = 1, \bm{X} = \bm{x}] – E[Y^{\text{obs}} \mid W = 0, \bm{X} = \bm{x}] $$
これを $\bm{X}$ の分布で周辺化すると ATE が得られます。
$$ \tau_{\text{ATE}} = E_{\bm{X}}[\tau(\bm{X})] $$
ランダム化実験がなぜ機能するか
ランダム化比較試験（RCT）では、処置の割り当て $W$ を乱数で決定します。これにより、
$$ W \perp\!\!\!\perp (Y(1), Y(0), \bm{X}) $$
が保証されます。無視可能性条件が共変量 $\bm{X}$ で条件づけなくても成立するため、
$$ \tau_{\text{ATE}} = E[Y^{\text{obs}} \mid W = 1] – E[Y^{\text{obs}} \mid W = 0] $$
をナイーブに計算してもバイアスが生じません。これがRCTが因果推論のゴールドスタンダードである理由です。
傾向スコア
定義
傾向スコア（propensity score） は、共変量 $\bm{X}$ が与えられたときに処置を受ける条件付き確率です。
$$ \begin{equation} e(\bm{x}) = P(W = 1 \mid \bm{X} = \bm{x}) \end{equation} $$
傾向スコアの性質（Rosenbaum-Rubin, 1983）
定理： 強い無視可能性が成り立つとき、
$$ (Y(1), Y(0)) \perp\!\!\!\perp W \mid e(\bm{X}) $$
すなわち、潜在結果と処置割り当ての独立性は、多次元の $\bm{X}$ で条件づける代わりに、1次元の傾向スコア $e(\bm{X})$ で条件づけるだけで十分です。
証明のスケッチ： 任意の $\bm{X}$ の関数 $h(\bm{X})$ に対して、
$$ \begin{align} E[h(\bm{X}) \cdot W \mid e(\bm{X})] &= E[E[h(\bm{X}) \cdot W \mid \bm{X}] \mid e(\bm{X})] \quad (\text{繰り返し期待値の法則}) \\ &= E[h(\bm{X}) \cdot E[W \mid \bm{X}] \mid e(\bm{X})] \quad (\bm{X}\text{で条件づけ}) \\ &= E[h(\bm{X}) \cdot e(\bm{X}) \mid e(\bm{X})] \\ &= e(\bm{X}) \cdot E[h(\bm{X}) \mid e(\bm{X})] \end{align} $$
一方、
$$ E[h(\bm{X}) \mid e(\bm{X})] \cdot P(W = 1 \mid e(\bm{X})) = e(\bm{X}) \cdot E[h(\bm{X}) \mid e(\bm{X})] $$
両者が等しいので $W \perp\!\!\!\perp \bm{X} \mid e(\bm{X})$ が成り立ちます。これと元の無視可能性条件を組み合わせると、$(Y(1), Y(0)) \perp\!\!\!\perp W \mid e(\bm{X})$ が得られます。
この定理の実用的意義は非常に大きいです。共変量が高次元であっても、傾向スコアという1つのスカラー値に集約できるため、次元の呪いを回避できます。
IPW推定量
逆確率重み付け（Inverse Probability Weighting）
傾向スコアを用いた IPW推定量（Horvitz-Thompson estimator） を導出します。
ATEの定義に戻ります。
$$ \tau_{\text{ATE}} = E[Y(1)] – E[Y(0)] $$
$E[Y(1)]$ を観測量から求めたいです。条件付き独立性と正値性の下で、
$$ \begin{align} E[Y(1)] &= E[E[Y(1) \mid \bm{X}]] \\ &= E[E[Y(1) \mid W = 1, \bm{X}]] \quad (\because \text{無視可能性}) \\ &= E[E[Y^{\text{obs}} \mid W = 1, \bm{X}]] \end{align} $$
ここで、次の恒等式を示します。
$$ E\left[\frac{W \cdot Y^{\text{obs}}}{e(\bm{X})}\right] = E[Y(1)] $$
導出：
$$ \begin{align} E\left[\frac{W \cdot Y^{\text{obs}}}{e(\bm{X})}\right] &= E\left[E\left[\frac{W \cdot Y^{\text{obs}}}{e(\bm{X})} \;\middle|\; \bm{X}\right]\right] \quad (\text{繰り返し期待値}) \\ &= E\left[\frac{1}{e(\bm{X})} E[W \cdot Y(1) \mid \bm{X}]\right] \quad (\because W=1\text{のとき}Y^{\text{obs}}=Y(1)) \\ &= E\left[\frac{1}{e(\bm{X})} E[W \mid \bm{X}] \cdot E[Y(1) \mid \bm{X}]\right] \quad (\because \text{無視可能性}) \\ &= E\left[\frac{e(\bm{X})}{e(\bm{X})} E[Y(1) \mid \bm{X}]\right] \\ &= E[E[Y(1) \mid \bm{X}]] = E[Y(1)] \end{align} $$
同様に、
$$ E\left[\frac{(1-W) \cdot Y^{\text{obs}}}{1 – e(\bm{X})}\right] = E[Y(0)] $$
したがって、IPW推定量は
$$ \begin{equation} \hat{\tau}_{\text{IPW}} = \frac{1}{n}\sum_{i=1}^{n}\frac{W_i Y_i^{\text{obs}}}{\hat{e}(\bm{X}_i)} – \frac{1}{n}\sum_{i=1}^{n}\frac{(1-W_i) Y_i^{\text{obs}}}{1-\hat{e}(\bm{X}_i)} \end{equation} $$
ここで $\hat{e}(\bm{X}_i)$ は傾向スコアの推定値（ロジスティック回帰などで推定）です。
IPWの直感
IPW推定量の直感は「代表性の補正」です。傾向スコアが低い（処置を受けにくい）のに処置を受けた人は、処置を受けなかった多くの類似の人々を「代表」しているので、大きな重みを与えます。逆に、傾向スコアが高い人は処置を受ける可能性が高いので、小さな重みで十分です。
安定化IPW推定量
基本的なIPW推定量は、傾向スコアが 0 や 1 に近いとき重みが極端に大きくなり、分散が増大するという問題があります。安定化IPW推定量（Hajek estimator） は、重みを正規化することで分散を低減します。
$$ \hat{\tau}_{\text{Hajek}} = \frac{\sum_{i=1}^{n}\frac{W_i Y_i^{\text{obs}}}{\hat{e}(\bm{X}_i)}}{\sum_{i=1}^{n}\frac{W_i}{\hat{e}(\bm{X}_i)}} – \frac{\sum_{i=1}^{n}\frac{(1-W_i) Y_i^{\text{obs}}}{1-\hat{e}(\bm{X}_i)}}{\sum_{i=1}^{n}\frac{(1-W_i)}{1-\hat{e}(\bm{X}_i)}} $$
傾向スコアマッチング
アイデア
傾向スコアマッチングは、傾向スコアが近い処置群個体と対照群個体をペアにして比較する方法です。
処置群の個体 $i$ に対して、対照群の中から傾向スコアが最も近い個体 $j(i)$ を選びます。
$$ j(i) = \underset{j: W_j = 0}{\arg\min}\, |e(\bm{X}_i) – e(\bm{X}_j)| $$
ATTの推定量は、
$$ \hat{\tau}_{\text{match}} = \frac{1}{n_1}\sum_{i: W_i = 1}(Y_i^{\text{obs}} – Y_{j(i)}^{\text{obs}}) $$
ここで $n_1 = \sum_{i=1}^n W_i$ は処置群の個数です。
マッチングの利点と注意点
マッチングの利点は直感的にわかりやすく、モデルの仮定が少ないことです。注意点としては、
- マッチングの質が傾向スコアの推定精度に依存する
- 1対1マッチングではサンプルが半減する
- キャリパー（許容距離）の設定が必要な場合がある
- マッチング後も残余交絡が残る可能性がある
Pythonでの実装
合成データの生成
交絡因子の影響がある合成データを作成し、真のATEと推定値を比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import expit # ロジスティック関数
np.random.seed(42)
n = 2000
# 共変量の生成（交絡因子）
X1 = np.random.normal(0, 1, n) # 連続変数
X2 = np.random.normal(0, 1, n) # 連続変数
X3 = np.random.binomial(1, 0.5, n) # 二値変数
# 処置割り当て（共変量に依存 → 交絡が存在）
logit_ps = -0.5 + 0.8 * X1 + 0.5 * X2 + 0.3 * X3
ps_true = expit(logit_ps) # 真の傾向スコア
W = np.random.binomial(1, ps_true, n)
# 潜在結果の生成
# Y(0) = 2 + X1 + 0.5*X2 + noise
# Y(1) = Y(0) + tau(X) where tau(X) = 3 + 0.5*X1 (異質な処置効果)
Y0 = 2 + X1 + 0.5 * X2 + np.random.normal(0, 1, n)
tau_true = 3 + 0.5 * X1 # 真の個人因果効果
Y1 = Y0 + tau_true
# 観測される結果
Y_obs = W * Y1 + (1 - W) * Y0
# 真のATE
ATE_true = np.mean(tau_true)
print(f"真のATE: {ATE_true:.4f}")
print(f"処置群の割合: {np.mean(W):.4f}")
print(f"処置群の平均Y: {np.mean(Y_obs[W == 1]):.4f}")
print(f"対照群の平均Y: {np.mean(Y_obs[W == 0]):.4f}")
ナイーブ比較
# ナイーブ比較（選択バイアスあり）
naive_estimate = np.mean(Y_obs[W == 1]) - np.mean(Y_obs[W == 0])
print(f"\nナイーブ推定量: {naive_estimate:.4f}")
print(f"真のATE: {ATE_true:.4f}")
print(f"バイアス: {naive_estimate - ATE_true:.4f}")
傾向スコアの推定
from sklearn.linear_model import LogisticRegression
# 傾向スコアの推定（ロジスティック回帰）
X_mat = np.column_stack([X1, X2, X3])
lr = LogisticRegression(max_iter=1000)
lr.fit(X_mat, W)
ps_hat = lr.predict_proba(X_mat)[:, 1]
# 推定された傾向スコアの分布を可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
ax = axes[0]
ax.hist(ps_hat[W == 1], bins=40, density=True, alpha=0.6,
color='coral', label='Treated (W=1)')
ax.hist(ps_hat[W == 0], bins=40, density=True, alpha=0.6,
color='steelblue', label='Control (W=0)')
ax.set_xlabel('Estimated Propensity Score', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title('Distribution of Propensity Scores', fontsize=14)
ax.legend(fontsize=11)
# 真の傾向スコアとの比較
ax = axes[1]
ax.scatter(ps_true, ps_hat, alpha=0.2, s=5)
ax.plot([0, 1], [0, 1], 'r--', linewidth=2)
ax.set_xlabel('True Propensity Score', fontsize=12)
ax.set_ylabel('Estimated Propensity Score', fontsize=12)
ax.set_title('True vs Estimated PS', fontsize=14)
plt.tight_layout()
plt.savefig("propensity_score.png", dpi=150, bbox_inches="tight")
plt.show()
IPW推定
# IPW推定量（Horvitz-Thompson型）
ipw_treated = np.sum(W * Y_obs / ps_hat) / n
ipw_control = np.sum((1 - W) * Y_obs / (1 - ps_hat)) / n
ate_ipw = ipw_treated - ipw_control
# 安定化IPW推定量（Hajek型）
w_treated = W / ps_hat
w_control = (1 - W) / (1 - ps_hat)
ate_hajek = np.sum(w_treated * Y_obs) / np.sum(w_treated) \
- np.sum(w_control * Y_obs) / np.sum(w_control)
print(f"IPW推定量 (HT): {ate_ipw:.4f}")
print(f"IPW推定量 (Hajek): {ate_hajek:.4f}")
print(f"真のATE: {ATE_true:.4f}")
傾向スコアマッチング
from scipy.spatial import KDTree
# 傾向スコアマッチング（1対1、最近傍）
treated_idx = np.where(W == 1)[0]
control_idx = np.where(W == 0)[0]
# KDTreeで最近傍検索
tree = KDTree(ps_hat[control_idx].reshape(-1, 1))
distances, indices = tree.query(ps_hat[treated_idx].reshape(-1, 1), k=1)
# マッチされた対照群のインデックス
matched_control_idx = control_idx[indices.ravel()]
# ATT推定
ate_matching = np.mean(Y_obs[treated_idx] - Y_obs[matched_control_idx])
# 真のATT
att_true = np.mean(tau_true[W == 1])
print(f"\nマッチング推定量 (ATT): {ate_matching:.4f}")
print(f"真のATT: {att_true:.4f}")
結果の比較と可視化
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import expit
from sklearn.linear_model import LogisticRegression
from scipy.spatial import KDTree
# --- 多数回のシミュレーションで各推定量のバイアスと分散を評価 ---
n = 1000
n_sim = 500
results = {'naive': [], 'ipw_ht': [], 'ipw_hajek': [], 'matching': []}
for sim in range(n_sim):
# データ生成
X1 = np.random.normal(0, 1, n)
X2 = np.random.normal(0, 1, n)
X3 = np.random.binomial(1, 0.5, n)
logit_ps = -0.5 + 0.8 * X1 + 0.5 * X2 + 0.3 * X3
ps_true = expit(logit_ps)
W = np.random.binomial(1, ps_true, n)
Y0 = 2 + X1 + 0.5 * X2 + np.random.normal(0, 1, n)
tau = 3 + 0.5 * X1
Y1 = Y0 + tau
Y_obs = W * Y1 + (1 - W) * Y0
ATE_true_i = np.mean(tau)
# ナイーブ
results['naive'].append(np.mean(Y_obs[W == 1]) - np.mean(Y_obs[W == 0]))
# 傾向スコア推定
X_mat = np.column_stack([X1, X2, X3])
lr = LogisticRegression(max_iter=1000)
lr.fit(X_mat, W)
ps_hat = lr.predict_proba(X_mat)[:, 1]
ps_hat = np.clip(ps_hat, 0.01, 0.99) # 極端な値のクリッピング
# IPW (HT)
ipw_t = np.sum(W * Y_obs / ps_hat) / n
ipw_c = np.sum((1 - W) * Y_obs / (1 - ps_hat)) / n
results['ipw_ht'].append(ipw_t - ipw_c)
# IPW (Hajek)
wt = W / ps_hat
wc = (1 - W) / (1 - ps_hat)
results['ipw_hajek'].append(
np.sum(wt * Y_obs) / np.sum(wt) - np.sum(wc * Y_obs) / np.sum(wc)
)
# マッチング
t_idx = np.where(W == 1)[0]
c_idx = np.where(W == 0)[0]
if len(c_idx) > 0 and len(t_idx) > 0:
tree = KDTree(ps_hat[c_idx].reshape(-1, 1))
_, idx = tree.query(ps_hat[t_idx].reshape(-1, 1), k=1)
matched_c = c_idx[idx.ravel()]
results['matching'].append(np.mean(Y_obs[t_idx] - Y_obs[matched_c]))
# 真のATE（母集団レベル）
ATE_true_pop = 3.0 # E[3 + 0.5*X1] = 3 (X1~N(0,1)なので E[X1]=0)
fig, ax = plt.subplots(figsize=(12, 6))
labels = ['Naive', 'IPW (HT)', 'IPW (Hajek)', 'Matching']
data = [results['naive'], results['ipw_ht'], results['ipw_hajek'], results['matching']]
colors = ['#ff6b6b', '#4ecdc4', '#45b7d1', '#96ceb4']
bp = ax.boxplot(data, labels=labels, patch_artist=True, widths=0.6)
for patch, color in zip(bp['boxes'], colors):
patch.set_facecolor(color)
patch.set_alpha(0.7)
ax.axhline(y=ATE_true_pop, color='red', linestyle='--', linewidth=2,
label=f'True ATE = {ATE_true_pop:.1f}')
ax.set_ylabel('Estimated ATE', fontsize=12)
ax.set_title(f'Comparison of ATE Estimators ({n_sim} simulations, n={n})',
fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3, axis='y')
# バイアスとRMSEの表示
for i, (label, vals) in enumerate(zip(labels, data)):
vals_arr = np.array(vals)
bias = np.mean(vals_arr) - ATE_true_pop
rmse = np.sqrt(np.mean((vals_arr - ATE_true_pop)**2))
ax.text(i + 1, ax.get_ylim()[0] + 0.1,
f'Bias={bias:.3f}\nRMSE={rmse:.3f}',
ha='center', fontsize=9, bbox=dict(boxstyle='round',
facecolor='wheat', alpha=0.5))
plt.tight_layout()
plt.savefig("ate_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
シミュレーション結果から、以下のことが確認できます。
- ナイーブ推定量 は交絡によるバイアスが大きく、ATEを過大推定します。これは処置を受けやすい人がもともと結果が良い傾向にあるためです。
- IPW推定量 はバイアスが大幅に低減されます。Hajek型はHT型より分散が小さい傾向があります。
- マッチング推定量 もバイアスの低減に効果的ですが、ATT寄りの推定値を返す傾向があります。
まとめ
本記事では、因果推論の基礎理論であるルービンの因果モデルについて解説しました。
- 因果推論の根本問題として、同一個体の反事実が観測できないことを説明しました
- 潜在結果 $Y(1), Y(0)$ を用いてATE、ATT、ATUを厳密に定義しました
- SUTVA仮定（干渉なし、処置の一貫性）の意味と重要性を説明しました
- ナイーブな比較では選択バイアスが生じることを数学的に示しました
- 強い無視可能性条件の下で因果効果が識別可能になることを証明しました
- ランダム化実験が因果推論で最も信頼できる理由を説明しました
- 傾向スコアの性質（次元削減）を証明し、IPW推定量を導出しました
- Pythonで合成データを用いて各推定量のバイアスとRMSEを比較検証しました
因果推論は観察データから介入効果を推定するための不可欠な枠組みです。次のステップとして、以下の記事も参考にしてください。