MLをそのまま使うと因果を誤る理由

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：二重頑健推定AIPW　|　識別の仮定　|　数理：正則化（Ridge・Lasso・Elastic Net）（機械学習）

要点（BLUF）

高次元の交絡を機械学習（Lasso 等）で「予測」して引き算する素朴なプラグインは、処置効果を系統的に外す。予測が上手いことと因果が正しいことは別物。
原因は推定式が Neyman 直交（orthogonal）でないこと。正則化で生じた撹乱母数（nuisance）のバイアスが、処置効果に一次（線形）で漏れる。
擬似データで、Lasso で結果だけ調整して生の処置に回帰すると真値 ATE_true=2.0 が 0.89 に潰れることを確認する。直交化（Double/Debiased Machine Learning（DML））で回収できる。

なぜ「ML で交絡を調整」したいのか

交絡変数 $X$ が高次元（数十〜数千次元）だと、回帰による調整とその限界で見た線形回帰では関数形を取り違えやすい。そこで「 $E[Y \mid X]$ や $E[T \mid X]$ のような撹乱母数（nuisance function）を柔軟な機械学習で当て、その分を取り除けば交絡を消せるはず」という発想が自然に出てくる。

しかしここに落とし穴がある。機械学習の予測器をそのまま因果推定の式に差し込む（プラグインする）と、推定値が真の効果から系統的にずれる。 その理由を、構造を完全に把握できる擬似データで突き止める。

部分線形モデルと因果構造

処置 $T$ と結果 $Y$ の関係を、次の**部分線形モデル（partially linear model）**で考える。

Y = \theta\, T + g(X) + \varepsilon, \qquad T = m(X) + \eta

$\theta$ ：知りたい処置効果（ここでは定数の ATE）。
$g(X) = E[Y \mid X, T=0]$ 相当の交絡の本体（ $X$ が $Y$ に与える複雑な影響）。
$m(X) = E[T \mid X]$ ： $X$ から処置への影響（傾向スコアの連続版）。
$\varepsilon, \eta$ ：それぞれ平均 0 のノイズ。 $X$ は処置と結果の共通原因だから交絡が生じる。

flowchart LR
    X["高次元交絡 X（多数の共変量）"] --> T["処置 T（用量）"]
    X --> Y["結果 Y"]
    T -->|"効果 θ"| Y

$X \to T$ と $X \to Y$ の両方が伸びているので、 $T$ と $Y$ の素朴な比較にはバックドアパス $T \leftarrow X \to Y$ が混ざる。

識別の仮定（ここを満たして初めて θ は因果）

推定の前に、 $\theta$ が因果効果として識別できる条件を明示する（識別の仮定）。

条件付き交換可能性： $X$ で条件づければ未観測交絡はない（ $\{Y(0),Y(1)\} \perp T \mid X$ ）。バックドアパスは $X$ で全て塞がる。
正値性：あらゆる $X$ の値で処置の割り当てに分散がある（ $0 < \Pr(T \mid X)$ 相当、連続処置なら $\mathrm{Var}(\eta) > 0$ ）。
SUTVA：個体間の干渉がなく処置の版が一つ。

この 3 つが成り立つ前提で初めて、 $\theta$ は構造方程式の係数＝因果効果になる。以下で問題にするのは「識別できた $\theta$ をどう推定するか」という推定の話であって、識別の失敗（未観測交絡）とは別の問題であることに注意する。

素朴な ML プラグインはなぜ外れるか

識別の仮定の下では、Frisch–Waugh–Lovell の関係から $\theta$ は両方を残差化した回帰係数として書ける。

\theta = \frac{E\big[(Y - E[Y\mid X])(T - E[T\mid X])\big]}{E\big[(T - E[T\mid X])^2\big]}

ここで「結果だけを $\ell(X)=E[Y\mid X]$ で残差化し、生の $T$ に回帰する」という素朴なプラグインを考えると何が起きるか。 $\ell(X) = \theta\, m(X) + g(X)$ なので、 $\hat\ell \to \ell$ と当てられたとして残差は

\tilde Y = Y - \ell(X) = \theta\big(T - m(X)\big) + \varepsilon = \theta\,\eta + \varepsilon

これを生の $T = m(X)+\eta$ に回帰すると、

\hat\theta_{\text{naive}} = \frac{E[\tilde Y\, T]}{E[T^2]} = \theta\,\frac{\mathrm{Var}(\eta)}{\mathrm{Var}(m(X)) + \mathrm{Var}(\eta)} \;<\; \theta

