Double Machine Learning

1.从线性回归说起

从观测数据获得因果效应的一个简单方式是使用线性回归，控制confounders的影响：
$$Sale{s}_i=\alpha +\tau Pirc{e}_i+{\beta }_{1}tem{p}_i+{\beta }_{2}Cos{t}_i+{\beta }_{3}Weekda{y}_i+e_i$$ $\tau$是我们唯一需要关注的，因为$\tau$是需要求的价格对销量的因果效应，其他变量我们并不关心，但是在线性模型中需要正确对待，否则得到的$\tau$是错误估计。比如，如果温度对销量的影响如果是二次的，根据上式估计出的treatment effect就是错的。

2.Frisch-Waugh-Lovell Law

假设有一系列特征$X_1$和$X_2$，

普通线性回归：$\hat{Y}={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2$

基于Frisch-Waugh-Lovell Law回归：

• 用$X_2$回归$Y$：${\hat{y}}^{\ast }={\gamma }_1{X}_2$
• 用$X_2$回归$X_1$：${\hat{X}}_1={\gamma }_2{X}_2$
• 分别计算残差：${\tilde{X}}_1=X_1-{\hat{X}}_1$，$\tilde{y}=y-{\hat{y}}^{\ast }$
• 用$X_1$的残差回归$y$的残差：${\hat{y}}^{\ast }={\gamma }_3{\tilde{X}}_1$

两种回归方式得到的系数是相同的，${\beta }_1={\gamma }_3$

Frisch-Waugh-Lovell Law的非线性推广：
$(Y-(Y\sim X))\sim (T-(T\sim X))⟺Y_i-E(Y_i|X_i)=\tau \ast (T_i-E(T_i|X_i))+ϵ$

Frisch-Waugh-Lovell Law解决了什么问题
Frisch-Waugh-Lovell Law可以借助机器学习模型获得$Y_i-E(Y_i|X_i)$和$T_i-E(T_i|X_i)$，无需像线性模型一样知道所有特征的正确形式，只需保证预测值${\hat{y}}^{\ast }$和$\tilde{X}$足够准确就可。

3.Partially Linear Model

• Y是实验影响的核心指标
• T是treatment，通常是0/1变量，代表样本进入实验组还是对照组，对随机AB实验𝑇⊥𝑋
• X是Confounder，可以简单理解为未被实验干预过的用户特征，通常是高维向量
  $$\begin{array}{rl}Y& ={\theta }_{0}T+g(X)+U,E(U|T,X)=0\\ T& =m(X)+V,E(V|T)=0\end{array}$$
     DML最终估计的是${\hat{\theta }}_0$，也就是实验对不同用户核心指标的不同影响。最直接的方法就是用X和T一起对Y建模，直接估计${\hat{\theta }}_0$。但这样估计出的${\hat{\theta }}_0$往往是有偏的，偏差部分来自于对样本的过拟合，部分来自于$\hat{g}(X)$估计的偏差。假定${\theta }_0$是参数的真实值，根据最小二乘公式$\theta ={(X^{T}X)}^{-1}{X}^{T}Y$，把$X=T$，$Y=Y-\hat{g}(X)$带入最小二乘公式，则估计值${\hat{\theta }}_0$为：
  θ ^ 0 = ( T T T ) − 1 T T ( Y − g ^ ( X ) ) = ∑ T i ( Y i − g ^ ( X i ) ) ∑ T i 2 = ∑ T i { [ Y i − g ( X i ) ] + [ g ( X i ) − g ^ ( X i ) ] } ∑ T i 2 = ∑ T i [ ( θ 0 T i + ϵ i ) + [ g ( X i ) − g ^ ( X i ) ] ∑ T i 2 = θ 0 + ∑ T i ϵ i + T i [ g ( X i ) − g ^ ( X i ) ] ∑ T i 2 \begin{aligned} \hat \theta_0 &= {(T^TT)}^{-1}T^T(Y-\widehat g(X)) \\ &= \frac {\sum T_i (Y_i-\widehat g(X_i))} {\sum T_i^2} \\ &= \frac {\sum T_i \{ [Y_i-g(X_i)]+ [g(X_i)-\widehat g(X_i)] \} } {\sum T_i^2} \\ &= \frac {\sum T_i [(\theta_0 T_i+\epsilon_i)+[g(X_i)-\widehat g(X_i)] } {\sum T_i^2} \\ &= \theta_0+\frac {\sum T_i \epsilon_i+ T_i [g(X_i)-\widehat g(X_i)] } {\sum T_i^2} \\ \end{aligned} θ^0​​=(TTT)−1TT(Y−g ​(X))=∑Ti2​∑Ti​(Yi​−g ​(Xi​))​=∑Ti2​∑Ti​{[Yi​−g(Xi​)]+[g(Xi​)−g ​(Xi​)]}​=∑Ti2​∑Ti​[(θ0​Ti​+ϵi​)+[g(Xi​)−g ​(Xi​)]​=θ0​+∑Ti2​∑Ti​ϵi​+Ti​[g(Xi​)−g ​(Xi​)]​​ 估计值的偏差为：
  $$\begin{array}{r}\sqrt{n}({\hat{\theta }}_0-{\theta }_0)=\sqrt{n}\ast \left(\underset{a}{\underbrace{\sum \frac{T_i{ϵ}_i}{T_i^2}}}+\underset{b}{\underbrace{\frac{T_i[g(X_i)-\hat{g}(X_i)]}{T_i^2}}}\right)\end{array}$$
     根据线性回归的假设$\sum T_i{ϵ}_i=0$，所以a项收敛为0。b项中，由于包含$T{\theta }_0$项，$g(X)\ne E[Y|X]$ ，因此，$\hat{g}(X_i)$是有偏的，得到的预测误差$g(X_i)-\hat{g}(X_i)$不会收敛为0。

4. Double Machine Learning

4.1 DML模型步骤

DML模型分为以下三个步骤：
步骤一.用任意ML模型拟合Y和T得到残差$\hat{U},\hat{V}$
$$\begin{array}{rl}\hat{U}& =Y-T{\theta }_0-\hat{g}(X_i)where\hat{g}(X_i)=E(Y|X,T)\\ \hat{V}& =T-\hat{m}(X_i)where\hat{m}(X_i)=E(T|X)\end{array}$$ 步骤二.用任意模型拟合$\hat{U},\hat{V}$的关系，得到${\hat{\theta }}_0$
${\theta }_0$的拟合可以是参数模型也可以是非参数模型，参数模型可以直接用线性模型拟合。而非参数模型因为只接受输入和输出所以需要再做如下变换，模型Target变为$\frac{\hat{U}}{\hat{V}}$, 样本权重为${\hat{V}}^2$
$$\begin{array}{rl}参数模型& ：\hat{U}={\theta }_0\ast \hat{V}+ϵargminE[(\hat{U}-{\theta }_0\ast \hat{V}{)}^2]=0\\ 非参模型& ：E[(\hat{U}-{\theta }_0\ast \hat{V}{)}^2]=E[{\hat{V}}^2\ast (\frac{\hat{U}}{\hat{V}}-{\theta }_0{)}^2]\end{array}$$ 步骤三. Cross-fitting
DML保证估计无偏很重要的一步就是Cross-fitting，用来降低overfitting带来的估计偏差。先把总样本分成两份：样本1，样本2。先用样本1估计残差，样本2估计𝜃̂ 1，再用样本2估计残差，样本1估计𝜃̂ 2，取平均得到最终的估计。当然也可以进一步使用K-Fold来增加估计的稳健性。

参数模型${\hat{\theta }}_0$的表示式：
根据两阶段最小二乘公式$\theta ={(Z^{T}X)}^{-1}{Z}^{T}Y$，把$Z=\hat{V}$，$X=T$，$Y=Y-\hat{g}(X)$带入最小二乘公式，则：
$$\begin{array}{rl}{\hat{\theta }}_0& ={({\hat{V}}^{T}T)}^{-1}{\hat{V}}^T(Y-\hat{g}(X))\\ & =\frac{\sum {\hat{V}}_i(Y_i-\hat{g}(X_i))}{\sum {\hat{V}}_i{T}_i}\\ & =\frac{\sum {\hat{V}}_i([Y_i-g(X_i)]+[g(X_i)-\hat{g}(X_i)])}{\sum {\hat{V}}_i{T}_i}\\ & =\frac{\sum {\hat{V}}_i([T_i{\theta }_0+U_i]+[g(X_i)-\hat{g}(X_i)])}{\sum {\hat{V}}_i{T}_i}\\ & ={\theta }_0+\frac{\sum {\hat{V}}_i{U}_i+\sum {\hat{V}}_i[g(X_i)-\hat{g}(X_i)])}{\sum {\hat{V}}_i{T}_i}\\ & ={\theta }_0+\frac{\sum {\hat{V}}_i{U}_i+\sum [V_i-(\hat{m}(X_i)-m(X_i))][g(X_i)-\hat{g}(X_i)])}{\sum {\hat{V}}_i{T}_i}\\ & ={\theta }_0+\underset{a}{\underbrace{\frac{\sum {\hat{V}}_i{U}_i}{\sum {\hat{V}}_i{T}_i}}}+\underset{b}{\underbrace{\frac{\sum [m(X_i)-\hat{m}(X_i)][g(X_i)-\hat{g}(X_i)]}{\sum {\hat{V}}_i{T}_i}}}+\underset{c}{\underbrace{\frac{\sum V_i[g(X_i)-\hat{g}(X_i)]}{\sum {\hat{V}}_i{T}_i}}}\end{array}$$ a项等于0；b项由于$\sum [m(X_i)-\hat{m}(X_i)]$=0，且$\sum [m(X_i)-\hat{m}(X_i)]$和$\sum [g(X_i)-\hat{g}(X_i)]$独立，所以b项等于0；c项中，$[g(X_i)-\hat{g}(X_i)]$是方程$Y=T{\hat{\theta }}_0+\hat{g}(X)+\hat{U}$的预测误差，如果$\hat{g}(x)$过拟合，预测误差中会包含噪声$U$，如果噪声$U和V$相关联，那么c项不会等于0（$V$和$[g(X_i)-\hat{g}(X_i)]$相关联，不独立）。为了解决这个问题，需要用到步骤三的样本分割方法。

如果有任何未观测的confounder，$V$和$[g(X_i)-\hat{g}(X_i)]$相关联，但我们假设观测到了所有confounders

另一种推导方式( g(X) 换成p(X) )：
根据Frisch-Waugh-Lovell Law，DML也可以表示成如下形式：
$$\begin{array}{rl}Y& =p(X)+W,p(X)=E(Y|X),\hat{W}=Y-\hat{p}(X)\\ T& =m(X)+V,m(X)=E(T|X),\hat{V}=T-\hat{m}(X)\\ \hat{W}& =\hat{V}{\theta }_0+\eta \end{array}$$ ${\hat{\theta }}_0$的表达式为：
$$\begin{array}{rl}{\hat{\theta }}_0& =\frac{\sum {\hat{V}}_i{\hat{W}}_i}{\sum {{\hat{V}}_i}^2}\\ & =\frac{\sum [V_i+(m(X_i)-\hat{m}(X_i))][W_i+(p(X_i)-\hat{p}(X_i))]}{\sum {{\hat{V}}_i}^2}\\ & =\frac{\sum V_i{W}_i+[m(X_i)-\hat{m}(X_i)][p(X_i)-\hat{p}(X_i)]+W_i[m(X_i)-\hat{m}(X_i)]+V_i[p(X_i)-\hat{p}(X_i)]}{\sum {{\hat{V}}_i}^2}\\ & =\underset{a}{\underbrace{\frac{\sum V_i{W}_i}{\sum {{\hat{V}}_i}^2}}}+\underset{b}{\underbrace{\frac{\sum [m(X_i)-\hat{m}(X_i)][p(X_i)-\hat{p}(X_i)]}{\sum {{\hat{V}}_i}^2}}}+\underset{c}{\underbrace{\frac{\sum W_i[m(X_i)-\hat{m}(X_i)]}{\sum {{\hat{V}}_i}^2}}}+\underset{{\theta }_0}{\underbrace{\frac{\sum V_i[p(X_i)-\hat{p}(X_i)]}{\sum {{\hat{V}}_i}^2}}}\end{array}$$
为什么$$\frac{\sum V_i[p(X_i)-\hat{p}(X_i)]}{\sum {{\hat{V}}_i}^2}={\theta }_0???$$

4.2 GMM角度理解DML

PLM（Partially Linear Model）的总体矩条件：
$$\begin{gathered} {\phi }_{plm}(W;{\theta }_0,\eta )=[Y-T{\theta }_0-g(X)]T \\ where\eta =g(X),{\eta }_0=g_0(X)\text{\hspace{1em}(1)} \end{gathered}$$ DML的总体矩条件：
$$\begin{array}{rl}{\phi }_{dml}(W;{\theta }_0,\eta )& =[Y-T{\theta }_0-g(X)][T-m(X)]\\ or：{\phi }_{dml}(W;{\theta }_0,\eta )& =[Y-p(X)-(T-m(X)){\theta }_0][T-m(X)]\\ where\eta & =(p(X),m(X)),{\eta }_0=(p_0(X),m_0(X))=(E(Y|X),E(T|X))\end{array}\text{\hspace{1em}(2)}$$
根据Neyman orthogonality：PML矩条件的求导结果不等于0，
而DML矩条件的加托求导(Gateaux derivative)结果是0：
∂ η E { [ Y − p ( X ) − ( T − m ( X ) ) θ 0 ] [ T − m ( X ) ] } ∣ η = η 0 = 0 \partial_\eta E\{ [ Y - p(X)- (T-m(X))\theta_0 ] [T-m(X)] \} |_{\eta=\eta_0} = 0 ∂η​E{[Y−p(X)−(T−m(X))θ0​][T−m(X)]}∣η=η0​​=0
这表明用于估计parameters of interest的矩条件, 对nuisance parameters 的导数等于0 , 换句话说，即矩条件在true value of nuisance parameters 处对 nuisance parameters 的导数与从true value of nuisance parameters指向nuisance parameters的向量正交(which indicates the terminology of orthogonality)。由此得到的estimator 也叫 doubly robust estimator，指对于parameters of insterest 和 nuisance parameters 都是稳健的。看不懂。。。。

4.3 DML正交化实例说明

步骤一.原始数据
  探索冰淇淋销量和价格之间的关系，可以看到weekday是confounder，周末时价格和销量均更高。如果要得到准确的causal effect，需要矫正confounders带来的偏差。此数据集中包含的confounders有温度、周几、成本

步骤二.regress T on X
  用X构建T的回归模型，得到残差$\hat{V}=T-\hat{m}(X)$，残差$\hat{V}$可以看作将X对T的作用从T中去除后剩下的量，此时残差$\hat{V}$正交于$X$。

步骤三.regress Y on X
  同样的，用X构建Y的回归模型，得到残差$\hat{W}=Y-\hat{p}(X)$，而模型$\hat{p}(X)$的作用在于去除Y的方差，即将X引起的Y的方差从Y中去除。和步骤二图相比，步骤三图Y的方差变小。

步骤四. 对比DML结果和PML结果
  从上图中可以看到price_res和sales_res呈明显负相关关系。分别构建price_res和sales_res线性模型以及price和sales线性模型，得到前者得到显著负相关关系，后者由于混淆变量的存在得到正相关关系。

4.4 DML模拟

假设
$$\begin{array}{rl}g(X)& =si{n}^2(X)\\ m(X)& =\frac{1}{2\pi }\frac{sinh(\gamma )}{cosh(\gamma )-cosh(x-v)}\end{array}$$ 令$v=0,\gamma =1$，$X$和$m(X)$的关系为

零.生成数据  ${\theta }_0$真实值=0.5

import numpy as np
from sklearn.datasets import make_spd_matrix
import math
import statsmodels.api as sm 
from sklearn.ensemble import RandomForestRegressor 
randomseednumber = 11022018
np.random.seed(randomseednumber)
import warnings 
warnings.filterwarnings('ignore')
import seaborn as sns
import matplotlib.pyplot as plt

N = 500 # No. obs 
k=10 # = No. variables in x_i 
theta=0.5 # Structural parameter,the true value of the causal parameter will be θ=0.5, and bk=1/k. 
b= [1/k for k in range(1,11)] # x weights 
# sigma = make_spd_matrix(k,randomseednumber) #shape(10, 10)
sigma = make_spd_matrix(k) #shape(10, 10)
MC_no = 500 #NUmber of simulations 

def g(x):
    return np.power(np.sin(x),2)

def m(x,nu=0.,gamma=1.):
    return 0.5/math.pi*(np.sinh(gamma))/(np.cosh(gamma)-np.cos(x-nu))

# Array of estimated thetas to store results 
theta_est = np.zeros(shape=[MC_no,3]) #shape(500, 3)

一. 真实值： ${\theta }_0$真实值=0.5 sm.OLS(Y-G,D)得到${\theta }_0$的真实值

theta_est_true = np.zeros(shape=[MC_no,1]) #shape(500, 1)

for i in range(MC_no):
    X = np.random.multivariate_normal(np.ones(k),sigma,size=[N,]) #shape(500, 10)
    G = g(np.dot(X,b)) #shape(500,)
    M = m(np.dot(X,b)) #shape(500,)
    D = M+np.random.standard_normal(size=[500,])
    Y = np.dot(theta,D)+G+np.random.standard_normal(size=[500,])

    # true value of the causal parameter is θ=0.5
    model = sm.OLS(Y-G,D)
    results = model.fit()
    theta_est_true[i][0] = results.params[0] # x1的coef,即theta0
#     theta_est_true[i][0] = np.sum(np.dot(D,Y-G))/np.sum(np.dot(D.T,D)) #两者等价
    
sns.distplot(theta_est_true[:,0],kde=True)
plt.axvline(0.5,0,5)

二. PLM方法估计值：

• 设定${\theta }_0$的初始值 ${\widehat{{\theta }_0}}^{(0)}=0.4$
• 建立$\mathbf{X}$ 与$Y-D{\widehat{{\theta }_0}}^{(0)}$ 之间的机器学习模型，以估计${\widehat{g_0}}^{(1)}$
• 建立D 与$Y-{\widehat{g_0}}^{(1)}(\mathbf{X})$之间的线性回归，以估计${\widehat{{\theta }_0}}^{(1)}$
• 重复上述方法，直到 $|{\widehat{{\theta }_0}}^{(i)}-{\widehat{{\theta }_0}}^{(i-1)}|<ϵ$ 收斂，得到最终的估计$\widehat{g_0}$与$\widehat{{\theta }_0}$

theta_est_plm = np.zeros(shape=[MC_no,1]) #shape(500, 3)

for i in range(MC_no):
    theta_0,theta_delta = 0.4,0.1
    X = np.random.multivariate_normal(np.ones(k),sigma,size=[N,]) #shape(500, 10)
    G = g(np.dot(X,b)) #shape(500,)
    M = m(np.dot(X,b)) #shape(500,)
    D = M+np.random.standard_normal(size=[500,])
    Y = np.dot(theta_0,D)+G+np.random.standard_normal(size=[500,])
    while theta_delta>0.01:
        theta_pre = theta_0

        # partial linear model
        naiveDMLg = RandomForestRegressor(max_depth=2)
        naiveDMLg.fit(X,Y-np.dot(D,theta_0)) #Compute ghat  
        Ghat = naiveDMLg.predict(X) #shape(500,)

        theta_0 = np.sum(np.dot(D,Y-Ghat))/np.sum(np.dot(D,D))
        theta_delta = abs(theta_pre-theta_0)

    theta_est_plm[i][0] = theta_0

sns.distplot(theta_est_plm[:,0],kde=True)
plt.axvline(0.5,0,5)

三. Naive DML方法估计值：np.mean(np.dot(Vhat,Y-Ghat))/np.mean(np.dot(Vhat,D))得到${\theta }_0$的预估值，有偏。

for i in range(MC_no):
#     Generate data: no. obs x no. variables in x_i    
#     xi has length K=10 and will be generated from a multivariate normal distribution,
    X = np.random.multivariate_normal(np.ones(k),sigma,size=[N,]) #shape(500, 10)
    G = g(np.dot(X,b)) #shape(500,)
    M = m(np.dot(X,b)) #shape(500,)
    D = M+np.random.standard_normal(size=[500,])
    Y = np.dot(theta,D)+G+np.random.standard_normal(size=[500,])

    # Naive double machine Learning
    naiveDMLg = RandomForestRegressor(max_depth=2)
    naiveDMLg.fit(X,Y) #Compute ghat  
    Ghat = naiveDMLg.predict(X) #shape(500,)
    naiveDMLm =RandomForestRegressor(max_depth=2)
    naiveDMLm.fit(X,D)
    Mhat= naiveDMLm.predict(X) #Compute Mhat,Vhat as residual
    Vhat = D-Mhat
    theta_est[i][1] = np.mean(np.dot(Vhat,Y-Ghat))/np.mean(np.dot(Vhat,D))
    
sns.distplot(theta_est[:,1],kde=True)
plt.axvline(0.5,0,5)    

四. cross-fitting DML方法估计值：sample spliting 无偏。

for i in range(MC_no):
#     Generate data: no. obs x no. variables in x_i    
#     xi has length K=10 and will be generated from a multivariate normal distribution,
    X = np.random.multivariate_normal(np.ones(k),sigma,size=[N,]) #shape(500, 10)
    G = g(np.dot(X,b)) #shape(500,)
    M = m(np.dot(X,b)) #shape(500,)
    D = M+np.random.standard_normal(size=[500,])
    Y = np.dot(theta,D)+G+np.random.standard_normal(size=[500,])

    # Cross-fitting DML
    I = np.random.choice(N,np.int(N/2),replace=False) # Split the sample 
    I_C = [x for x in np.arange(N) if x not in I]
    # Ghat for both 
    Ghat_1 = RandomForestRegressor(max_depth=2).fit(X[I],Y[I]).predict(X[I_C])
    Ghat_2 = RandomForestRegressor(max_depth=2).fit(X[I_C],Y[I_C]).predict(X[I])
    # Mhat and vhat for both     
    Mhat_1 = RandomForestRegressor(max_depth=2).fit(X[I],D[I]).predict(X[I_C])
    Mhat_2 = RandomForestRegressor(max_depth=2).fit(X[I_C],D[I_C]).predict(X[I])
    Vhat_1,Vhat_2 = D[I_C]-Mhat_1,D[I] - Mhat_2
    theta_1 = np.mean(np.dot(Vhat_1,(Y[I_C]-Ghat_1)))/np.mean(np.dot(Vhat_1,D[I_C]))
    theta_2 = np.mean(np.dot(Vhat_2,(Y[I]-Ghat_2)))/np.mean(np.dot(Vhat_2,D[I]))
    theta_est[i][2] = 0.5*(theta_1+theta_2)

sns.distplot(theta_est[:,2],kde=True)
plt.axvline(0.5,0,5) 

五. 结果对比  下图是四种方法的估计值分布。

4.5 非参数DML 累了 写不动

参考资料

DML原始论文
DML详细ppt
DML详细讲解论文
DML训练过程
DML 例子
AB实验人群定向HTE模型4 - Double Machine Learning
DML中文推导
因果推断笔记——DML ：Double Machine Learning案例学习（十六）