理解Double/debiased machine learning

基础：线性回归

考虑一个经典线性高斯模型：

$$y=ax+u$$

其中U服从标准高斯分布，a是回归系数，那么回归的目的是找到一个a，使得x与u独立，即

$$\begin{array}{rl}& cov(y-ax,x)=0\\ ⟹& cov(y,x)-acov(x,x)=0\\ ⟹& a=cov(y,x)/cov(x,x)\end{array}$$

这些知识在后面会反复用到。

Double/debiased machine learning

我们首先考虑这个因果模型(Partially Linear Regression)：

$$\begin{array}{cl}Y=D{\theta }_0+g_0(X)+U,& \mathrm{E}[U\mid X,D]=0,\\ D=m_0(X)+V,& \mathrm{E}[V\mid X]=0,\end{array}$$

其中$Y$ 是outcome,$D$ 是policy/treatment,且X是个高维的变量

$$X=(X_1,\dots ,X_p)$$

在这里，我们最关心的是${\theta }_0$，因为当我们给定某个X的时候，${\theta }_0$就表示了在这个X的群体中，D这个treatment的因果效应。那么估计${\theta }_0$最简单的方法就是，用机器学习去估计，这里我们先将数据随机分成两份，分别是$I,I^c$，不妨假设${\hat{g}}_0$是通过ML估计的函数，于是给定${\hat{g}}_0$，${\hat{\theta }}_0$可以用线性回归得到
$${\hat{\theta }}_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}\text{\hspace{1em}(1)}$$
这样的估计量其实很容易有bias，主要原因的它非常依赖于${\hat{g}}_0$的准确度，万一它有一点偏差就会产生很大的影响，我们来分析一下：

$$\begin{array}{rl}\sqrt{n}({\hat{\theta }}_0-{\theta }_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}-\left(\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}\right)\\ & =\underset{:=a}{\underbrace{{\left(\frac{1}{n}\sum _{i\in I}{D}_i^2\right)}^{-1}\frac{1}{\sqrt{n}}\sum _{i\in I}{D}_i{U}_i}}+\underset{:=b}{\underbrace{{\left(\frac{1}{n}\sum _{i\in I}{D}_i^2\right)}^{-1}\frac{1}{\sqrt{n}}\sum _{i\in I}{D}_i(g_0(X_i)-{\hat{g}}_0(X_i))}}.\end{array}$$

误差分为两部分，一部分是U和D的独立性，如果不独立的话会形成误差，即

$${\theta }_0=\frac{cov(D,Y-g_0(X))}{var(D)}-\underset{\ne 0}{\underbrace{\frac{cov(D,U)}{var(D)}}}$$

当两者独立的话，第一项的收敛其实很快，服从$a\sim N(0,\overline{\Sigma })$. 问题在$b$ 项，如果X维度很高，那么我们可能会在拟合g时加点正则项（如lasso），就会形成正则化误差，使得g有错误的估计，于是导致b项发散，

$$b={\left(E\left[D_i^2\right]\right)}^{-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)$$

b项可以进一步写为上式，可以看到主要问题出在$m_0(X_i)(g_0(X_i)-{\hat{g}}_0(X_i))$这个地方，$m_0(X_i)$的大小将会极大的影响误差的大小，而$m_0$是数据的性质，无法改变，就会导致我们的估计很不稳定（数据会决定bias的大小）。

有没有可能消掉这一项呢？我们先分析下$m_0(X_i)$出现的最主要的原因是我们用D来回归，然而D包含了X的信息，于是才会有$m_0(X_i)$这一项出现。那有没有可能不用D来回归呢？一个重要的观察是，V其实可以看做工具变量，

于是D和Y之前的因果效应可以被以下式子估计：

$${\theta }_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{\theta }_0+U)}{cov(D-m_0(X),D)}$$

因此，为了求V，我们可以求：

$$\hat{V}=D-{\hat{m}}_0(X)$$

这个${\hat{m}}_0(X)$可以用X对D回归来得到，于是得到了新的一种估计
$${\breve{\theta }}_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}\widehat{V_i}{D}_i}\text{\hspace{1em}(2)}$$
直观来看，因为V不包含X的信息，所以在上面b项中$m_0(X_i)$带来的误差可以被消去，事实上，在这个新的估计量下，这个新的b项将变为

$$b^{\ast }={\left(E\left[D_i^2\right]\right)}^{-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))$$

显然，这将能得到一个更为robust估计，因为现在bias只取决于回归误差了。

PS：从${\theta }_0=\frac{cov(D-m_0(X),D{\theta }_0+U)}{cov(D-m_0(X),D)}$这个式子来看，其实我们完全可以将下面的D也换成${\theta }_0=\frac{cov(D-m_0(X),D{\theta }_0+U)}{cov(D-m_0(X),D-m_0(X))}$，这并不影响我们的结果，这样的估计量最早由Robinson (1988) [1] 提出.

