双重机器学习（DML）：原理、代码和几个坑

双重机器学习（double/debiased machine learning, DML）这几年在国内期刊里越来越常见，多半出现在稳健性检验里，用来回应“控制变量的函数形式设定错了怎么办”。这篇记下它解决什么、不解决什么。

要解决的问题

部分线性模型（partially linear regression, PLR）：

$\theta$ 是关心的处理效应，$g_0$ 和 $m_0$ 是未知函数，可能是非线性的，$X$ 的维度也可能很高。OLS 把 $X$ 线性地放进回归，等于假定 $g_0$ 是线性的。设定错了，$X$ 里没控制干净的部分同时影响 $D$ 和 $Y$，$\hat\theta$ 就有偏。

注意 PLR 假定处理效应是一个常数。如果效应因人而异，PLR 估出的 $\theta$ 是按 $D$ 的条件方差加权的平均，不等于平均处理效应（ATE）。二元处理想要 ATE，应该用交互模型（ddml 里是 interactive，DoubleML 里是 DoubleMLIRM）。

那直接用机器学习估计 $g_0$，再拿 $Y-\hat g(X)$ 对 $D$ 回归行不行？不行。机器学习靠正则化换预测精度，$\hat g$ 的偏差会直接传到 $\hat\theta$ 里，收敛速度也到不了 $\sqrt n$。Chernozhukov et al.（2018）一开头就用模拟展示了这种“朴素”估计量的分布是怎么偏掉的。

两个关键：正交化和交叉拟合

正交化（partialling out）。 记 $\ell_0(X)=\E[Y\mid X]$，两式相减得到

所以只要分别用机器学习预测 $Y$ 和 $D$，取残差 $\hat W=Y-\hat\ell(X)$、$\hat V=D-\hat m(X)$，再做残差对残差的回归：

这和 Frisch–Waugh–Lovell 定理、Robinson（1988）的半参数估计是同一个思路，只是把条件期望换成了机器学习的预测。它的好处叫 Neyman 正交性：估计方程对 $\ell$ 和 $m$ 的一阶导数为零，两个预测的误差只以乘积的形式进入 $\hat\theta$ 的偏差。所以每个预测只要以快于 $n^{-1/4}$ 的速度收敛就够了，很多机器学习方法在一定条件下能做到。

交叉拟合（cross-fitting）。 如果用同一批样本既训练模型又算残差，过拟合会让残差系统性地偏小，带来新的偏差。做法是把样本分成 $K$ 份，每一份的残差都用其余 $K-1$ 份训练出的模型来算，再把所有残差合在一起估计 $\theta$。切分有随机性，所以一般重复切分几次，取估计值的中位数。

只切分一次时，标准误就是残差回归的异方差稳健标准误；重复切分取中位数时，还会把不同切分之间的差异计入标准误。

模拟：非线性混杂

