因果推断笔记——DML ：Double Machine Learning案例学习（十六）
文章目录
1 背景知识
参考：
Hetergeneous Treatment Effect旨在量化实验对不同人群的差异影响，进而通过人群定向/数值策略的方式进行差异化实验，或者对实验进行调整。Double Machine Learning把Treatment作为特征，通过估计特征对目标的影响来计算实验的差异效果。
Machine Learning擅长给出精准的预测，而经济学更注重特征对目标影响的无偏估计。DML把经济学的方法和机器学习相结合，在经济学框架下用任意的ML模型给出特征对目标影响的无偏估计，有知乎大佬的一个总结：
HTE问题可以用以下的notation进行简单的抽象
- Y是实验影响的核心指标
- T是treatment，通常是0/1变量，代表样本进入实验组还是对照组，对随机AB实验T⊥X
- X是Confounder，可以简单理解为未被实验干预过的用户特征，通常是高维向量
- DML最终估计的是θ(x)，也就是实验对不同用户核心指标的不同影响
最直接的方法就是用X和T一起对Y建模，直接估计θ(x)。 但这样估计出的θ(x)往往是有偏的，偏差部分来自于对样本的过拟合，部分来自于g(X)估计的偏差，假定θ0是参数的真实值，则偏差如下：
2 DML Dynamic Double Machine Learning
【因果推断/uplift建模】Double Machine Learning(DML)
2.1 DML训练过程
参考： AB实验人群定向HTE模型4 - Double Machine Learning 训练的过程： 步骤一. 用任意ML模型拟合Y和T得到残差Y ,T
比较常见的是使用lasso / RF进行E(Y|x)的估计
步骤二. 对Y,T用任意ML模型拟合θ
θ(X)的拟合可以是参数模型也可以是非参数模型，参数模型可以直接拟合。而非参数模型因为只接受输入和输出所以需要再做如下变换，模型Target变为Y/T, 样本权重为T^2
步骤三. Cross-fitting
这个Cross-fitting步骤非常重要，DML保证估计无偏很重要的一步就是Cross-fitting，用来降低overfitting带来的估计偏差。先把总样本分成两份：样本1，样本2。 先用样本1估计残差，样本2估计θ1，再用样本2估计残差，样本1估计θ2，取平均得到最终的估计。当然也可以进一步使用K-Fold来增加估计的稳健性。 Jonas在他的博客里比较了不使用DML，使用DML但是不用Cross-fitting，以及使用Cross-fitting的估计效果如下：
2.2 HTL无偏估计
直观角度上，这里对回归比较熟悉的朋友可以知道，线性回归是拟合Y在特征空间X的最佳投影（既误差最小化），所以残差是垂直样本空间X的，既最大限度消除了（独立）X的相关性！如下图所示。
关于残差正交化可得到无偏差因果效应的数学原理：https://zhuanlan.zhihu.com/p/41993542
2.3 使用DML估计ATE
具体到因果推断的例子上，我们只关心Treatment T 对 outcome Y的影响，因此我们可以首先使用X回归T，得到一个T的残差（实际T - 预测T），然后使用X回归Y，得到一个Y的残差（实际Y - 预测Y），最后使用T的残差回归Y的残差，估计的参数即我们想要的ATE。
于是乎， 具体的DML方法也就出来了，其核心思想即分别用机器学习算法基于X预测T和Y，然后使用T的残差回归Y的残差：
Y_i - \hat{M}_y(X_i) = \tau \cdot (T_i - \hat{M}_t(X_i)) + \epsilon 其中， \hat{M}_y(X_i) 建模 E[Y|X] ， \hat{M}_t(X_i) 建模 E[T|X] 。那么问题来了，为什么说DML能去偏呢？其关键在于 \hat{M}_t(X_i) ，它对Treatment实施了去偏。T的残差可以看作将X对T的作用从T中去除后剩下的量，此时T的残差独立于X。而 \hat{M}_y(X_i) 的作用在于去除Y的方差，即将X引起的Y的方差从Y中去除。最后，我们再对残差建模lr，即得到ATE。
2.4 使用DML估计CATE
同样地，我们首先基于X使用ML获得T的残差和Y的残差，之后使用lr拟合残差，不同的是，这次我们把X和T的交互项加进来，即
然后我们就可以计算CATE的值了： \hat{\mu}(\partial Sales_i, X_i) = M(Price=1, X_i) - M(Price=0, X_i) 其中，M即最后的lr模型。
2.5 直接预测反事实的Y
在处理非线性CATE时，另一种方案是我们将不再尝试估计CATE的线性近似。相反，我们将做出反事实的预测。 (这种方案在实际中使用也很多，但并没有严格的理论证明！)
其流程分为两个步骤：
第一步，依然是估计T和Y的残差： \tilde{Y}_i = \tau(X_i) \tilde{T}_i + e_i
第二步，基于X和T的残差使用S-learner预测Y的残差： \tilde{Y}_i = \tau(X_i, \tilde{T}_i) + e_i 最后在 \hat{M}_y(X_i) 预测的 \hat{Y} 上加上 \tilde{Y}_i ，即得到最后的y值。
3 相关案例
3.1 一个DML实现的价格弹性的案例
3.1.1 价格需求弹性
经济学课程里谈到价格需求弹性，描述需求数量随商品价格的变动而变化的弹性。价格一般不直接影响需求，而是被用户决策相关的中间变量所中介作用。假设 Q 为某个商品的需求的数量，P 为该商品的价格，则计算需求的价格弹性为，
通过上式可以简单知道，价格改变 1 元比价格改变 100 元，会导致更大的需求改变。比如以 5 元的价格每日可以卖 100 单位产品，如果价格需求弹性为 -3 ，那供应商将价格提升 5%（dp /P，从 5 元-> 5.25 元），需求将下降 15%（dQ/Q ，从 100->85）。那么收入将减少 100 5-5.25 85=53.75。 如果单价降低 5%，那么收入同理将提升 46.25。如果供应商知道了产品的价格弹性，那无须反复测试，即可清楚为提升收入到底应该是提价还是降价。
3.1.2 DML（Double Machine Learning） 求解过程
最好的方式，当然是直接进行 A/B 实验测试不同价格对用户的需求反应，但是价格这类的外生因素在同一产品同一阶段上，对不同用户展示不同的价格会直接损坏用户体验。因此从观察历史数据进行因果推断，但混杂因素（季节性、产品质量等）如何控制是因果推断的挑战。
这里采用 DML（Double Machine Learning） 方法进行因果推断，该方法主要解决两个问题：
- 第一，通过正则化挑拣重要控制变量；
- 第二，对比传统的线性回归模型，用非参数推断可以解决非线性问题。
DML 先应用机器学习算法去分别通过特征变量 X, W 拟合结果变量 Y 和处理变量 T，然后通过线性模型，使用处理变量的残差拟合出结果变量的残差。 目标是估计 ，这里的 Y 函数构成为 T 的因果作用和 X、W 的协同作用之和。
将数据分为两部分，
- 一部分样本选用随机森林等模型，用混杂变量预测处理变量（价格 P），得到 E[P|X]；
- 另外的样本同样可选择随机森林模型，用混杂变量预测结果变量（需求量 Q），得到 E[Q|X]。 计算残差，得到不受混杂变量影响的价格 P 和 需求量 Q，即为：
关于残差正交化可得到无偏差因果效应的数学原理：https://zhuanlan.zhihu.com/p/41993542
因此直接将 P,Q 进行 log-log 回归就能得到弹性系数 ：
关于为何用 log-log 回归可求解弹性系数：
3.2 因果分析工具在快手的应用：DML篇
文章来源地址： 因果分析工具在快手的应用
快手这篇我觉得比较有意思的是他们给出了适用例子：遗漏变量是否可量化
- 如果可量化，适合DML
具体PPT见：
4 econml的notebook案例
4.1 其他理论介绍
参考： 期刊：因果推理
4.1.1 surrogate indices 代替指标
试图利用一个近期的或容易得到的中间指标来替代远期的或难以得到的终点指标。
抑制心律失常能降低发生心脏骤停的可能性，曾将[抑制心律失常]作为评价[治疗猝死药物]的替代指标。
因为混杂因子而造成了很多的悖论：
- Yule-Simpson 悖论:吸烟分别对男性和对女性都有害，但是吸烟对人类有益
- 替代指标悖论:临床试验将心律失常作为猝死的替代指标。但是，一些能有效纠正心律失常的药物，后来发现不但不能减少猝死，反而导致数万人过早死亡
4.1.2 动态因果效应的估计
https://www.jianshu.com/p/871fbf457a6d
在这个框架内，单个对象在不同时间点上分别扮演了处理组和对照组两种角色。
当X_t严格外生时，有两个可选的动态因果效应估计量：
- 第一个方法涉及自回归分布滞后模型（ADL）
- 第二种方法利用广义最小二乘法（GLS）
4.2 Double Machine Learning Notebook
参考官方： Double Machine Learning Notebook
该案例有非常多小的案例：
- Example Usage with Single Continuous Treatment Synthetic Data and Model Selection
- Example Usage with Single Binary Treatment Synthetic Data and Confidence Intervals
- Example Usage with Multiple Continuous Treatment Synthetic Data
- Example Usage with Single Continuous Treatment Observational Data
- Example Usage with Multiple Continuous Treatment, Multiple Outcome Observational Data
其中包括了：连续数值的T，二分类01的T，多个T，多个Y
# 连续数值的T
array([ 0.41083324, 0.57893171, 1.08130798, ..., -0.43799495,
1.61941775, 1.64209826])
# 二分类01的T
array([0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0])
# Multi-T
array([[1.35325451, 1.15373159, 0.46373402],
[1.35325451, 1.15373159, 0.98954119],
[1.35325451, 0.87129337, 0.73716407],
[1.13462273, 1.03673688, 0.46373402],
[1.0260416 , 0.95165788, 0.39877612],
[1.05779029, 0.78390154, 0.55961579]])
# Multi-Y
array([[ 9.01869549, 8.40737833, 9.26482856],
[ 8.72323127, 8.44934252, 8.98719682],
[ 8.25322765, 9.91145572, 8.83171192],
[ 9.93653541, 9.51546936, 9.50599061],
[10.85591733, 9.31793838, 10.9273763 ],
[10.23652533, 11.20553036, 8.85936345]])方法包括：
- LinearDML：默认
- SparseLinearDML：Polynomial Features for Heterogeneity
- DML：Polynomial Features with regularization
- CausalForestDML： Non-Parametric Heterogeneity with Causal Forest
包括数据生成：
import econml
## Ignore warnings
import warnings
warnings.filterwarnings("ignore")
# Main imports
from econml.dml import DML, LinearDML, SparseLinearDML, CausalForestDML
# Helper imports
import numpy as np
from itertools import product
from sklearn.linear_model import (Lasso, LassoCV, LogisticRegression,
LogisticRegressionCV,LinearRegression,
MultiTaskElasticNet,MultiTaskElasticNetCV)
from sklearn.ensemble import RandomForestRegressor,RandomForestClassifier
from sklearn.preprocessing import PolynomialFeatures
import matplotlib.pyplot as plt
import matplotlib
from sklearn.model_selection import train_test_split
%matplotlib inline
# Treatment effect function
def exp_te(x):
return np.exp(2*x[0])
# DGP constants
np.random.seed(123)
n = 2000
n_w = 30
support_size = 5
n_x = 1
# Outcome support
support_Y = np.random.choice(np.arange(n_w), size=support_size, replace=False)
coefs_Y = np.random.uniform(0, 1, size=support_size)
epsilon_sample = lambda n: np.random.uniform(-1, 1, size=n)
# Treatment support
support_T = support_Y
coefs_T = np.random.uniform(0, 1, size=support_size)
eta_sample = lambda n: np.random.uniform(-1, 1, size=n)
# Generate controls, covariates, treatments and outcomes
W = np.random.normal(0, 1, size=(n, n_w))
X = np.random.uniform(0, 1, size=(n, n_x))
# Heterogeneous treatment effects
TE = np.array([exp_te(x_i) for x_i in X])
T = np.dot(W[:, support_T], coefs_T) + eta_sample(n)
Y = TE * T + np.dot(W[:, support_Y], coefs_Y) + epsilon_sample(n)
Y_train, Y_val, T_train, T_val, X_train, X_val, W_train, W_val = train_test_split(Y, T, X, W, test_size=.2)
# Generate test data
X_test = np.array(list(product(np.arange(0, 1, 0.01), repeat=n_x)))模型：
# 模型
est = LinearDML(model_y=RandomForestRegressor(),
model_t=RandomForestRegressor(),
random_state=123)
est.fit(Y_train, T_train, X=X_train, W=W_train)
te_pred = est.effect(X_test)
est1 = SparseLinearDML(model_y=RandomForestRegressor(),
model_t=RandomForestRegressor(),
featurizer=PolynomialFeatures(degree=3),
random_state=123)
est1.fit(Y_train, T_train, X=X_train, W=W_train)
te_pred1 = est1.effect(X_test)
est2 = DML(model_y=RandomForestRegressor(),
model_t=RandomForestRegressor(),
model_final=Lasso(alpha=0.1, fit_intercept=False),
featurizer=PolynomialFeatures(degree=10),
random_state=123)
est2.fit(Y_train, T_train, X=X_train, W=W_train)
te_pred2 = est2.effect(X_test)
est3 = CausalForestDML(model_y=RandomForestRegressor(),
model_t=RandomForestRegressor(),
criterion='mse', n_estimators=1000,
min_impurity_decrease=0.001,
random_state=123)
est3.tune(Y_train, T_train, X=X_train, W=W_train)
est3.fit(Y_train, T_train, X=X_train, W=W_train)
te_pred3 = est3.effect(X_test)
plt.figure(figsize=(10,6))
plt.plot(X_test, te_pred, label='DML default')
plt.plot(X_test, te_pred1, label='DML polynomial degree=3')
plt.plot(X_test, te_pred2, label='DML polynomial degree=10 with Lasso')
plt.plot(X_test, te_pred3, label='ForestDML')
expected_te = np.array([exp_te(x_i) for x_i in X_test])
plt.plot(X_test, expected_te, 'b--', label='True effect')
plt.ylabel('Treatment Effect')
plt.xlabel('x')
plt.legend()
plt.show()
后面包括模型选择：
score={}
score["DML default"] = est.score(Y_val, T_val, X_val, W_val)
score["DML polynomial degree=3"] = est1.score(Y_val, T_val, X_val, W_val)
score["DML polynomial degree=10 with Lasso"] = est2.score(Y_val, T_val, X_val, W_val)
score["ForestDML"] = est3.score(Y_val, T_val, X_val, W_val)
mse_te={}
mse_te["DML default"] = ((expected_te - te_pred)**2).mean()
mse_te["DML polynomial degree=3"] = ((expected_te - te_pred1)**2).mean()
mse_te["DML polynomial degree=10 with Lasso"] = ((expected_te - te_pred2)**2).mean()
mse_te["ForestDML"] = ((expected_te - te_pred3)**2).mean()
print("best model selected by MSE of TE: ", min(mse_te, key=lambda x: mse_te.get(x)))score指的是MSE，越小越好
4.3 Dynamic Double Machine Learning Examples
参考： Dynamic Double Machine Learning Examples.ipynb
Dynamic DoubleML is an extension of the Double ML approach for treatments assigned sequentially over time periods. This estimator will account for treatments that can have causal effects on future outcomes.
相当于DML + 时序效应，所以一般都是面板数据才会需要 生成数据
import econml
# Main imports
from econml.dynamic.dml import DynamicDML
from econml.tests.dgp import DynamicPanelDGP, add_vlines
# Helper imports
import numpy as np
from sklearn.linear_model import Lasso, LassoCV, LogisticRegression, LogisticRegressionCV, MultiTaskLassoCV
import matplotlib.pyplot as plt
%matplotlib inline
# Define DGP parameters
np.random.seed(123)
n_panels = 5000 # number of panels
n_periods = 3 # number of time periods in each panel
n_treatments = 2 # number of treatments in each period
n_x = 100 # number of features + controls
s_x = 10 # number of controls (endogeneous variables)
s_t = 10 # treatment support size
# Generate data
dgp = DynamicPanelDGP(n_periods, n_treatments, n_x).create_instance(
s_x, random_seed=12345)
Y, T, X, W, groups = dgp.observational_data(n_panels, s_t=s_t, random_seed=12345)
true_effect = dgp.true_effect可以看到面板数据，有三年数据，5000个样本，两个处理T
模型训练
# Train Estimator
est = DynamicDML(
model_y=LassoCV(cv=3, max_iter=1000),
model_t=MultiTaskLassoCV(cv=3, max_iter=1000),
cv=3)
est.fit(Y, T, X=None, W=W, groups=groups)
# Effect of target policy over baseline policy
# Must specify a treatment for each period
baseline_policy = np.zeros((1, n_periods * n_treatments))
target_policy = np.ones((1, n_periods * n_treatments))
eff = est.effect(T0=baseline_policy, T1=target_policy)
print(f"Effect of target policy over baseline policy: {eff[0]:0.2f}")
# Period treatment effects + interpretation
for i, theta in enumerate(est.intercept_.reshape(-1, n_treatments)):
print(f"Marginal effect of a treatments in period {i+1} on period {n_periods} outcome: {theta}")
# Period treatment effects with confidence intervals
est.summary()
conf_ints = est.intercept__interval(alpha=0.05)可以知道CATE 的情况：
画图：
# Some plotting boilerplate code
plt.figure(figsize=(15, 5))
plt.errorbar(np.arange(n_periods*n_treatments)-.04, est.intercept_, yerr=(conf_ints[1] - est.intercept_,
est.intercept_ - conf_ints[0]), fmt='o', label='DynamicDML')
plt.errorbar(np.arange(n_periods*n_treatments), true_effect.flatten(), fmt='o', alpha=.6, label='Ground truth')
for t in np.arange(1, n_periods):
plt.axvline(x=t * n_treatments - .5, linestyle='--', alpha=.4)
plt.xticks([t * n_treatments - .5 + n_treatments/2 for t in range(n_periods)],
["$\\theta_{}$".format(t) for t in range(n_periods)])
plt.gca().set_xlim([-.5, n_periods*n_treatments - .5])
plt.ylabel("Effect")
plt.legend()
plt.show()
4.4 Long-Term Return-on-Investment at Microsoft via Short-Term Proxies
Case Study - Long-Term Return-on-Investment at Microsoft via Short-Term Proxies.ipynb
这里有一个关于DynamicDML的非常难懂的案例，笔者自己看的时候，觉得非常费解。这里也涉及到了【4.1】里面的代替指标。 在这简单根据笔者的理解，进行一下总结，生成数据：
# imports
from econml.data.dynamic_panel_dgp import SemiSynthetic
from sklearn.linear_model import LassoCV, MultiTaskLassoCV
import numpy as np
import matplotlib.pyplot as plt
# generate historical dataset (training purpose)
np.random.seed(43)
dgp = SemiSynthetic()
dgp.create_instance()
n_periods = 4 # Y包含了几个时期，默认为4个时间段
n_units = 5000 # 样本量
n_treatments = dgp.n_treatments # 3 - 多干预
random_seed = 43
thetas = np.random.uniform(0, 2, size=(dgp.n_proxies, n_treatments)) # 误差项
panelX, panelT, panelY, panelGroups, true_effect = dgp.gen_data(
n_units, n_periods, thetas, random_seed
)这里有四年（时期）的数据，每个时期5k样本量
4.4.1 通过DynamicDML获得无处理的控制Y
DynamicDML 结果为uplift增量，Y真实值 - DynamicDML增量 = Y常规趋势
这里有点像时间序列分解中，拆分为：长期趋势，季节变动，循环波动，不规则波动
EconML’s DynamicDML estimator is an extension of Double Machine Learning approach to dynamically estimate the period effect of treatments assigned sequentially over time period. In this scenario, it could help us to adjust the cumulative revenue by subtracting the period effect of all of the investments after the target investment.
# on historical data construct adjusted outcomes
from econml.dynamic.dml import DynamicDML
panelYadj = panelY.copy()
est = DynamicDML(
model_y=LassoCV(max_iter=2000), model_t=MultiTaskLassoCV(max_iter=2000), cv=2
for t in range(1, n_periods): # for each target period 1...m
# learn period effect for each period treatment on target period t
est.fit(
long(panelY[:, 1 : t + 1]),
long(panelT[:, 1 : t + 1, :]), # reshape data to long format
X=None,
W=long(panelX[:, 1 : t + 1, :]),
groups=long(panelGroups[:, 1 : t + 1]),
# print(f'y的shape {long(panelY[:, 1 : t + 1]).shape} , t 为 {t}， x为{long(panelX[:, 1 : t + 1, :]).shape}')
# remove effect of observed treatments
T1 = wide(panelT[:, 1 : t + 1, :])
panelYadj[:, t] = panelY[:, t] - est.effect(
T0=np.zeros_like(T1), T1=T1
) # reshape data to wide format其中，修复每一期的过程中，用数据训练的DynamicDML 也不一致； 比如第一期periods，只有5000样本，就先拿5000样本； 第二期periods是把第一期+第二期的5000样本，一共10000样本。
整个案例非常费解，是因为这些数据没有现实意义，笔者也只能按照理解， 假设这个数据是，预测银行加息、降息、新增衍生品（T1/T2/T3，都是连续值）以及人口、经济（X/W）等因素，在4年内，对5000个国家GDP（Y）的影响 第一个剔除的操作，就是把 银行加息、降息、新增衍生品的影响去掉，剩下了单独人口/经济等因素，然后来求解 调整后的GDP
4.4.2 短期时序作为长期时序的代替指标
# 单独训练一个Y ~ X的回归模型进行预测
# train surrogate index on historical dataset
XS = np.hstack(
[panelX[:, 1], panelYadj[:, :1]]
) # concatenate controls and surrogates from historical dataset
# 有四期，第2个时期的X / Y
# (5000, 71) + (5000, 1) => (5000, 72)
TotalYadj = np.sum(panelYadj, axis=1) # total revenue from historical dataset
adjusted_proxy_model = LassoCV().fit(
XS, TotalYadj
) # train proxy model from historical dataset
# predict new long term revenue
XSnew = np.hstack(
[panelXnew[:, 1], panelYnew[:, :1]]