処置の予測可能な分散 $\mathrm{Var}(m(X))$ の分だけ、効果が縮小（減衰）する。 処置を残差化し損ねた $m(X)$ が、そのまま $\hat\theta$ のバイアスとして「漏れた」のである。次のコードで、真値 2.0 が理論上 $2 \times \frac{1}{\mathrm{Var}(m)+1} \approx 0.77$ まで潰れることを確かめる。

コード：高次元交絡で素朴な調整が外れる

下のコードは、20 次元の交絡を仕込んだ部分線形モデルから、(1) 無調整の回帰と (2) 「Lasso で結果だけ調整して生の処置に回帰」する非直交プラグインを計算する。ATE_true = 2.0 を当てられるかを見る。

# === 高次元交絡を仕込み、無調整と「非直交プラグイン」が外すことを確かめる ===
import numpy as np
from sklearn.linear_model import LassoCV, LinearRegression
from sklearn.model_selection import KFold

rng = np.random.default_rng(42)
ATE_true = 2.0
n, p = 500, 20
X = rng.standard_normal((n, p))
beta  = 1.0 / np.arange(1, p + 1)        # g(X) の係数（X→結果Y の交絡経路）
gamma = 1.0 / np.arange(1, p + 1)        # m(X) の係数（X→処置T の交絡経路）
m = X @ gamma
T = m + rng.standard_normal(n)                         # 連続処置（用量）: T = m(X) + η
Y = ATE_true * T + X @ beta + rng.standard_normal(n)   # Y = θ·T + g(X) + ε

# (1) 無調整：Y を T だけに回帰（交絡を放置）
b_unadj = LinearRegression().fit(T.reshape(-1, 1), Y).coef_[0]

# (2) 非直交プラグイン：Lasso で E[Y|X] を当て、Y から引いてから「生の T」に回帰
kf = KFold(n_splits=5, shuffle=True, random_state=0)
y_hat = np.zeros(n)
for idx_tr, idx_te in kf.split(X):
    g_model = LassoCV(cv=5, random_state=0).fit(X[idx_tr], Y[idx_tr])
    y_hat[idx_te] = g_model.predict(X[idx_te])
Y_resid = Y - y_hat
b_plugin = np.mean(T * Y_resid) / np.mean(T * T)

var_m = gamma @ gamma                      # Var(m(X)) = Σ γ_j²（X は標準正規・独立）
shrink = 1.0 / (var_m + 1.0)               # 理論上の減衰率 Var(η)/Var(T)
print("ATE_true              :", ATE_true)
print("(1) 無調整 Y~T        :", round(b_unadj, 3))
print("(2) 非直交プラグイン  :", round(b_plugin, 3))
print("理論減衰率 Var(η)/Var(T):", round(shrink, 3),
      " → 予測値 θ×減衰 =", round(ATE_true * shrink, 3))

出力は次のとおり。

ATE_true              : 2.0
(1) 無調整 Y~T        : 2.636
(2) 非直交プラグイン  : 0.894
理論減衰率 Var(η)/Var(T): 0.385  → 予測値 θ×減衰 = 0.77

無調整は 2.64 と過大（交絡 $X$ が処置と結果を同方向に押すため上振れ）。一方、Lasso でちゃんと結果を調整したはずの非直交プラグインは 0.89 と、真値 2.0 を大きく下回る。理論減衰 0.77 とほぼ一致しており、ずれが偶然ではなく構造的なバイアスであることが分かる。「ML で予測してから引く」だけでは因果は当たらない。

コード：100 回反復で分布の中心を見る

1 回の値が偶然でないことを示すため、データ生成を 100 回繰り返し、各推定量の分布の中心が真値 2.0 に来るかを見る。あわせて、次ノートで扱う DML（処置も残差化して交差適合） を先取りで重ねる。

# === 100回反復して、各推定量の分布（中心が真値2.0か）を見る ===
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
from sklearn.linear_model import LassoCV, LinearRegression
from sklearn.model_selection import KFold

ATE_true = 2.0
n, p = 400, 20
beta = 1.0 / np.arange(1, p + 1)
gamma = 1.0 / np.arange(1, p + 1)

def generate(rng):
    X = rng.standard_normal((n, p))
    T = X @ gamma + rng.standard_normal(n)
    Y = ATE_true * T + X @ beta + rng.standard_normal(n)
    return X, T, Y

def crossfit_resid(X, Z):                  # E[Z|X] を交差適合で当てて残差を返す
    kf = KFold(n_splits=5, shuffle=True, random_state=0)
    pred = np.zeros(len(Z))
    for idx_tr, idx_te in kf.split(X):
        model = LassoCV(cv=3, random_state=0).fit(X[idx_tr], Z[idx_tr])
        pred[idx_te] = model.predict(X[idx_te])
    return Z - pred

