因果推断：双重机器学习-ddml

发布时间：2023-05-17 阅读 735

Stata连享会主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：李金桐 (中山大学)
邮箱：lijt65@mail2.sysu.edu.cn

1. 理论基础

1.1 为什么我们需要 DDML

实证研究往往会面临一个质疑：模型设定是正确的吗？例如，研究在班级中加入助教对教学质量的影响，常见的方法是构造回归方程：

S c o r e = β_{0} + β_{1} A s s i s t a n t + \sum_{j} β_{j} \times C o n t r o l + ε

其中， $S c o r e$ 代表成绩； $A s s i s t a n t$ 是标志是否加入助教的二元变量； $C o n t r o l$ 为控制变量，可能包括：每天学习时间、作业完成率、出勤率等。那么这些特征的关系真的是线性的吗？显然不是。例如随着学习时间增加，成绩自然会提高，然而学习时间过长很可能导致疲惫、睡眠不足等，进而造成学习效率下降，反而使得成绩下降。

需要注意的是，我们实际上并不关心学习时间对成绩的影响，我们只希望研究 $β_{1}$ ，我们只是需要处理控制变量对 $β_{1}$ 造成的影响。

接下来用更严谨的方法描述上述问题。考虑因果模型：

Y = D θ_{0} + g_{0} (X) + U, E [U ∣ X, D] = 0

D = m_{0} (X) + V, E [V ∣ X] = 0

其中 $Y$ 是模型的 Outcome， $D$ 是因果模型的 treatment。这里，我们关注 $θ_{0}$ ，即 treatment 的因果效应。一种常见的思路是，通过假设 (例如常见的线性假设)，或者利用一定方法 (通常是机器学习) 估计，得到 ${\hat{g}}_{0}$ ，随后就可以利用线性回归得到 ${\hat{θ}}_{0}$ ：

{\hat{θ}}_{0} = \frac{cov (D, Y - {\hat{g}}_{0} (X))}{var (D)} = \frac{\frac{1}{n} \sum_{i \in I} D_{i} (Y_{i} - {\hat{g}}_{0} (X_{i}))}{\frac{1}{n} \sum_{i \in I} D_{i}^{2}} (1)

接下来，很自然的想要研究这个估计量是否无偏。遗憾的是 ${\hat{θ}}_{0}$ 往往是有偏的：

\begin{aligned} \sqrt{n} ({\hat{θ}}_{0} - θ_{0}) & = \sqrt{n} \frac{\frac{1}{n} \sum_{i \in I} D_{i} (Y_{i} - {\hat{g}}_{0} (X_{i}))}{\frac{1}{n} \sum_{i \in I} D_{i}^{2}} \\ - (\sqrt{n} \frac{\frac{1}{n} \sum_{i \in I} D_{i} (Y_{i} - g_{0} (X_{i}))}{\frac{1}{n} \sum_{i \in I} D_{i}^{2}} - \sqrt{n} \frac{\frac{1}{n} \sum_{i \in I} D_{i} U_{i}}{\frac{1}{n} \sum_{i \in I} D_{i}^{2}}) \\ = \underset{:= a}{\underset{⏟}{{(\frac{1}{n} \sum_{i \in I} D_{i}^{2})}^{- 1} \frac{1}{\sqrt{n}} \sum_{i \in I} D_{i} U_{i}}} + \underset{:= b}{\underset{⏟}{{(\frac{1}{n} \sum_{i \in I} D_{i}^{2})}^{- 1} \frac{1}{\sqrt{n}} \sum_{i \in I} D_{i} (g_{0} (X_{i}) - {\hat{g}}_{0} (X_{i}))}} \end{aligned}

可以看出误差分为两项。 $a$ 项来自于 $U$ 和 $D$ 的独立性，即 $\frac{cov (D, U)}{var (D)}$ ，若二者不独立则会造成偏误。然而问题来源于 $b$ 项，我们将其展开为以下形式：

b = {(E [D_{i}^{2}])}^{- 1} \frac{1}{\sqrt{n}} \sum_{i \in I} m_{0} (X_{i}) (g_{0} (X_{i}) - {\hat{g}}_{0} (X_{i})) + o_{P} (1)

注意到 $m_{0} (X_{i}) (g_{0} (X_{i}) - {\hat{g}}_{0} (X_{i}))$ 项。首先， $g_{0}$ 的估计往往存在误差，例如对于高维数据，往往会采用正则项处理，造成正则化误差，此时 $b$ 项发散。此外， $m_{0} (X_{i})$ 是数据本身的性质，因此数据会决定偏误的大小而无法改变，导致估计非常不稳健。

综合以上推论，可以说因果模型 treatment effect 的传统估计方法并不完美。因此，我们引入Double/Debiased Machine Learning (DDML) 的概念，为因果估计提供更为稳健的方法。

1.2 DDML 的本质

