双重机器学习

我们之前就提到，我们如何控制混杂因素，在线性回归里，我们将其纳入模型里进行控制。我们仍然用的例子是雪糕销量与定价，在这个例子里，我们还考虑到了假期，温度，以及成本的混杂因素

注意我们只对参数 $\tau$ 感兴趣，因为那是我们的处理效应。其他参数称为干扰（nuisance）参数，因为我们并不关心它们。不过，事实证明，即使我们不关心它们，也必须正确估计它们，否则处理效应就会偏离真实值

又想想，在温度上，随着气温升高，其雪糕的销量会随之升高，但如果温度太高了，那么由于天气太热，人们出门买雪糕的机会就减少了，所以销量又下降了，因此温度这个变量可能不是线性的，而是一个开口向下的二次函数。

所以，对于干扰参数的建模需要很谨慎的考虑，当我们的协变量的数量随之升高时候，这个时候就会特别复杂了！

Frisch-Waugh-Lovell 定理

这一章，我们在面板数据那里介绍过。这个定理是理解正交机器学习的关键：

我们回忆一下步骤：

假设现在有一个普通的线性回归：

1. 用 $X_2$ 去预测结果变量 $Y$：

这一步会得到一个预测值 $\hat{Y}=X_2\hat{\gamma}_y$。它表示的是：只根据 $X_2$，我们能解释 $Y$ 的哪一部分。

2. 用 $X_2$ 去预测我们关心的变量 $X_1$：

这一步会得到一个预测值 $\hat{X}_1=X_2\hat{\gamma}_x$。它表示的是：只根据 $X_2$，我们能解释 $X_1$ 的哪一部分。

3. 分别计算 $Y$ 和 $X_1$ 中不能被 $X_2$ 解释的残差：

其中 $r_y$ 表示排除了 $X_2$ 之后，$Y$ 剩下的变化；$r_x$ 表示排除了 $X_2$ 之后，$X_1$ 剩下的变化。

4. 用 $r_x$ 去解释 $r_y$：

在这里，我们需要有一个小的设定：我们设立 $X_1$ 为处理变量，则 $X_2$ 为干扰变量，意味着我们可以单独估计所有干扰参数。

my = smf.ols("sales~temp+C(weekday)+cost", data=train).fit()
mt = smf.ols("price~temp+C(weekday)+cost", data=train).fit()
smf.ols("sales_res~price_res", 
        data=traindata.assign(sales_res=my.resid, # sales residuals
                              price_res=mt.resid) # price residuals
       ).fit().summary().tables[1]

换句话说，处理效应可以通过残差回归来得到，其中我们先将 $Y$ 对 $X$ 回归得到残差，再将残差回归于将 $T$ 对 $X$ 回归得到的残差。用符号表示为：$(Y - (Y \sim X)) \sim (T - (T \sim X))$，这样子我们最后求出来的系数就是我们想要的 $\tau$.

双重/去偏机器学习DML

在这个方法里，我们构建结果和处理残差时使用机器学习模型就可以很好的处理干扰变量。因为我们之前就提到了，我们的干扰变量很有可能不是线性关系的，而恰好，机器学习天然可以处理非线性效应和捕捉交互。

我们不必要在对干扰变量进行线性假设，在不存在未观测混杂的情况下，我们可以通过以下正交化步骤估计 ATE：

1. 使用机器学习模型 $M_y$ 用特征 $X$ 去预测结果 $Y$
2. 使用机器学习模型 $M_t$ 用特征 $X$ 去预测结果 $T$
3. 计算残差$r_y = Y - M_y(X)$ 和 $r_t = T - M_t(X)$
4. 将结果残差回归于处理残差 $r_y = \tau*r_x + e$

我们可以使用机器学习来更灵活地去拟合模型，但同样就面临着过拟合的问题，因此我们需要引入交叉预测和折外残差 (Out-of-fold Residuals)。

