この連載について
因果推論を「本を読んだ」で終わらせず、自分の言葉で説明でき、コードで再現できる状態まで落とす30日連載です。直前のDay 23では、S/T/X-learnerという3つのメタラーナーで、手持ちの機械学習を組み合わせて「誰に効くか」τ(x) を推定しました。今日はPhase 4のもう1つの系譜です。Day 22の地図で言うと、Day 23・25・26は「誰に効くか（CATE）」の系譜、今日のDay 24は「平均効果を高次元の交絡の下で頑健に出す」系譜。主役は Double Machine Learning（DML）です。図はすべて合成データで、真の効果 θ=0.5 を私たちが握った状態で見せます。
TL;DR（3行）
- 機械学習に丸ごと任せると正則化バイアス。 Lassoに処置wと共変量200列を渡してwの係数を読むと、真値0.5に対し 0.581。モンテカルロ100回では平均0.563、100回中100回が上振れします。正則化が交絡の係数を縮める → 縮んだ分の交絡がwの係数へ漏れる、が機構です。
- DMLはFWL（Day 16）のML版。 wとyをそれぞれMLで予測して残差を取り、残差×残差の単回帰。同じデータで 0.510、100回の平均 0.512 と真値を囲みます。ニュイサンスを少し間違えても1次では効かない、これが Neyman直交性 の直感です。
- 交差適合（cross-fitting）が生命線。 自分のデータで当てはめた残差は過学習で消滅し（残差分散 0.0001）、推定値は1.523へ暴走。5分割の交差適合で0.536に戻り、econml LinearDML の95%CI [0.466, 0.607] が真値を含みます。
今日の問い
交絡らしき共変量が200列ある。素のOLSに全部積むのは心もとない。予測が得意な機械学習に交絡の後始末を任せたい。
では、機械学習が推定した「処置の係数」を、そのまま因果効果と読んでよいか？
答えは「そのままでは読めない」です。機械学習は予測を当てるために、係数を歪める癖（正則化）を意図的に持っています。今日のゴールは、この癖がなぜ因果推定に漏れるのかを機構で言語化し、それを打ち消す2つの手当て、残差化と交差適合をコードで組むことです。組み上がったものがDMLで、正体はDay 16のFWL定理の現代版だと分かります。
概念① 正則化バイアス：予測のための縮みが、係数への漏れになる
機械学習が予測に強い理由のひとつは、わざとバイアスを入れることです。Lassoは係数をゼロへ縮め、木は枝を刈り、ブースティングは早めに止まる。分散を減らして予測誤差の合計を下げる、理にかなった設計です。ところがこの縮みは、係数を因果として読んだ瞬間に毒になります。
今日の主役の1人目を「素朴プラグイン」と呼びます。処置 $w$ と共変量 $X$ をまとめて機械学習に渡して $y$ を予測させ、出てきた $w$ の係数（や $w$ を動かしたときの予測差）を効果として読む方法です。何がまずいのか。Lassoで言えば、罰則は交絡因子の係数も縮めます。縮んだ分の交絡は消えてなくなるわけではなく、モデルの説明し残しとして残る。そして交絡因子は定義上 $w$ と相関しているので、説明し残した交絡の受け皿に、罰則のかかっていない方向で最も相関の高い $w$ の係数がなるのです。交絡を消すために入れたはずのXの係数が縮んだツケを、θの推定値が肩代わりする。これが正則化バイアスです。
たちが悪いのは、このズレが偶然の誤差ではなく系統的なバイアスであることです。標本を増やしても、機械学習のニュイサンス推定はOLSよりゆっくりとしか真に近づかないため、バイアスは信頼区間の縮む速さより遅くしか消えません。つまり、標本が増えるほど「狭い区間で自信満々に間違える」方向へ進みます。手を動かす②の図1で、分布ごと真値から外れる様子を見ます。
概念② DML：残差×残差、FWLの機械学習版
処方箋はDay 16で伏線を張ったとおりです。部分線形モデル
$$
y = \theta w + g(X) + \varepsilon, \qquad w = m(X) + v
$$
を考えます。$g(X)$ が結果側、$m(X)$ が処置側のニュイサンス（迷惑関数）です。処置が2値なら $m(X)$ は傾向スコア $e(X)$（Day 7）そのもので、残差化は傾向スコアの親戚にあたります。今日は仕組みが見やすいように処置を連続量（投資の強度）にしますが、考え方は同じです。FWL定理（Day 16）は「重回帰の係数＝統制変数で残差化した残差同士の単回帰」でした。close_collegeでは educ の係数 0.0746 が残差×残差の単回帰と差 1.5e-16 で厳密に一致しましたね。DMLはこの手順の「残差を作る回帰」を、任意の機械学習に置き換えたものです。
なぜ残差×残差だと正則化バイアスが消えるのか。ここが今日いちばんの難所、Neyman直交性です。数式は原典に譲り、直感だけ言葉にします。
素朴プラグインでは、$\hat g$ の誤りがそのまま1次でθの推定誤差に乗ります。ニュイサンスを1割間違えれば、θもそれに比例して動く構造です。一方、残差×残差の推定量では、θの誤差に効くのは「$\hat m$ の誤り × $\hat\ell$ の誤り」の積だけになります。片方のニュイサンスを少しだけ間違えた状況を想像してください。その間違いは、もう片方の残差（平均ゼロで、間違いの方向と無相関）と掛け合わされて消えます。つまりどちらか一方を少し間違えても、1次では推定値が動かない。動くのは両方の誤りの積という2次の項だけで、それぞれが並の精度でも積は十分速く小さくなります。並の予測器2本を掛け算で1本の高精度な推定量に変える設計、と言ってもいいと思います。
ただしこの理屈には前提がひとつあります。残差を作るモデルが、残差化される当のデータを暗記していないことです。ここを守る仕掛けが交差適合で、手を動かす③で壊してから直します。
手を動かす①：素朴プラグインは分布ごと上振れる
まず高次元の交絡を合成します。共変量200列のうち5列だけが本物の交絡で、wにもyにも係数0.5で効きます。真の効果は θ=0.5。列がここまで増えると、実務ではまず正則化つきのモデルに手が伸びます。その手癖のまま処置の係数を読むと何が起きるかを測ります。
import numpy as np
from sklearn.linear_model import LassoCV
rng = np.random.default_rng(24)
N, P = 2000, 200
X = rng.normal(0, 1, (N, P)) # 共変量200列
beta = np.zeros(P)
beta[:5] = 0.5 # うち5列だけが交絡（wにもyにも効く）
w = X @ beta + rng.normal(0, 1, N) # 処置の強度（例：点検投資）
y = 0.5 * w + X @ beta + rng.normal(0, 1, N) # 真の効果 θ = 0.5
# (1) 何も調整しない単回帰
print("単回帰 y ~ w の傾き = %.3f" % np.polyfit(w, y, 1)[0])
# (2) 素朴プラグイン：Lasso に w と X をまとめて渡し、w の係数を読む
plugin = LassoCV(cv=3, random_state=0, alphas=30).fit(np.column_stack([w, X]), y)
print("素朴プラグインの θ = %.3f (真値 0.5)" % plugin.coef_[0])
print("交絡5列の係数 =", np.round(plugin.coef_[1:6], 2), "(真値 0.5)")
単回帰 y ~ w の傾き = 1.059
素朴プラグインの θ = 0.581 (真値 0.5)
交絡5列の係数 = [0.44 0.38 0.39 0.41 0.43] (真値 0.5)
単回帰の1.06は交絡まみれなので論外として、注目は素朴プラグインの 0.581 です。200列をちゃんとモデルに入れ、交差検証で罰則も選んだのに、真値0.5から系統的に上振れしています。3行目が種明かしで、交絡5列の係数は0.38〜0.44と、真値0.5からそれぞれ縮められています。縮んだ量は5列合計でおよそ0.45。交絡1列あたりの $w$ との共分散が0.5、$w$ の分散が2.25なので、ざっくり 0.45 × 0.5 / 2.25 ≒ 0.10 が $w$ の係数へ流れ込む勘定です。実際の上振れ+0.08とほぼ合います（θ自身にも罰則がかかって少し縮むぶん、勘定ぴったりにはなりません）。概念①の「縮んだツケの肩代わり」が、そのまま数字で見えました。
手を動かす②：残差×残差の手組みDML、FWLの絵が戻ってくる
同じデータをDMLで推定します。やることは3手順です。wをXから予測して残差を取る、yをXから予測して残差を取る、残差×残差で単回帰する。Day 16の手組みFWLの「回帰」を LassoCV に置き換えただけです。あわせて、後半の主題になる交差適合（ここでは2分割）も先に仕込んでおきます。半分のデータで学んだモデルで、残り半分を残差化する構造です。
from sklearn.model_selection import KFold
def dml_theta(X, w, y, seed=0):
"""手組みDML：2分割の交差適合で残差を作り、残差×残差の傾きを返す"""
rw, ry = np.zeros(len(w)), np.zeros(len(y))
for tr, te in KFold(2, shuffle=True, random_state=seed).split(X):
rw[te] = w[te] - LassoCV(cv=3, random_state=0, alphas=30).fit(X[tr], w[tr]).predict(X[te])
ry[te] = y[te] - LassoCV(cv=3, random_state=0, alphas=30).fit(X[tr], y[tr]).predict(X[te])
return (rw @ ry) / (rw @ rw)
print("手組みDMLの θ = %.3f (真値 0.5)" % dml_theta(X, w, y, seed=24))
手組みDMLの θ = 0.510 (真値 0.5)
素朴プラグインの0.581が 0.510 まで戻りました。使っている道具はまったく同じ LassoCV です。違うのは読み出し方の設計だけ。下の図の右が、200列の交絡を2次元に還元した残差×残差の散布図で、Day 16のFWLの絵がそのまま帰ってきています。
1回の比較では偶然かもしれないので、データ生成から100回繰り返して分布で見ます。
def make_data(seed):
rng = np.random.default_rng(seed)
X = rng.normal(0, 1, (N, P))
w = X @ beta + rng.normal(0, 1, N)
y = 0.5 * w + X @ beta + rng.normal(0, 1, N)
return X, w, y
naive_mc, dml_mc = [], []
for r in range(100): # 100回の反復（30秒ほどかかります）
Xr, wr, yr = make_data(r)
naive_mc.append(LassoCV(cv=3, random_state=0, alphas=30)
.fit(np.column_stack([wr, Xr]), yr).coef_[0])
dml_mc.append(dml_theta(Xr, wr, yr, seed=r))
print("素朴プラグイン: 平均 %.3f ± %.3f（100回中%d回が真値より上）"
% (np.mean(naive_mc), np.std(naive_mc), (np.array(naive_mc) > 0.5).sum()))
print("手組みDML : 平均 %.3f ± %.3f" % (np.mean(dml_mc), np.std(dml_mc)))
素朴プラグイン: 平均 0.563 ± 0.023（100回中100回が真値より上）
手組みDML : 平均 0.512 ± 0.023
ばらつき（±0.023）は両者ほぼ同じで、違いはバイアスだけです。素朴プラグインは平均0.563、真値から約2.7標準偏差ぶん離れたところに分布ごと引っ越していて、100回中100回が上振れ。この分布から計算した信頼区間は、狭い顔をして真値を外し続けます。手組みDMLの平均は 0.512 で、分布が真値0.5を囲みます。残る+0.012は概念②で言った「両方の誤りの積」の有限標本での残りで、標本を増やせば区間より速く縮んでいく性質のズレです。
手を動かす③：交差適合を外すと、残差は信号ごと消える
概念②の最後に置いた前提を、今度はわざと壊します。交差適合を省いて、全データで当てはめたモデルで、その同じデータを残差化したらどうなるか。極端な思考実験から入ると分かりやすいです。1近傍法（自分と同じ点を返すモデル）で残差を作れば、残差は全点で厳密にゼロ。残差×残差の回帰は0/0になり、何も推定できません。過学習はこの劣化版で、モデルが交絡だけでなくノイズまで暗記し、残差から本来残るべき信号 $v$ と $\varepsilon$ を削り取ってしまいます。
今度は非線形のニュイサンス（sinとtanh）を仕込み、学習器も勾配ブースティングに替えて確かめます。
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import cross_val_predict
rng = np.random.default_rng(1)
n = 1000
Z = rng.normal(0, 1, (n, 10)) # 共変量10列（今度は非線形に効く）
m = np.sin(2 * Z[:, 0]) + 1.5 * np.tanh(Z[:, 1]) # 処置側のニュイサンス m(x)
g = 2 * np.sin(2 * Z[:, 0]) + 3.0 * np.tanh(Z[:, 1]) # 結果側のニュイサンス g(x)
w2 = m + rng.normal(0, 1, n)
y2 = 0.5 * w2 + g + rng.normal(0, 1, n) # 真の効果 θ = 0.5
def gb():
return HistGradientBoostingRegressor(max_iter=500, random_state=0)
# (1) 交差適合なし：自分のデータで当てはめたモデルで、自分を残差化する
rw = w2 - gb().fit(Z, w2).predict(Z)
ry = y2 - gb().fit(Z, y2).predict(Z)
print("交差適合なし: θ = %.3f 処置残差の分散 = %.4f（本来は1.0）"
% ((rw @ ry) / (rw @ rw), rw.var()))
# (2) 5分割の交差適合：自分以外のfoldで学んだモデルで残差化する
kf = KFold(5, shuffle=True, random_state=0)
rw = w2 - cross_val_predict(gb(), Z, w2, cv=kf)
ry = y2 - cross_val_predict(gb(), Z, y2, cv=kf)
print("交差適合あり: θ = %.3f 処置残差の分散 = %.3f"
% ((rw @ ry) / (rw @ rw), rw.var()))
交差適合なし: θ = 1.523 処置残差の分散 = 0.0001（本来は1.0）
交差適合あり: θ = 0.536 処置残差の分散 = 1.291
交差適合なしでは、処置残差の分散が 0.0001 まで潰れています。設計上、残差には分散1.0のノイズ $v$ が残るはずなので、ブースティングがノイズをほぼ完全に暗記した、ということです。残った0.01%の屑どうしの比を取った推定値が1.523で、これはもう何も測っていません。実際、シードを変えて同じ推定を16回繰り返すと、交差適合なしは0.10から1.52まで暴れます（下の図の左）。たまたま0.5に近い値が出ることもあるのが、かえって怖いところです。交差適合ありは0.43〜0.59に収まり、真値の周りに素直に集まります。
分割の仕組みはこうです。fold Bの残差はfold Aで学んだモデルで作るので、どのモデルも残差化する相手のデータを見ていません。それでいて、全データが最終的に残差として使われるので、サンプルを捨てる無駄もありません。
最後に実務版です。econml の LinearDML は、今日手組みした「MLで2本残差化 → 残差×残差」をそのまま実装していて、cv 引数が交差適合の分割数、X=None は「効果は定数」の指定です（ここに共変量を渡すと効果の異質性のモデル化に進みます。Day 25の話です）。
from econml.dml import LinearDML
for cv in [1, 5]:
est = LinearDML(model_y=gb(), model_t=gb(), cv=cv, random_state=0)
est.fit(y2, w2, X=None, W=Z)
lo, hi = est.effect_interval()
print("LinearDML cv=%d: θ = %.3f 95%%CI [%.3f, %.3f]"
% (cv, est.effect()[0], lo[0], hi[0]))
LinearDML cv=1: θ = 1.523 95%CI [0.538, 2.509]
LinearDML cv=5: θ = 0.536 95%CI [0.466, 0.607]
cv=1（交差適合なし）は手組みの1.523をそのまま再現し、区間も[0.538, 2.509]と間延びした上に真値を含みません。cv=5は 0.536、95%CI [0.466, 0.607] で真値0.5を狭く囲みます。DMLの漸近理論（普通の $\sqrt{n}$ の信頼区間が使えること）は、直交性と交差適合の2つセットで初めて成立します。片方でも欠くと、点推定だけでなく区間も信用できなくなる、というのが今日の締めの教訓です。
つまづき・誤解しやすい点
-
「予測精度が高いモデルほど、因果推定も正確」ではありません。 手を動かす①と②で、ニュイサンスの予測に使った道具は同じ
LassoCVでした。差を生んだのは精度ではなく係数の読み出し方の設計です。正則化バイアスはR²やRMSEにいっさい顔を出さないので、予測指標でモデルを磨いても検出できません。ニュイサンスのモデル選択は予測基準でよく、読み出しは直交化した推定式で行う。役割分担で覚えるのが安全です。 - DMLは交絡を自動で見つけてはくれません。 消してくれるのは「共変量に入れた変数を経由するバックドア」だけで、非交絡性が仮定であることはDay 6から1ミリも変わっていません。何を入れるかはバックドア基準（Day 11）で決め、合流点や処置後の変数を200列に混ぜれば、DMLは高精度にバイアスを推定するだけです。隠れ交絡への耐性は感度分析（Day 9）で別途語ります。
- DMLのθは「効果は1つ」を仮定した平均です。 部分線形モデルは全員に同じθが効く書き方で、効果が人によって違うとき、θ̂は処置残差の分散で重みづけた加重平均になります。Day 16で見た「回帰＝分散加重平均」の親戚です。誰に効くかが論点なら、メタラーナー（Day 23）やCausal Forest（Day 25）のCATE系へ進みます。
GISデータ実務での使い方
- 共変量が数百列ある施策評価は、DMLの主戦場です。 管路の更新投資や点検頻度の効果を測ろうとすると、管路属性・土壌・気象・衛星由来の指標で説明変数はすぐ数百列になります。全列OLSは分散が暴れ、変数選択は恣意性を疑われる。DMLなら「処置を予測するモデル」と「結果を予測するモデル」を2本立てるだけで、変数選択は正則化に任せられます。
-
手持ちの予測パイプラインがそのままニュイサンスになります。 漏水リスクや劣化予測で運用しているLightGBMのような資産を、
model_yとmodel_tに差し替えるだけでDMLに転用できます。ただし交差適合は必ずセットで。econmlはcvを内蔵しているので、既定のまま使えば手を動かす③の事故は起きません。 - 空間データではfoldの切り方が本丸になります。 隣接する管路や自治体は互いに似ているので、ランダムなKFoldでは訓練foldの情報が検証foldへ漏れ、「交差適合したつもりで暗記している」状態が起きます。空間ブロックで切る交差適合が処方箋で、これは空間機械学習の空間CVと同じ話です。空間交絡・空間干渉ごとまとめてDay 28で掘ります。
参考（本棚）
- Chernozhukov et al. (2018) Double/Debiased Machine Learning for Treatment and Structural Parameters — DMLの原典。正則化バイアス・Neyman直交性・交差適合の3点セットはこの論文の骨格
- Robert Osazuwa Ness『Causal AI』 — 機械学習と因果推論の接続の見取り図。今日の下敷き
- Matheus Facure『Causal Inference in Python』 — 残差化（直交化）とDebiasingの実装視点
-
econml ドキュメント（PyWhy） —
LinearDMLの引数と推定量の仕様 - 図はすべて合成データです。真の効果 θ=0.5 と交絡の構造を私たちが握ることで、素朴プラグインの上振れと交差適合の有無を安全に再現しました。生成スクリプトは元ノート側に添付しています
次回予告（Day 25）
今日のDMLは、効果を1つの定数θとして推定しました。でもDay 22で見たとおり、実務の問いはしばしば「平均」ではなく「誰に効くか」です。残差化した世界で木を育て、効果の異質性 τ(x) を信頼区間つきで切り出すのが次回の Causal Forest（econml CausalForestDML）です。分割の規準を予測誤差から「効果の差」に替えるところ、そして信頼区間の生命線になる honesty（標本分割）が今日の交差適合と同じ精神であるところまで、手を動かして確かめます。