実データ事例と推定パイプライン

🎓 レベル：発展　|　重要度：B（標準）

📎 前提：傾向スコア　|　二重頑健推定AIPW　|　因果推論のチェックリスト

要点（BLUF）

因果分析は単発の推定ではなく 問い → DAG → 識別 → 推定 → 反証/感度 の一気通貫パイプライン。この順序を 1 つの擬似データで最初から最後まで通す。
真の効果 ATE_true=30（職業訓練が年収を +30 万円）を仕込み、素朴比較は 71.7 と大きく上振れするが、回帰調整・IPW・AIPW(DML) はいずれも約 30 を回収する。識別（バックドア）が先、推定はその後。
推定で止めず反証する：プラセボ処置（効果が 0 に消える）・部分集団の安定性・E-value（未観測交絡への感度）。仕上げに dowhy の model→identify→estimate→refute も紹介（要最新確認）。

1. パイプライン全体図

因果推論の章を貫く背骨は「識別を済ませてから推定し、推定したら反証する」。手順を図にする。

flowchart TB
  Q["①問い: 訓練Tは年収Yを上げるか"] --> D["②DAG: 交絡C(技能)・ageを描く"]
  D --> I["③識別: バックドア基準で調整集合を確定"]
  I --> E["④推定: 回帰調整 / IPW / AIPW(DML)"]
  E --> R["⑤反証: プラセボ・部分集団・E-value"]
  R --> C["結論 + 残る不確実性の明示"]

各段は本テキストの章に対応する：②③はバックドア基準と識別、④は回帰による調整とその限界・逆確率重み付けIPW・二重頑健推定AIPW・Double/Debiased Machine Learning（DML）、⑤は未観測交絡への感度分析・因果推論のチェックリスト。

2. 問いとデータ生成構造

問い：職業訓練プログラム $T$ （受講 1／非受講 0）は年収 $Y$ （万円）を上げるか。真の効果は $\mathrm{ATE}=30$ 万円と仕込む。

交絡構造：訓練前の 技能スコア $C$ が高い人ほど (a) 訓練に応募しやすく、(b) もともと年収も高い。年齢 $age$ も両方に効く。だから $C,age$ は共通原因＝交絡で、放置すると素朴比較が上振れする。

flowchart LR
  C["交絡C: 技能"] --> T["処置T: 訓練受講"]
  C --> Y["結果Y: 年収"]
  AGE["交絡: 年齢"] --> T
  AGE --> Y
  T -->|"真の効果 +30"| Y

3. 識別——バックドア基準で調整集合を確定する

推定の前に識別を済ませる。 $T$ から $Y$ へのバックドアパスは $T\leftarrow C\to Y$ と $T\leftarrow age\to Y$ の 2 本。調整集合 $\{C, age\}$ が バックドア基準を満たす（すべてのバックドアパスを塞ぎ、子孫を含まない）。よって識別の仮定の条件付き交換可能性のもとで ATE はバックドア調整公式で識別される。

\mathrm{ATE} = E_{C,age}\!\Big[\,E[Y\mid T=1, C, age] - E[Y\mid T=0, C, age]\,\Big]

加えて正値性（どの $C,age$ でも受講・非受講の両方が起こりうる）を後でデータで確認する。これが満たされて初めて推定に進める——識別なき推定は数字遊びである。

4. 推定——素朴比較は外し、調整は真値を当てる

擬似データを作り、(1) 素朴な平均差、(2) 回帰調整、(3) 傾向スコア IPW、(4) AIPW（二重頑健＝ATE に対する DML）の 4 つを並べる。IPW と AIPW の推定量は

\hat{\mathrm{ATE}}_{\mathrm{IPW}} = \frac{\sum_i T_i Y_i/\hat e(X_i)}{\sum_i T_i/\hat e(X_i)} - \frac{\sum_i (1-T_i) Y_i/(1-\hat e(X_i))}{\sum_i (1-T_i)/(1-\hat e(X_i))}