折外残差： 将数据集划分为 $K$ 折。模型在 $K-1$ 折上进行训练，并对未参与训练的那“一折”（Hold-out 折）进行预测。重复 $K$ 次，使每一折都有对应的预测值。OOF 残差就是这些基于未见过数据的预测值与真实标签的差值，其目的主要是为了防止过拟合。

估计ATE

我们首先从处理模型 $M_t$ 开始分析，我们依旧使用LightGBM框架：

from lightgbm import LGBMRegressor
from sklearn.model_selection import cross_val_predict

y = "sales"
T = "price"
X = ["temp", "weekday", "cost"]

debias_m = LGBMRegressor(max_depth=3, force_col_wise=True, verbose=-1)

train_pred = train.assign(price_res =  train[T] -
                          cross_val_predict(debias_m, train[X], train[T], cv=5) # K-fold 交叉验证
                          )
# 对于一些情况，如可视化，可以在最后加上均值，便于观察

在这里， $M_t$ 又被叫做去偏模型，此时处理这个变量已经转变为处理残差，因为其中来自 $X$ 的混杂偏差已经都被模型所剔除， 换句话来说，我们知道 $r_t$ 与 $X$ 正交，因为 $r_t$ 已经不能用 $X$ 进行解释。

同样地，我们开始结果模型 $M_y$：

denoise_m = LGBMRegressor(max_depth=3, force_col_wise=True, verbose=-1)

train_pred = train_pred.assign(sales_res =  train[y] -
                               cross_val_predict(denoise_m, train[X], train[y], cv=5)
                               )

这个模型的主要作用，是为了消除 $Y$ 的方差，也称为降噪模型 （denoising model）。直观来说，我们用 $r_y$ 去表示我们模型里的 $Y$，其中来自 $X$ 的所有方差都被解释掉了。

最后，为了估计ATE，我们对残差进行回归：

final_model = smf.ols(formula='sales_res ~ price_res', data=train_pred).fit()
final_model.summary().tables[1]

估计CATE （R-Learner）

现在，在CATE的分析里，我们知道上述的因果参数 $\tau$ 会随着单位协变量变化而变化。

为了估计这个模型，我们将使用同样的价格和销售残差，但现在我们将价格残差与其他协变量相互作用。

也就是说我们在建立去偏和降噪模型的步骤是一样的，只有最后对残差建立回归有区别。

应注意，为了预测 $CATE$， 我们将使用随机处理的测试集。

final_model_cate = smf.ols(formula='sales_res ~ price_res * (temp + C(weekday) + cost)', 	
                           data=train_pred).fit()

cate_test = test.assign(cate=final_model_cate.predict(test.assign(price_res=1))
                        - final_model_cate.predict(test.assign(price_res=0)))

非参数DML

无论上面我们估计ATE还是CATE，最后的残差回归我们都是用到线性回归去建模，这也意味着，我们这里有一个潜在的假设：假设价格对销售的影响是线性的。但真实世界的微观经济学里，我们知道其边缘获益很多时候并不是线性的。

一样的，我们可以将最后一步的对残差建模换成机器学习模型。让机器学习去捕捉非线性和复杂的函数形式。

与之前一样，我们先建立去偏模型和降噪模型;

y = "sales"
T = "price"
X = ["temp", "weekday", "cost"]

debias_m = LGBMRegressor(max_depth=3, force_col_wise=True, verbose=-1)
denoise_m = LGBMRegressor(max_depth=3, force_col_wise=True, verbose=-1)

train_pred = train.assign(price_res =  train[T] - cross_val_predict(debias_m, train[X], train[T], cv=5),
                          sales_res =  train[y] - cross_val_predict(denoise_m, train[X], train[y], cv=5))

我们在DML估计CATE里，最后我们估计的是：

我们注意到，这里有个误差项，换句话说，如果误差项最小，我们估计的 $\tau(X_i)$ 就是最为准确的，因此，我们要使得这个误差项最小，我们回到高中的时候，把感兴趣的变量放到一侧，进行分离变量：