根据上述推论，传统因果估计方法不稳健的核心在于 $m_{0} (X_{i}) (g_{0} (X_{i}) - {\hat{g}}_{0} (X_{i}))$ 项，其中 $g_{0} (X_{i}) - {\hat{g}}_{0} (X_{i})$ 部分是 $g_{0}$ 的估计误差，显然是难以避免的。因此更为实际的考虑是消除 $m_{0} (X_{i})$ 。其出现原因在于用于回归的 $D$ 实际上包含了 $X$ 的信息。

如图，注意到 $V$ 实际上可以看作工具变量，此时可以构造估计：

\begin{aligned} θ_{0} & = \frac{cov (V, Y - g_{0} (X))}{cov (V, D)} \\ = \frac{cov (D - m_{0} (X), Y - g_{0} (X))}{cov (D - m_{0} (X), D)} \\ = \frac{cov (D - m_{0} (X), D θ_{0} + U)}{cov (D - m_{0} (X), D)} \end{aligned}

为了求 $V$ ，可以采用：

\hat{V} = D - {\hat{m}}_{0} (X)

其中 ${\hat{m}}_{0} (X)$ 可以通过 $X$ 对 $D$ 回归得到，因此我们得到一种新的估计：

{\overset{˘}{θ}}_{0} = \frac{\frac{1}{n} \sum_{i \in I} {\hat{V}}_{i} (Y_{i} - {\hat{g}}_{0} (X_{i}))}{\frac{1}{n} \sum_{i \in I} \hat{V_{i}} D_{i}} (2)

在这个估计下，新的 $b$ 项变为：

b^{*} = {(E [D_{i}^{2}])}^{- 1} \frac{1}{\sqrt{n}} \sum_{i \in I} ({\hat{m}}_{0} (X_{i}) - m_{0} (X_{i})) (g_{0} (X_{i}) - {\hat{g}}_{0} (X_{i}))

此时偏误仅仅取决于回归误差，因此这个估计更为稳健。直观角度上，线性回归是拟合 $Y$ 在特征空间 $X$ 的最佳投影 (既误差最小化)，所以残差垂直于样本空间 $X$ ，最大限度消除了 $X$ 的相关性。

此时我们已经在因果模型中构造了稳健的估计量。接下来，我们希望从中总结出一套普适性更强的方法论，以用于更多场景。我们对上述过程进行总结：

存在目标函数 $θ_{0}$ 和其他不太关心的回归函数 $g_{0}$ 和 $m_{0}$ ；
$θ_{0}$ 的传统估计方法中，如果拟合的函数 ${\hat{g}}_{0}$ 存在误差，则误差项不收敛；
需要构造估计量使得 $g_{0}$ 和 $m_{0}$ 出现误差时，偏误依然足够小。

这里我们用 Gateaux Dervative 描述 $g_{0}$ 估计误差对偏误的影响：

\partial_{f} E [φ (W; θ_{0}, f_{0})] = lim_{r \to 0^{+}} \frac{E [φ (W; θ_{0}, f_{0} + r (f - f_{0}))] - E [φ (W; θ_{0}, f_{0})]}{r}

我们将 (1) 式的传统估计方法重新写成：

\frac{1}{n} \sum_{i \in I} φ (W; {\hat{θ}}_{0}, {\hat{g}}_{0}) = 0

其问题在于 Gateaux Dervative 不为零：

\partial_{g} E [φ (W; θ_{0}, g_{0})] [g - g_{0}] \neq 0

换言之， $g_{0}$ 的微小扰动会导致 $ϕ$ 发生较大的变化，因此估计不准时存在较大误差。

我们将 (2) 式的估计方法重新写成：

\frac{1}{n} \sum_{i \in I} ψ (W; {\overset{˘}{θ}}_{0}, {\hat{η}}_{0}) = 0

其中 ${\hat{η}}_{0} = ({\hat{m}}_{0}, {\hat{g}}_{0})$ 称为 Nuisance Parameter。可以证明这个估计量的 Gateaux Dervative 为零：

\partial_{η} E [φ (W; θ_{0}, η_{0})] [η - η_{0}] = 0

直观来说，Gateaux Dervative 事实上刻画了 $ϕ$ 对 $g_{0}$ 变化的敏感程度。如果 Gateaux Dervative 等于零，意味着这个估计量更为稳健。我们称这个性质为 Neyman Orthogonality。这种构造就是 Double/Debiased Machine Learning 的思想基础。

此外，估计无偏重要的一步是 Cross-fitting，用来降低 overfitting 带来的估计偏差。将总样本分成样本 1，样本 2 两份。先用样本 1 估计残差，样本 2 估计 $θ_{1}$ ，再用样本 2 估计残差，样本 1 估计 $θ_{2}$ ，取平均得到最终的估计 $θ = (θ_{1} + θ_{2}) / 2$ 。当然也可以进一步使用 K-Fold 来增加估计的稳健性。