Neyman orthogonality and moment conditions.

上文只是介绍了一个Partially Linear Regression构造robust估计量的方法，那这套方法能不能被抽象出来，使其适用于更多的场景？

那我们先总结一些性质，

1. 我们有一个目标参数${\theta }_0$，然后有其他不太关心的但可以用各种机器学习方法得到的回归函数，$g,m$
2. 估计${\theta }_0$有多种不同的方法，有的好有的坏，而坏的原因在于某个用机器学习拟合的$\hat{g}$如果出现误差，则误差项的均值将不再为0，从而不收敛。
3. 为此，我们需要构造一种“正交”的估计量，使得即使g或者m出现误差也能保证误差足够小。

针对第一个性质，我们可以将公式1的估计量写成：

$$\frac{1}{n}\sum _{i\in I}\phi (W;{\hat{\theta }}_0,{\hat{g}}_0)=0$$

在上面的例子，就是线性回归中，找到噪声与x协方差独立，即

$$\begin{gathered} cov(Y-\theta D-g(X),D)=0 \\ ⟹\theta =cov(Y-g(X),D)/var(D) \end{gathered}$$

因此
$$\phi (W;\theta ,g)=(Y-\theta D-g(X))D$$

相当于$cov(Y-\theta D-g(X),D)=0$.
然而在这个估计量的问题是，当g不太准的时候，将会有较大的误差，这个误差理解为g一些微小的扰动而导致$\phi$发生较大的变化。那么怎么衡量这一变化呢？其实就是Gateaux dervative导数：
$${\partial }_{g}E[\phi (W;{\theta }_0,g_0)][g-g_0]\ne 0$$
不等于0。这个Gateaux dervative的定义是

$${\partial }_{g}E[\phi (W;{\theta }_0,g_0)]=\underset{r\to 0^{+}}{\lim }\frac{E[\phi (W;{\theta }_0,g_0+r(g-g_0))]-E[\phi (W;{\theta }_0,g_0)]}{r}$$

直观上，它刻画了$\phi$对$g$变化的敏感度，那如果这个导数等于0，则意味着我们找到了一个更为robust的估计量，事实上，我们上文中的第二种方法正是这样的估计量，公式2可以写成：

$$\frac{1}{n}\sum _{i\in I}\psi (W;{\breve{\theta }}_0,{\hat{\eta }}_0)=0$$

其中${\hat{\eta }}_0=({\hat{m}}_0,{\hat{g}}_0)$，被称为nuisance parameter，并且

$$\psi (W;\theta ,\eta )=(Y-\theta D-g(x))(D-m(x))$$

相当于$cov(Y-\theta D-g(x),D-m(x))$。可以证明这个估计量的Gateaux dervative为0：
$${\partial }_{\eta }E[\phi (W;{\theta }_0,{\eta }_0)][\eta -{\eta }_0]=0\text{\hspace{1em}(3)}$$
我们把这个性质称为Neyman orthogonality.

那有没有一种通用的方法来构造满足Neyman orthogonality性质的这样一个$\psi$?

不妨设

$$\underset{\theta ,\beta }{\max }{E}_p[l(W;\theta ,\beta )]$$

我们的目标可以通过最大化某个函数得到（通常为对数似然度），于是在极值下的导数为0：

$$E_p[{\partial }_{\theta }l(W;{\theta }_0,{\beta }_0)]=0,E_p[{\partial }_{\beta }l(W;{\theta }_0,{\beta }_0)]=0$$

因此，自然地我们可以令$\varphi (W;\theta ,\beta )={\partial }_{\theta }l(W;\theta ,\beta )$，然而这样的函数不一定满足Neyman orthogonality，所以我们可以自己构造一个估计量

$$\psi (W;\theta ,\beta )={\partial }_{\theta }l(W;\theta ,\beta )-\mu {\partial }_{\beta }l(W;\theta ,\beta )$$

只要能找到一个合适的$\mu$，使得$\psi$，满足Neyman orthogonality就可以了，即

$${\partial }_{\beta }\psi (W;\theta ,\beta )={\partial }_{\theta \beta }l(W;\theta ,\beta )-\mu {\partial }_{\beta \beta }l(W;\theta ,\beta )=0$$

我们可以得出解析解，$\mu =J_{\theta \beta }{J}_{\beta \beta }^{-1}$

参考资料

[1] Robinson, Peter M. “Root-N-consistent semiparametric regression.” Econometrica: Journal of the Econometric Society (1988): 931-954.
[2] Double/debiased machine learning for treatment and structural parameters
[3] https://observablehq.com/@herbps10/one-step-estimators-and-pathwise-derivatives
[4] Visually Communicating and Teaching Intuition for Influence Functions
[5] Approximation Theorems of Mathematical Statistics.
[6] https://vanderlaan-lab.org/2019/12/24/cv-tmle-and-double-machine-learning/
[7] https://academic.oup.com/biostatistics/article/21/2/353/5631845