这个我们跟训练神经网络得时候一样，这看着就像损失函数的样子，所以我们可以称作因果损失函数。如果我们最小化该损失的平方，我们将估计 $\tau(X_i)$ 的期望，即CATE。

这个损失也称为 R-Loss，现在我们对他进行最小化操作：

最小化上述损失等价于最小化括号内的表达式，我们用 $r_{Ti}^2$ 作为权重，用以构建非参数因果损失。然后我们就可以使用任意预测方法在权重下预测目标 ${r_{Yi}}/{r_{Ti}}$

model_final = LGBMRegressor(max_depth=3, force_col_wise=True, verbose=-1)
 
# 权重（处理残差的平方）
w = train_pred["price_res"] ** 2 
 
# 训练对象 （其实是在每一个点上的瞬时变化率）其实就是CATE
y_star = (train_pred["sales_res"] / train_pred["price_res"])
 
model_final.fit(X=train[X], y=y_star, sample_weight=w)

# 估计测试集的CATE
cate_test_non_param = test.assign(cate=model_final.predict(test[X]))

在这里，“非参数”主要指我们不预先规定 CATE 必须服从某个固定函数形式。对于价格这类连续处理变量，前面的 R-Loss 更像是在估计当前价格附近的边际斜率：

如果真实的 $\mu(t,x)$ 是弯曲的，那么某个价格点附近的斜率不一定能代表另一个价格区间。把低价区间估计出的边际效应用到高价区间，本质上仍然是在外推，可能高估也可能低估反事实销量。

那么，我们应该怎么办？一个更稳妥的思路是：让最终模型学习一条关于“处理残差”的非线性响应曲线，而不是只保留某一点上的局部斜率。

交叉估计

交叉估计（cross-fitting）的核心作用是给每个样本构造样本外残差。假设我们要估计两个干扰函数：

如果直接在全体训练集上拟合 $\hat m_y,\hat m_t$，再用同一批样本计算残差，机器学习模型很容易把训练样本记住。这样得到的残差会过于乐观，最后一步因果模型会把一部分过拟合误差误认为处理效应。交叉估计的做法是把训练集分成 $K$ 折。对于第 $k$ 折样本，只用其余 $K-1$ 折训练干扰模型，再回到第 $k$ 折上做预测，使用交叉估计是为了让残差化这一步尽量不被过拟合所影响：

因此每个样本的残差都是由“没见过它”的模型产生的：

这一步是 DML 里的关键。正交化负责降低干扰函数估计误差对目标参数的影响，而交叉估计则让这个正交化更可信，因此，我们最终训练的模型看到的是样本外残差，而不是训练集内的拟合残差。

from sklearn.model_selection import KFold

def cv_estimate(data, n_splits, model_cls, model_params, X, target):
    cv = KFold(n_splits=n_splits, shuffle=True, random_state=123)
    oof_pred = pd.Series(np.nan, index=data.index)
    models = []

    for train_idx, valid_idx in cv.split(data):
        model = model_cls(**model_params)
        model.fit(data[X].iloc[train_idx], data[target].iloc[train_idx])
        oof_pred.iloc[valid_idx] = model.predict(data[X].iloc[valid_idx])
        models.append(model)

    return oof_pred, models

接着我们分别对结果变量和处理变量做交叉预测，并构造残差：

y = "sales"
T = "price"
X = ["temp", "weekday", "cost"]

lgbm_params = dict(max_depth=3, force_col_wise=True, verbose=-1, random_state=123)

y_hat_oof, models_y = cv_estimate(train, 5, LGBMRegressor, lgbm_params, X, y)
t_hat_oof, models_t = cv_estimate(train, 5, LGBMRegressor, lgbm_params, X, T)

train_res = train.assign(
    y_res=train[y] - y_hat_oof,
    t_res=train[T] - t_hat_oof,
)