\hat{\mathrm{ATE}}_{\mathrm{AIPW}} = \frac1n\sum_i \Big[\hat\mu_1(X_i)-\hat\mu_0(X_i) + \frac{T_i\,(Y_i-\hat\mu_1)}{\hat e(X_i)} - \frac{(1-T_i)(Y_i-\hat\mu_0)}{1-\hat e(X_i)}\Big]

ここで $\hat e$ は傾向スコア、 $\hat\mu_t$ は結果モデル。AIPW は結果モデルか傾向スコアのどちらかが当たれば一致（二重頑健・二重頑健推定AIPW）。交差適合を入れるとDouble/Debiased Machine Learning（DML）の ATE 版になる。

# === 推定パイプライン（擬似データ生成→識別→4推定量） ===
import warnings
warnings.filterwarnings("ignore")   # scipy/sklearnの無害なソルバ警告を抑制
import numpy as np
import pandas as pd

rng = np.random.default_rng(20)
ATE_true = 30.0      # 職業訓練の真の効果（年収 万円）
n = 4000

# 交絡 C(技能)・age が処置にも結果にも効く構造
skill = rng.normal(0, 1, n)
age = rng.uniform(22, 60, n)
logit = -0.4 + 1.2 * skill + 0.03 * (age - 40)
prop_true = 1 / (1 + np.exp(-logit))
T = (rng.uniform(0, 1, n) < prop_true).astype(int)
Y = 300 + ATE_true * T + 40 * skill + 1.5 * (age - 40) + rng.normal(0, 20, n)
df = pd.DataFrame({"T": T, "skill": skill, "age": age, "Y": Y})

# 正値性チェック（傾向が0/1に張り付いていないか）
print("処置率 P(T=1)      :", round(T.mean(), 3))
print("傾向スコアの範囲   : [", round(prop_true.min(), 3), ",", round(prop_true.max(), 3), "]")

# --- 推定1：素朴な平均差（交絡で上振れするはず） ---
naive = Y[T == 1].mean() - Y[T == 0].mean()
print("ATE_true           :", ATE_true)
print("(1)素朴な平均差    :", round(naive, 2))

# --- 推定2：回帰調整（OLSで交絡を共変量に） ---
import statsmodels.api as sm
Xmat = sm.add_constant(df[["T", "skill", "age"]])
ols = sm.OLS(df["Y"], Xmat).fit()
print("(2)回帰調整(OLS)   :", round(ols.params["T"], 2))

# --- 推定3：傾向スコアによるIPW ---
from sklearn.linear_model import LogisticRegression
covars = df[["skill", "age"]].values
ps_model = LogisticRegression().fit(covars, df["T"].values)
e_hat = ps_model.predict_proba(covars)[:, 1]
e_hat = np.clip(e_hat, 0.02, 0.98)
w = np.where(df["T"] == 1, 1 / e_hat, 1 / (1 - e_hat))
ipw = np.sum(w * (df["T"] == 1) * df["Y"]) / np.sum(w * (df["T"] == 1)) \
    - np.sum(w * (df["T"] == 0) * df["Y"]) / np.sum(w * (df["T"] == 0))
print("(3)IPW(傾向スコア) :", round(ipw, 2))

# --- 推定4：AIPW＋交差適合（二重頑健＝ATEに対するDML） ---
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold
kf = KFold(n_splits=5, shuffle=True, random_state=0)
mu1 = np.zeros(n); mu0 = np.zeros(n); e_cf = np.zeros(n)
for tr, te in kf.split(covars):
    out1 = LinearRegression().fit(covars[tr][df["T"].values[tr] == 1], Y[tr][df["T"].values[tr] == 1])
    out0 = LinearRegression().fit(covars[tr][df["T"].values[tr] == 0], Y[tr][df["T"].values[tr] == 0])
    mu1[te] = out1.predict(covars[te])
    mu0[te] = out0.predict(covars[te])
    psm = LogisticRegression().fit(covars[tr], df["T"].values[tr])
    e_cf[te] = np.clip(psm.predict_proba(covars[te])[:, 1], 0.02, 0.98)