results = {"無調整": [], "非直交プラグイン": [], "DML(直交+交差適合)": []}
for rep in range(100):
    rng = np.random.default_rng(1000 + rep)
    X, T, Y = generate(rng)
    results["無調整"].append(LinearRegression().fit(T.reshape(-1, 1), Y).coef_[0])
    Y_res = crossfit_resid(X, Y)
    results["非直交プラグイン"].append(np.mean(T * Y_res) / np.mean(T * T))
    T_res = crossfit_resid(X, T)           # 処置も残差化するのが直交化（次ノート）
    results["DML(直交+交差適合)"].append(np.mean(T_res * Y_res) / np.mean(T_res * T_res))

for name, vals in results.items():
    vals = np.array(vals)
    print(f"{name:18s} 平均={vals.mean():.3f}  バイアス={vals.mean()-ATE_true:+.3f}")

plt.figure(figsize=(8, 4))
for name, vals in results.items():
    plt.hist(vals, bins=20, alpha=0.6, label=name)
plt.axvline(ATE_true, color="black", linestyle="--", label="真値 ATE_true=2.0")
plt.xlabel("処置効果の推定値"); plt.ylabel("頻度")
plt.title("ML調整の素朴なプラグインは中心がずれる（100反復）")
plt.legend(); plt.tight_layout(); plt.show()

出力は次のとおり。

無調整                平均=2.621  バイアス=+0.621
非直交プラグイン           平均=0.812  バイアス=-1.188
DML(直交+交差適合)       平均=2.000  バイアス=+0.000

ヒストグラムでは、無調整は右に、非直交プラグインは左に外れ、どちらも真値 2.0 に中心が乗らない。一方で DML だけが 2.000（バイアス +0.000）で真値に重なる。「処置も残差化する」というたった一手の違いが、系統的バイアスを消すことを次ノートで掘り下げる。

直観：なぜ「直交性」が決定的か

撹乱母数 $\eta = (\ell, m)$ を含む推定式（モーメント条件） $\psi(W;\theta,\eta)$ が Neyman 直交であるとは、真値 $(\theta_0,\eta_0)$ における $\eta$ 方向の微分（ガトー微分）がゼロであること、すなわち撹乱母数の推定誤差に対して一次で鈍感であることを言う。

機械学習の予測器は、バイアス–バリアンスのトレードオフ（汎化と過学習・バイアスバリアンス分解）と正則化（正則化（Ridge・Lasso・Elastic Net））により、わざと縮小バイアスを残す。その収束は $\|\hat\eta-\eta_0\| = O(n^{-c})$ で、ふつう $c < 1/2$ （ $\sqrt n$ より遅い）。
非直交な式では、この誤差が $\hat\theta-\theta$ に一次でそのまま乗る： $\;|\hat\theta-\theta| \sim O(n^{-c})$ 。上の減衰バイアスはまさにこれで、Lasso の縮小が $\hat\theta$ に漏れた姿である。
直交な式にすると、誤差は二次でしか効かない： $\;|\hat\theta-\theta| \lesssim \|\hat\ell-\ell\|\cdot\|\hat m-m\| = O(n^{-2c})$ 。 $c>1/4$ なら $\sqrt n$ より速く消えるので、遅い ML 撹乱母数を使っても $\hat\theta$ は $\sqrt n$ 一致になる。

つまり「結果も処置も残差化する」直交モーメントにすれば、二つの遅い誤差の積になって無視できる。これが次ノートの DML の核心である。

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

「予測精度が高ければ因果も正しい」は誤り。 $E[Y\mid X]$ をどれだけ正確に当てても、推定式が非直交なら処置効果はずれる。予測の良さ（ $R^2$ ）と因果の正しさは別の評価軸。
同じ標本で撹乱母数と $\theta$ を推定すると過学習バイアスが乗る。 上のコードは交差適合（out-of-fold 予測）で過学習分はすでに除いてある。それでも残るのが「非直交性」由来のバイアスで、両者は別物（過学習バイアスは訓練・検証・テストと交差検証の発想で、非直交性は直交化で消す）。
ここでの議論は識別が成立している前提。 高次元 $X$ を入れても未観測交絡があれば交換可能性は崩れ、どんな手法でも因果は出ない。「ML で変数をたくさん入れたから大丈夫」は識別の保証にならない。
過剰調整に注意。 $X$ に媒介変数（ $T \to X \to Y$ の中間）や衝突点を混ぜると、機械学習で「上手く」予測できても効果はバイアスする。何を調整集合に入れるかは DAG で決める問題で、予測精度では決まらない。