前面 R-Loss 的写法相当于假设：

也就是说，在给定 $X_i$ 之后，处理残差 $r_{Ti}$ 对结果残差 $r_{Yi}$ 的影响是局部线性的。非参数的设定可以把这个假设放松为：

其中 $g(\cdot)$ 可以是 LightGBM、随机森林、神经网络等任意预测模型。它学习的不是原始销量函数，而是“在排除 $X$ 能解释的部分之后，价格偏离其正常水平多少，会让销量偏离其正常水平多少”。如果把 $g$ 限制成 $g(x,r)=\tau(x)r$，即线性假设，就退回到了前面的 CATE/R-Learner 写法；如果允许 $g$ 随 $r$ 非线性变化，就可以刻画更弯曲的价格响应。

final_X = X + ["t_res"]
monotone_constraints = [0 for _ in X] + [-1]

model_final = LGBMRegressor(
    max_depth=3,
    monotone_constraints=monotone_constraints,
    force_col_wise=True,
    verbose=-1,
    random_state=123,
)

model_final.fit(train_res[final_X], train_res["y_res"])

这里设置 monotone_constraints=-1 的含义是：在其他变量不变时，价格残差越高，预测的销量残差应该下降。这个约束不是 DML 必需的，只是在价格和销量的例子里符合经济直觉，可以减少最终模型学出明显反常的曲线。

训练完成后，我们可以把某个样本的价格替换成一组候选价格，并预测这些价格下的反事实销量。设候选处理水平为 $t$，则预测公式是：

其中

也就是说，训练集内构造残差时使用 out-of-fold 预测；到了测试集或反事实网格上，因为这些样本本来就没有参与任何 nuisance 模型训练，所以可以把 $K$ 个模型的预测取平均，作为集成预测。

def ensemble_pred(df, models, X):
    return np.mean([m.predict(df[X]) for m in models], axis=0)

price_grid = np.linspace(train[T].quantile(0.05), train[T].quantile(0.95), 9)

pred_grid = (
    test
    .reset_index()
    .rename(columns={"index": "unit_id", T: "factual_price"})
    .assign(_key=1)
    .merge(pd.DataFrame({"_key": 1, T: price_grid}), on="_key")
    .drop(columns="_key")
)

pred_grid["t_res_cf"] = pred_grid[T] - ensemble_pred(pred_grid, models_t, X)
pred_grid["y_base"] = ensemble_pred(pred_grid, models_y, X)
pred_grid[f"{y}_pred"] = (
    pred_grid["y_base"]
    + model_final.predict(pred_grid[X].assign(t_res=pred_grid["t_res_cf"]))
)

pred_grid.loc[pred_grid["unit_id"].eq(0), ["unit_id", "factual_price", T, f"{y}_pred"]]

例如，对同一天、同样的天气和成本，我们可以只改变价格，得到一条反事实需求曲线：

这张表的重点不是某一个价格点的预测，而是整条曲线的形状：价格升高时销量是否下降，下降是否越来越快，还是在某个区间趋于平缓。若我们关心某个价格水平 $t$ 附近的边际效应，可以把它理解为响应曲线在该点的斜率：

如果最终模型是树模型，这个导数通常不会直接计算，而是用相邻网格点的差分近似：

我们通过生存多个处理水平，把模型学到的局部斜率扩展成一条可检查的响应曲线。这样我们可以发现模型是否给出了违反直觉的非线性形状，同时也可以避免只拿某一点的 CATE 去线性外推所有价格。

不过，这个方法并不能解决所有外推问题。反事实价格最好落在训练数据有足够支持的范围内，因此上面的例子使用了训练集价格的 5% 到 95% 分位数。如果把价格设到训练集中几乎没出现过的区域，$t-\bar{\hat m}_t(x)$ 也会落到最终模型不熟悉的残差范围里，此时预测仍然是在外推。