Tarr = df["T"].values
psi = (mu1 - mu0) + Tarr * (Y - mu1) / e_cf - (1 - Tarr) * (Y - mu0) / (1 - e_cf)
aipw = psi.mean()
aipw_se = psi.std(ddof=1) / np.sqrt(n)
print("(4)AIPW/DML        :", round(aipw, 2), " (SE", round(aipw_se, 2), ")")

出力：

処置率 P(T=1)      : 0.436
傾向スコアの範囲   : [ 0.009 , 0.974 ]
ATE_true           : 30.0
(1)素朴な平均差    : 71.69
(2)回帰調整(OLS)   : 29.95
(3)IPW(傾向スコア) : 30.66
(4)AIPW/DML        : 30.69  (SE 0.87 )

出力の意味：傾向スコアの範囲が $[0.009, 0.974]$ で 0/1 に張り付いておらず正値性は概ね良好（識別の前提が成立）。素朴な平均差 71.69 は真値 30 の倍以上に上振れ——高技能者ほど受講しかつ高年収という交絡がそのまま乗っている。一方 回帰調整 29.95・IPW 30.66・AIPW/DML 30.69 はいずれも真値 30 を回収した。3 つの別経路（結果モデル・傾向スコア・両者の併用）が一致することが信頼の根拠になる。AIPW は標準誤差も出るので 95% 区間は約 $30.7\pm1.7$ 。

5. 反証——推定で止めない

点推定が出ても終わりではない。わざと壊れる操作を加えて、壊れるべき時に壊れるかを確かめる（反証）。3 つを自己完結のブロックで行う。E-value は連続アウトカムを標準化効果 $d=\mathrm{ATE}/\mathrm{SD}(Y)$ からリスク比へ近似変換し、

\mathrm{RR}\approx \exp(0.91\,d), \qquad E\text{-value} = \mathrm{RR} + \sqrt{\mathrm{RR}\,(\mathrm{RR}-1)}

で計算する（VanderWeele & Ding 2017・未観測交絡への感度分析）。

# === 反証（プラセボ処置・部分集団の安定性・E-value感度分析） ===
import numpy as np
import pandas as pd
import statsmodels.api as sm

rng = np.random.default_rng(20)
ATE_true = 30.0
n = 4000
skill = rng.normal(0, 1, n)
age = rng.uniform(22, 60, n)
logit = -0.4 + 1.2 * skill + 0.03 * (age - 40)
T = (rng.uniform(0, 1, n) < 1 / (1 + np.exp(-logit))).astype(int)
Y = 300 + ATE_true * T + 40 * skill + 1.5 * (age - 40) + rng.normal(0, 20, n)
df = pd.DataFrame({"T": T, "skill": skill, "age": age, "Y": Y})


def adjusted_effect(data):
    # 回帰調整での処置効果（Tの係数）を返す
    X = sm.add_constant(data[["T", "skill", "age"]])
    fit = sm.OLS(data["Y"], X).fit()
    return fit.params["T"]


point = adjusted_effect(df)
print("調整済み効果(点推定):", round(point, 2))

# --- 反証1：プラセボ処置（Tをシャッフル＝因果を持たない偽の処置）→ 効果は0近傍のはず ---
df_placebo = df.copy()
df_placebo["T"] = rng.permutation(df["T"].values)
placebo_effect = adjusted_effect(df_placebo)
print("反証1 プラセボ効果  :", round(placebo_effect, 2), "(0近傍なら健全)")