1.3 如何构造 DDML

接下来，我们希望找到通用的方法构造 $ψ$ 满足Neyman Orthogonality。不妨设 $max_{θ, β} E_{p} [l (W; θ, β)]$ ，我们的目标可以通过最大化某个函数得到，通常为对数似然，故在极值处导数为 0：

E_{p} [\partial_{θ} l (W; θ_{0}, β_{0})] = 0, E_{p} [\partial_{β} l (W; θ_{0}, β_{0})] = 0

因此，一个自然的想法是令 $ϕ (W; θ, β) = \partial_{θ} l (W; θ, β)$ 。然而，为了满足 Neyman Orthogonality，我们还需要自己构造估计量：

ψ (W; θ, β) = \partial_{θ} l (W; θ, β) - μ \partial_{β} l (W; θ, β)

此时只需要找到合适的 $μ$ 使 $ψ$ 满足 Neyman Orthogonality：

\partial_{β} ψ (W; θ, β) = \partial_{θ β} l (W; θ, β) - μ \partial_{β β} l (W; θ, β) = 0

可以得到解析解： $μ = J_{θ β} J_{β β}^{- 1}$ 。

2. Stata 实现

2.1 命令介绍

ddml 是用于 Double Debiased Machine Learning 的 Stata 包，支持五种不同的模型，允许二元或连续 treatment 变量，支持内生性、高维控制，并支持工具变量，且支持多种不同的机器学习程序。

. ssc install ddml, replace

利用 ddml 包进行 Double Debiased Machine Learning 分为四步进行：

初始化并选择模型：

ddml init model [if] [in] [ , mname(name) kfolds(integer) 
    fcluster(varname) foldvar(varlist) reps(integer) 
    norandom tabfold vars(varlist) ]

添加监督机器学习以估计条件期望：

ddml eq [ , mname(name) vname(varname) learner(varname) 
    vtype(string) predopt(string) ] :  command depvar vars [ , cmdopt ]

Cross-fitting：

ddml crossfit [ , mname(name) shortstack ]

估计因果效应：

ddml estimate [ , mname(name) robust cluster(varname) vce(type) att trim(real) ]

2.2 模型设定

Partial linear model [partial]

Y = D θ_{0} + g_{0} (X) + ζ, E (ζ ∣ D, X) = 0

D = m_{0} (X) + V, E (V ∣ X) = 0

其中目标是在控制 $X$ 的同时估计 $θ_{0}$ 。为此，我们使用有监督的机器学习估计条件期望 $E [Y | X]$ 和E $[D | X]$ 。

Partial linear IV model [iv]

Y - D θ_{0} = g_{0} (X) + ζ, E (ζ ∣ Z, X) = 0

Z = m_{0} (X) + V, E (V ∣ X) = 0

在这个模型中， $Y$ 和 $D$ 之间存在结构或因果关系，其关系由 $θ$ 来表示。 $g_{0} (X) + ζ$ 可以看作是随机误差，其中 $g_{0} (X)$ 可以被协变量 $X$ 所解释。由于 $Y$ 和 $D$ 是共同确定的，所以我们引入一个外部因素，称为工具 $Z$ ，以解释 $D$ 的外生变化。

注意， $Z$ 应该影响 $D$ 。这里的 $X$ 再次作为 confounding factors，因此我们可以认为 $Z$ 中的变化是外生的，仅以 $X$ 为条件。

一个来自生物统计学的例子： $Y$ 为健康指标， $D$ 为吸烟的指标，则 $θ_{0}$ 代表着抽烟对健康的影响，健康指标 $Y$ 和吸烟行为 $D$ 同时被确定。 $X$ 代表患者特征， $Z$ 是医生建议不要吸烟 (或其他行为治疗)，这一行为由 $X$ 决定， $Z$ 只能通过改变行为 $D$ 去影响结果 $Y$ 。

Interactive model [interactive]

Y = g_{0} (D, X) + U, E (U ∣ X, D) = 0

D = m_{0} (X) + V, E (V ∣ X) = 0

在该模型设定中，treatment effect 完全异质且 $D$ 为二元变量。由于 $D$ 与 $X$ 是不可线性分离的，所以在 $D$ 为二元变量的情况下，该模型比部分线性模型的应用广泛性更好。

我们在模型中感兴趣的目标参数是平均 treatment effect(ATE)：

θ_{0} = E [g_{0} (1, X) - g_{0} (0, X)]

另一个常见的目标参数是接受 treat 者的平均 treatment effect(ATTE)：

θ_{0} = E [g_{0} (1, X) - g_{0} (0, X) ∣ D = 1]

Interactive IV model [interactiveiv]

Y = ℓ_{0} (D, X) + ζ, E (ζ ∣ Z, X) = 0