set seed 20260929
set obs 2000
forvalues j = 1/10 {
    gen x`j' = rnormal()
}
gen d = 0.8*sin(x1) + 0.6*(x2^2 - 1) + 0.5*(x3 > 0) + rnormal()
gen y = 0.5*d + 1.2*(x2^2 - 1) + 0.8*cos(x1) + 0.5*x3*x4 + rnormal()

真实的 $\theta=0.5$。$x_2^2$ 同时进入 $D$ 和 $Y$，把 $x_1$ 到 $x_{10}$ 线性放进回归是吸收不了它的。

Stata 的结果（ddml，5 折，重复切分 3 次）：

堆叠这一行是按下文“代码”一节原样跑出来的；三次切分的估计分别是 0.526、0.504、0.502，取中位数。换一个种子，结果会在 0.50 上下小幅变动。

同样的数据生成过程在 Python 里再跑一遍（DoubleML 0.11.4，5 折，重复 3 次）。Python 和 Stata 的随机数不同，所以这是另一份数据：

几点体会：

1. 学习器是线性的，DML 和 OLS 差不多。 用线性回归或只放原始变量的 Lasso 做学习器，结果和 OLS 几乎一样（两者不是严格相等，交叉拟合和正则化还是会带来一点差别）。DML 能不能纠偏，取决于学习器能不能逼近 $g_0$ 和 $m_0$。
2. 不同学习器的结果可能差不少。 这次模拟里随机森林的置信区间连真实值 0.5 都没覆盖到。这只是一次模拟，不能说明随机森林一定不行，但说明学习器的选择要当回事：学习器组合最好事先定好，比如像 pystacked 那样事先给定一组学习器做堆叠，由交叉验证决定权重；换学习器的结果放在稳健性检验里报告，而不是试了很多个只挑最好看的那个。这个例子里堆叠权重大部分给了梯度提升（$D$ 方程 0.77，$Y$ 方程 0.90），其余几乎都给了随机森林，OLS 和 Lasso 的权重接近 0。
3. DML 解决的是函数形式，不是内生性。 它的识别仍然靠“给定 $X$ 后处理近似随机”（条件独立，selection on observables），此外还要求给定 $X$ 后 $D$ 仍有足够的变化（重叠），以及 $g_0$、$m_0$ 能被估得足够好。遗漏了观测不到的混杂，换多复杂的学习器也没用；把受处理影响的变量放进 $X$（坏控制），同样会出问题。

代码

Stata（ddml，Ahrens et al., 2024）：

* 需要 Stata 16 以上，以及装了 scikit-learn 的 Python（pystacked 要用）
ssc install ddml
ssc install pystacked

global X x1-x10
set seed 2026
ddml init partial, kfolds(5) reps(3)
ddml E[Y|X]: pystacked y $X, type(reg) method(ols lassocv rf gradboost)
ddml E[D|X]: pystacked d $X, type(reg) method(ols lassocv rf gradboost)
ddml crossfit
ddml estimate, robust
ddml extract, show(stweights)     // 看各学习器的堆叠权重

Stata 找不到装了 scikit-learn 的 Python 时，用 set python_exec 指定 python.exe；Python 包装在虚拟环境里的，再用 set python_userpath 指向虚拟环境的 site-packages 目录。

Python（DoubleML，Bach et al., 2022）：

import numpy as np
import doubleml as dml
from sklearn.ensemble import HistGradientBoostingRegressor

# df 是 pandas.DataFrame，含 y、d 和 x1…x10
x_cols = [f"x{j}" for j in range(1, 11)]
data = dml.DoubleMLData(df, y_col="y", d_cols="d", x_cols=x_cols)
ml = HistGradientBoostingRegressor(random_state=1)

np.random.seed(42)                 # 样本切分用的是 numpy 的全局随机数
plr = dml.DoubleMLPLR(data, ml_l=ml, ml_m=ml, n_folds=5, n_rep=3)
plr.fit()
print(plr.summary)

ml_l 预测的是 $\E[Y\mid X]$，ml_m 预测的是 $\E[D\mid X]$。

面板数据

国内论文多数用企业—年份面板，直接套 PLR 有两个问题。

固定效应。 PLR 本身不处理个体固定效应。常见的做法是把个体、年份虚拟变量放进 $X$，或者先对变量做组内去均值再做 DML。这两种都只是权宜之计：前者对树模型不太友好；后者和非线性的 $g_0$ 放在一起，去均值之后未必还是原来那个模型。Clarke and Polselli（2026）专门讨论了带固定效应的静态面板怎么做 DML，要用的话值得先读。

聚类。 同一家企业的观测不应该被切到不同的折里，否则训练集和验证集之间并不独立。ddml 可以按企业切分，并报告聚类标准误：

ddml init partial, kfolds(5) reps(3) fcluster(id)
* 指定学习器、crossfit 的步骤同上
ddml estimate, cluster(id)

参考文献