# --- 反証2：部分集団の安定性（ランダムに半分→再推定）→ 点推定の近傍で安定のはず ---
idx = rng.permutation(n)[: n // 2]
subset_effect = adjusted_effect(df.iloc[idx])
print("反証2 部分集団効果  :", round(subset_effect, 2), "(点推定の近傍なら健全)")

# --- 反証3：E-value（未観測交絡への感度） ---
# 連続アウトカムは標準化効果 d を RR に近似変換：RR=exp(0.91*d)、E=RR+sqrt(RR*(RR-1))
d = point / df["Y"].std()
RR = np.exp(0.91 * d)
if RR < 1:
    RR = 1 / RR
e_value = RR + np.sqrt(RR * (RR - 1))
print("標準化効果 d        :", round(d, 3))
print("近似リスク比 RR     :", round(RR, 3))
print("E-value             :", round(e_value, 2))

出力：

調整済み効果(点推定): 29.95
反証1 プラセボ効果  : 0.41 (0近傍なら健全)
反証2 部分集団効果  : 31.06 (点推定の近傍なら健全)
標準化効果 d        : 0.539
近似リスク比 RR     : 1.634
E-value             : 2.65

出力の意味：

反証1 プラセボ：処置 $T$ を無作為シャッフルし因果を消した偽処置に同じ推定をかけると、効果は 0.41 ≈ 0。もしここで大きな「効果」が出たら、推定手続きが交絡やバグで偽陽性を作っている証拠になる。0 近傍なので手続きは健全。
反証2 部分集団：データを半分に絞って再推定しても 31.06 と点推定 29.95 の近傍で安定。特定の部分集団に依存した不安定な結果ではない。
反証3 E-value ＝ 2.65：意味は「測定済み共変量を超えて、未観測交絡が処置とも結果とも少なくともリスク比 2.65 で結びついていなければ、この効果を完全には説明し去れない」。2.65 はそこそこ大きく、弱い交絡では覆らない＝頑健性は中程度と読む。逆に E-value が 1.2 程度なら、ありふれた未観測交絡で消えうる脆い結果ということ。

6. dowhy で 4 ステップを自動化（要最新確認）

dowhy（py-why）は model → identify → estimate → refute の 4 ステップを統一 API にしたライブラリ。本ノートで手作業した流れをそのまま自動化できる。API は動きやすいので要最新確認、未導入なら pip install dowhy（依存が重く最新環境では衝突しうるので隔離環境推奨）。本ノートの主経路は上の手実装で完結しており、dowhy は同じ流れの確認用。

# 要 pip install dowhy・要最新確認（本ノートの主経路は上の手実装で完結。下は流れの対応図）
from dowhy import CausalModel

model = CausalModel(data=df, treatment="T", outcome="Y",
                    common_causes=["skill", "age"])          # ①model: DAG/交絡を宣言
estimand = model.identify_effect()                           # ②identify: バックドアで識別
estimate = model.estimate_effect(
    estimand, method_name="backdoor.propensity_score_weighting")  # ③estimate: IPWで推定
refute = model.refute_estimate(
    estimand, estimate, method_name="placebo_treatment_refuter")  # ④refute: プラセボで反証

identify_effect() が本ノート §3 の識別、estimate_effect(method_name=...) が §4 の推定、refute_estimate(method_name=...) が §5 の反証に対応する。ライブラリは流れを固定してくれるが、識別の妥当性（DAG が正しいか）は人間の責任で、そこは自動化されない。

⚠️ よくある誤解・落とし穴

識別を飛ばして推定から入らない。DAG とバックドアの確認（§3）を省くと、どの共変量を入れるべきか（あるいは入れてはいけないか＝過剰調整とMバイアス）が決まらず、数字は出ても意味が定まらない。
複数推定量の一致は「正しさ」でなく「同じ仮定の整合」。回帰・IPW・AIPW が揃っても、全員が同じ未観測交絡で同じ方向に外れているかもしれない。だから §5 の感度分析（E-value）が要る。
プラセボで効果が消えないのは赤信号。手続き・データ整形・交絡調整のどこかが偽の関連を作っている。点推定の前にプラセボを通す癖をつける。
E-value は「未観測交絡の不在」を証明しない。「どれだけ強い交絡なら結論が覆るか」を示すだけ。小さな E-value は脆さの警告、大きくても交絡が無い保証にはならない。
正値性の確認を省かない。傾向スコアが 0/1 に張り付く層では IPW の重みが爆発し、推定が不安定になる（逆確率重み付けIPW）。範囲の印字は儀式でなく必須。