第4章 线性回归的不可思议有效性
在本章中，你将在因果推断武器库中增添首个重要的去偏技术：线性回归或普通最小二乘法\(OLS\)以及正交化。你将了解到，在估计处理变量与结果变量之间的关系时，线性回归如何对混杂因素进行调整。但不止于此，我希望让你掌握 处理正交化 这一强大概念。这个源于线性回归的理念，在你日后开始使用机器学习模型进行因果推断时会大有用处。
你所需的一切不过是线性回归
别急着翻到下一章，心里想着：“哦，回归太简单了！我成为数据科学家所学的第一个模型就是它” 之类的话。我向你保证，不，实际上你并不了解线性回归。事实上，回归是因果推断中最引人入胜且功能强大但又极具危险性的模型之一。诚然，它已经有一百多年的历史了。但直至今日，即便是最优秀的因果推断研究人员也常常会被它搞得措手不及。
普通最小二乘法研究
不信我？看看最近发表的一些关于这个主题的论文，你就明白了。一个不错的切入点是
\(Andrew \space Goodman-Bacon\)
撰写的《Difference-in-Difference with Variation in Treatment
Timing》或者 \(Tymon \space
Sloczynski\) 的《Interpreting OLS Estimands when Treatment
Effects Are Heterogeneous》，甚至还有 \(Goldsmith-Pinkham \space et \space al.\)
的《Contamination Bias in Linear Regressions》。
我向你保证：回归不仅是因果推断的主力军，而且会是你使用最为频繁的方法。回归也是更高级技术的主要基石，比如大多数面板数据方法—双重差分法和双向固定效应、机器学习方法—双重/去偏机器学习以及其他识别技术—工具变量法或断点回归设计。
为何我们需要模型
希望我已经成功说服你留下来，现在咱们言归正传。为了说明使用回归的必要性，让我们来思考银行业及整个借贷行业中一个颇具挑战性的问题：弄清楚贷款金额或信用卡额度对违约率的影响。很自然地，提高某人的信用卡额度会增加或者至少不会降低他们拖欠信用卡账单的几率。但是，如果你查看任何银行数据，就会发现信用额度与违约率之间呈负相关。显然，这并不意味着更高的额度会让客户违约的可能性降低。相反，这仅仅反映了处理分配机制：银行和贷款公司会依据其承保模型将更多信贷提供给违约可能性较低的客户。你所看到的这种负相关，其实是混杂偏差的影响：
import warnings
warnings.filterwarnings('ignore')
# 用于绘制因果图、流程图等结构化图形
import graphviz as gr
from matplotlib import style
import seaborn as sns
from matplotlib import pyplot as plt
style.use("ggplot")
import statsmodels.formula.api as smf
import matplotlib
# 用于定义 matplotlib 中的样式循环（如颜色、线条的循环规则）
from cycler import cycler
default_cycler = (cycler(color=['0.3', '0.5', '0.7', '0.5']))
color=['0.3', '0.5', '0.7', '0.9']
linestyle=['-', '--', ':', '-.']
marker=['o', 'v', 'd', 'p']
plt.rc('axes', prop_cycle=default_cycler)
gr.set_default_format("png")## 'pdf'
g_risk = gr.Digraph(graph_attr={"rankdir":"LR"})
g_risk.edge("Risk", "X")
g_risk.edge("X", "Credit Limit")
g_risk.edge("X", "Default")
g_risk.edge("Credit Limit", "Default")
# g_risk.render(filename='risk', format='png', view=False)
当然，银行并不了解违约的内在风险，但它可以使用代理变量\(X\)（比如收入或信用评分）来对其进行估算。在前几章中，你已经了解了如何对变量进行调整，以使处理干预看起来像是随机分配的。具体来说，你已经了解了调整公式：\(ATE = E_X \{ E[Y \mid T = 1, X = x] - E[Y \mid T = 0, X = x] \}\) 。结合条件独立假设 \((Y_0, Y_1) \perp T \mid X\)，这个公式能让你识别出因果效应。
然而，如果你真的要照本宣科地应用这个调整公式，情况很快就会变得难以掌控。首先，你需要根据特征\(X\)将数据划分成不同的分段。要是你的离散特征非常少，这样做倒也没问题。但要是特征数量众多，而且其中一些还是连续的，那该怎么办呢？比如说，假设你知道银行用了\(10\)个变量，每个变量有\(3\)个分组，以此来审核客户并分配信用额度。这看起来好像不算多是吧？可实际上这会产生 \(59049\) 个，也就是 \(3^{10}\) 个单元格。要在每个单元格中估算\(ATE\)，然后对结果求平均，只有在你拥有海量数据的情况下才有可能实现。这就是维度诅咒，大多数数据科学家对这个问题都很熟悉。在因果推断的情境中，这种诅咒带来的一个影响是如果你有很多协变量，那么单纯套用调整公式就会遭遇数据稀疏性问题。
因果推断术语与机器学习术语
大多数数据科学家所熟知的机器学习文献和通常源自计量经济学或流行病学的因果推断文献使用的术语有所不同。所以，要是你需要在两者之间进行转换，以下是你会遇到的一些主要对应关系：
特征：Feature 协变量：Covariates 或自变量：independent variables
权重：Weights 参数：Parameters 或系数：coefficients
目标：Target 结果：Outcome 或因变量：dependent variable
解决这种维度问题的一种方法是假设潜在结果可以用类似线性回归的模型来建模，线性回归能够对由众多不同\(X\)定义的单元格进行插值和外推。在这种情境下，你可以把线性回归看作一种降维算法。它将结果变量投影到\(X\)变量上，并在这个投影上对处理组和对照组进行比较。这相当巧妙。不过我有点操之过急了。要真正（我是说打从心底里）理解回归，你得从小处入手：\(A/B\) 测试情境中的回归。
A/B 测试中的回归分析
假设你在一家在线流媒体公司工作负责优化其推荐系统。你的团队刚刚完成了该系统的新版本开发，采用了尖端技术和机器学习领域的最新理念。尽管这些成果令人印象深刻，但你的经理真正关心的是：这个新版本系统能否提高流媒体服务的观看时长？为了验证这一点，你决定开展一项\(A/B\)测试，具体步骤如下：
首先，从你的客户群体中抽取一部分具有代表性的样本（样本规模较小）；
然后，将新版本推荐系统随机部署到该样本中 \(\frac{1}{3}\) 的用户群体，其余用户则继续使用旧版本推荐系统；
\(1\) 个月后，你收集了测试结果，结果以日均观看时长为衡量指标：
import pandas as pd
import numpy as np
data = pd.read_csv("./data/rec_ab_test.csv")
data.head()## recommender age tenure watch_time
## 0 challenger 15 1 2.39
## 1 challenger 27 1 2.32
## 2 benchmark 17 0 2.74
## 3 benchmark 34 1 1.92
## 4 benchmark 14 1 2.47
由于推荐系统版本是随机分配的，对不同版本间的平均观看时长进行简单比较，就能得出 \(ATE\)。但之后为了检验统计显著性，你还得费尽周折地计算标准误以得到置信区间。那么，要是我告诉你，你可以用回归分析来解读\(A/B\)测试的结果，而且这样还能免费获得你所需的所有推断统计量，你觉得怎样呢？回归分析背后的思路是，你要估算以下方程或模型：\(WatchTime_i = \beta_0 + \beta_1 challenger_i + e_i\) 。其中，若客户属于使用新版本推荐系统的组，\(challenger\) 取值为 \(1\)，否则为 \(0\)。要是你估算这个模型，新版本的影响会由 \(\beta_1\) 的估计值 \(\widehat{\beta_1}\) 体现。
要在 \(Python\) 中运行这个回归模型，你可以使用 \(statsmodels\) 的公式应用程序接口 。它能让你用 \(R\) 风格的公式简洁地表示线性模型。比如，你可以用公式 \(watch\_time \sim C(recommender)\) 来表示前面的模型。要估算该模型只需调用 \(.fit()\) 方法，要查看结果就对之前拟合好的模型调用 \(.summary()\) 方法 ：
import statsmodels.formula.api as smf
result = smf.ols('watch_time ~ C(recommender)', data=data).fit()
print(result.summary().tables[1])## ================================================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------------------------
## Intercept 2.0491 0.058 35.367 0.000 1.935 2.163
## C(recommender)[T.challenger] 0.1427 0.095 1.501 0.134 -0.044 0.330
## ================================================================================================
在那种\(R\)风格的公式里，结果变量放在前面，接着是波浪号 \(\sim\)。然后，加上解释变量。在这种情况下，你只会用到推荐系统变量，它是有两个类别的分类变量。你可以用 \(C(...)\) 把那个变量括起来，明确表明该列是分类变量。
Patsy
这种公式语法糖是进行特征工程的一种极其便捷的方式。你可以在 \(pasty\) 库中了解更多相关内容 。
接下来查看结果。首先有截距项。它是模型中 \(\beta_0\) 参数的估计值。它告诉你，当模型中其他变量为 \(0\) 时结果变量的预期值。由于这里唯一的其他变量是 \(challenger\) 指示变量，你可以将截距项解读为使用旧版推荐系统的用户的预期观看时长。在这里，这意味着使用旧版推荐系统时，客户平均每天花费 \(2.04\) 小时观看你们的流媒体内容。
最后，看与新版本推荐系统相关的参数估计值 \(\widehat{\beta_1}\)，你可以了解到因这个新版本而带来的观看时长增加量。如果 \(\widehat{\beta_0}\) 是旧版推荐系统下观看时长的估计值，那么 \(\widehat{\beta_0} + \widehat{\beta_1}\) 能告诉你使用新版本推荐系统的用户的预期观看时长。换句话说，\(\widehat{\beta_1}\) 是 \(ATE\) 的一个估计值。由于是随机分配（进行 \(A/B\) 测试时的随机分组 ），你可以给这个估计值赋予因果意义：你可以说，新推荐系统平均每天将观看时长增加了 \(0.14\) 小时。不过，这个结果在统计上并不显著。
先别管这个不显著的结果，因为你刚刚完成的事相当了不起。你不仅估计出了 \(ATE\)，还免费得到了置信区间和 \(p\) 值！不仅如此，你自己也能明白，回归分析正在做它该做的事 —— 为每种处理估计 \(E[Y \mid T]\)（给定处理下结果的期望 ） ：
data.groupby("recommender")["watch_time"].mean()## recommender
## benchmark 2.049064
## challenger 2.191750
## Name: watch_time, dtype: float64就像我所说的，从数学层面来讲，截距项等同于对照组的平均观看时长。这些数值是一样的，因为在这种情况下，回归分析在数学上等同于直接对平均值进行比较。这也意味着 \(\widehat{\beta_1}\) 是两组之间的平均差值：\(2.191 - 2.049 = 0.1427\)
好啦，你确实实实在在地用回归分析重现了分组平均值。但那又怎样呢？又不是说你之前没法做这种简单比较，所以这么做真正的好处到底是什么呢？
用回归进行调整
为了体会回归的强大作用，让我把你带回到最初的例子：估算信用额度对违约的影响。银行数据通常是这样的，有一堆列呈现客户特征，这些特征可能预示着信用可靠性，比如月工资、信用机构提供的各类信用评分、在当前公司的工作年限等等。然后，有给该客户的信用额度（在这个例子里就是处理干预），还有一列告诉你客户是否违约 —— 也就是结果变量：
risk_data = pd.read_csv("./data/risk_data.csv")
risk_data.head()## wage educ exper ... credit_score2 credit_limit default
## 0 950.0 11 16 ... 518.0 3200.0 0
## 1 780.0 11 7 ... 429.0 1700.0 0
## 2 1230.0 14 9 ... 571.0 4200.0 0
## 3 1040.0 15 8 ... 411.0 1500.0 0
## 4 1000.0 16 1 ... 518.0 1800.0 0
##
## [5 rows x 8 columns]模拟数据
我再次以真实世界的数据为基础进行构建，并对其加以修改以满足本章的需求。这次，我使用的是 \(wage1\) 数据集，该数据集由 \(Jeffrey \space M.Wooldridge\) 教授整理，可在 \(woolsridge\) 包中获取 。
这里，处理变量 \(credit\_limit\) 的类别数量过多。在这种情况下将其视为连续变量而非分类变量会更妥当。你可以把平均处理效应表示为结果的期望对处理变量的导数而非处理变量多个水平之间的差值：\(ATE = \frac{\partial}{\partial t} E [y \mid t]\)
如果这听起来很复杂，别担心。它的意思很简单，就是在处理变量增加一个单位时，你预期结果会发生的变化量。在这个例子中，它代表当信用额度增加 \(1\) 美元时，你预期违约率会发生的变化程度 。
估计这类变量的一种方法是进行回归分析。具体来说，你可以估计以下模型：\(Default_i = \beta_0 + \beta_1 limit_i + e_i\)。其中，估计值 \(\widehat{\beta_1}\) 可解读为：当信用额度增加 \(1\) 美元时，你预期风险发生的变化量。要是信用额度是随机分配的，这个参数就具有因果解释意义。但你很清楚，实际情况并非如此，因为银行往往会给风险较低的客户更高的信用额度。实际上，要是你运行上述模型，你会得到 \(\beta_1\) 的负估计值：
model = smf.ols('default ~ credit_limit', data=risk_data).fit()
print(model.summary().tables[1])## ================================================================================
## coef std err t P>|t| [0.025 0.975]
## --------------------------------------------------------------------------------
## Intercept 0.2192 0.004 59.715 0.000 0.212 0.226
## credit_limit -2.402e-05 1.16e-06 -20.689 0.000 -2.63e-05 -2.17e-05
## ================================================================================鉴于风险和信用额度之间的关系因混杂因素而呈负相关，这一点也不奇怪。要是你把拟合的回归直线和按信用额度划分的平均违约率画在一起，就能清楚地看到这种负向趋势：
plt_df = (risk_data
.assign(size=1)
.groupby("credit_limit")
.agg({"default": "mean", "size":sum})
.reset_index()
.assign(prediction = lambda d: model.predict(d)))
plt.figure(figsize=(10,4))
sns.scatterplot(data=plt_df,
x="credit_limit",
y="default",
size="size",
sizes=(1,100))
plt.plot(plt_df["credit_limit"], plt_df["prediction"], color="C1")
plt.title("Default Rate by Credit Limit")
# plt.show()
为了调整这种偏差，理论上，你可以按所有混杂因素对数据进行分段，在每个分段内对违约情况关于信用额度进行回归分析，提取斜率参数，然后对结果求平均。然而，由于维度诅咒，即使你尝试针对数量不算多的混杂因素（比如信用评分）这样做，你也会发现有些单元格里只有一个样本，这使得你无法进行回归分析。更不用说还有许多单元格完全是空的。
risk_data.groupby(["credit_score1", "credit_score2"]).size().head()## credit_score1 credit_score2
## 34.0 339.0 1
## 500.0 1
## 52.0 518.0 1
## 69.0 214.0 1
## 357.0 1
## dtype: int64幸好再一次回归在这里能帮上你的忙。你无需手动对混杂因素进行调整，只需把它们添加到你要用\(OLS\)估计的模型里就行：\(Default_{i}=\beta_{0}+\beta_{1}limit_{i}+\theta X_{i}+\epsilon_{i},\) 这里\(\boldsymbol{X}\) 是混杂变量的向量，\(\theta\) 是与这些混杂因素相关的参数向量。\(\theta\) 参数没什么特别的，它们的表现和 \(\beta_{1}\) 完全一样。我把它们表示得不一样，是因为它们只是用来帮你得到对 \(\beta_{1}\) 的无偏估计。也就是说，你其实并不关心它们的因果解释（严格来说，它们被称为干扰参数）。
在信用的例子中，你可以把信用分数和工资这些混杂因素加到模型里。模型会是这样的： \(Default_{i}=\beta_{0}+\beta_{1}limit_{i}+\theta_{1}wage_{i}+\theta_{2}creditScore1_{i}+\theta_{3}creditScore2_{i}+\epsilon_{i}.\) 我会更详细地讲模型中纳入变量是如何对混杂因素进行调整的，但现在有个很简单的方法能看明白。前面的模型是关于 \(E[y|t,X]\) 的模型。回想一下，你想要的是 \(\frac{\partial}{\partial t}E[y|t,X]\)。要是你对处理变量—信用额度求这个模型的偏导，会发生什么呢？嗯，你直接得到 \(\beta_{1}\)！从某种意义上说，\(\beta_{1}\)可以看作是违约的期望值对信用额度的偏导数。或者更直观地说，它可以被看作是：在模型中所有其他变量都保持不变的情况下，信用额度有一个小幅度增加时，你预计违约会有多大变化。这种解释已经能让你了解一点回归是如何对混杂因素进行调整的：在估计处理变量和结果变量之间的关系时，它会将混杂因素固定住。
要实际验证这一点，你可以估计上述模型。只需加入一些混杂因素，奇妙的是信用额度和违约之间的关系就会变成正向的！
formula = 'default ~ credit_limit + wage+credit_score1+credit_score2'
model = smf.ols(formula, data=risk_data).fit()
print(model.summary().tables[1])## =================================================================================
## coef std err t P>|t| [0.025 0.975]
## ---------------------------------------------------------------------------------
## Intercept 0.4037 0.009 46.939 0.000 0.387 0.421
## credit_limit 3.063e-06 1.54e-06 1.987 0.047 4.16e-08 6.08e-06
## wage -8.822e-05 6.07e-06 -14.541 0.000 -0.000 -7.63e-05
## credit_score1 -4.175e-05 1.83e-05 -2.278 0.023 -7.77e-05 -5.82e-06
## credit_score2 -0.0003 1.52e-05 -20.055 0.000 -0.000 -0.000
## =================================================================================别被 \(\beta_{1}\) 的小估计值骗了。回想一下，\(cresut\_limit\) 的单位是千美元，而 \(default\) 是 \(0\) 要么是 \(1\)。所以，信用额度每增加 \(1\) 美元，预期违约率只会有很小的增长，这并不奇怪。不过这个数值在统计上是显著的，它告诉你，随着信用额度的提高风险也会增加，这更符合你对世界运行方式的直觉。
先记住这个想法，因为你马上要更正式地探究它了。终于到了学习所有因果推断工具中最强大的工具之一的时候了：\(FWL\) 定理。这是一种消除偏差的绝妙方法，可惜数据科学家们大多不了解它。\(FWL\) 是理解更先进去偏方法的前提，但我觉得它最有用的地方在于，它可以作为一个去偏的预处理步骤。还是用银行业的例子来说，假设这家银行里很多数据科学家和分析师都在试图弄清楚信用额度如何影响很多不同的业务指标而不只是风险。但只有你知道信用额度是如何分配的，这意味着只有你清楚信用额度这个处理变量受到了什么样的偏差影响。借助 \(FWL\)，你可以利用这些知识对信用额度数据进行去偏，这样不管其他人对哪个结果变量感兴趣，都能使用这些数据。\(FWL\)定理能让你把去偏步骤和影响估计步骤分离开来。但要学习它，你得先快速复习一点回归理论。
回归理论
我不打算深入探讨线性回归是如何构建和估计的。不过，一点理论知识对解释它在因果推断中的强大作用会大有帮助。首先，回归解决了最佳线性预测问题。设\(\beta^*\)为参数向量：\(\beta^* = \underset{\beta}{argmin} \ E\left[(Y_i -
X_i'\beta)^2\right]\)
线性回归会找到使均方误差最小的参数。如果你对其求导并令导数为 \(0\)，就会发现这个问题的线性解由下式给出：\(\beta^* =
E[X'X]^{-1}E[X'Y]\)。你可以用样本等效式来估计这个 \(\beta\)：\(\widehat{\beta} =
(X'X)^{-1}X'Y\)
但别只听我这么说。要是你属于那种更懂代码而非公式的人，自己试试吧。在下面的代码中，我用 \(OLS\) 的代数解法来估计你刚看到的模型的参数（我把截距作为最后一个变量加入，所以第一个参数估计值会是 \(\widehat{\beta_1}\)）：
X_cols = ["credit_limit", "wage", "credit_score1", "credit_score2"]
X = risk_data[X_cols].assign(intercep=1)
y = risk_data["default"]
def regress(y, X):
return np.linalg.inv(X.T.dot(X)).dot(X.T.dot(y))
beta = regress(y, X)
beta## array([ 3.06252773e-06, -8.82159125e-05, -4.17472814e-05, -3.03928359e-04,
## 4.03661277e-01])如果你往回看一看，就会发现这些数字和你之前用 \(statsmodels\) 里的 \(ols\) 函数估计模型时得到的数字完全一样。
我经常使用 \(pandas\) 中的 \(.assgin()\) 方法。如果你不熟悉这个方法，简单来说，它会返回一个新的 \(DataFrame\)，并且这个 \(DataFrame\) 中包含你传递给该方法的新创建的列。
# new_df = df.assign(new_col_1 = 1, new_col_2 = df["old_col"] + 1)
# new_df[["old_col", "new_col_1", "new_col_2"]].head()单变量线性回归
上一节中的 \(\widehat{\beta}\) 公式具有相当的普遍性。不过，研究只有一个回归变量的情况是很有价值的。在因果推断中，你常常想要估计变量 \(T\) 对结果 \(Y\) 的因果影响。所以，你会用这个单变量的回归来估计这种影响。
对于单个回归变量 \(T\)，与之相关的参数由下式给出： \(\widehat{\tau} = \frac{Cov(Y_i, T_i)}{Var(T_i)} = \frac{E\left[(T_i - \overline{T})\left(Y_i - \overline{Y}\right)\right]}{E\left[(T_i - \overline{T})^2\right]}\)
如果 \(T\) 是随机分配的，\(\beta_1\)就是平均处理效应。重要的是，通过这个简单的公式，你能明白回归在做什么：它在找出处理变量和结果变量是如何共同变化的（如分子中的协方差所表达的），并通过处理变量的单位来缩放这种共同变化，这是通过除以处理变量的方差来实现的。
注意
你也可以把这和通用公式联系起来。协方差与点积密切相关。所以你几乎可以说，在协方差/方差公式中，\(X'X\)起到了分母的作用，而\(X'y\)起到了分子的作用。
多元线性回归
事实证明，除了您之前看到的一般公式外，还有另一种看待多元线性回归的方式。这种另一种方式让我们更清楚回归在做什么。
如果您有不止一个回归变量，您可以扩展单变量回归公式以适应这种情况。假设其他变量只是辅助性的，而您真正感兴趣的是估计与 \(T\) 相关的参数 \(\tau\)： \(y_i = \beta_0 + \tau T_i + \beta_1 X_{1i} + \dots + \beta_k X_{ki} + u_i\), \(\tau\) 可以用以下公式估计： \(\widehat{\tau} = \frac{Cov(Y_i, \widetilde{T}_i)}{Var(\widetilde{T}_i)}\)。其中 \(\widetilde{T}_i\) 是 \(T_i\) 对所有其他协变量 \(X_{1i} + \dots + X_{ki}\) 进行回归得到的残差。
现在让我们来体会一下这有多酷。这意味着，多元回归的系数是在考虑了模型中其他变量的影响后同一回归变量的二元系数。在因果推断术语中，\(\tau\) 是在使用所有其他变量对 \(T\) 进行预测后 \(T\) 的二元系数。
这背后有一个很好的直觉。如果您可以使用其他变量来预测 \(T\)，这意味着它不是随机的。然而，一旦您控制了所有的混杂变量 \(X\)，您可以让 \(T\)看起来尽可能随机。要做到这一点，您可以使用线性回归从混杂因素中预测它，然后取该回归的残差 \(\widetilde{T}\)。根据定义，\(\widetilde{T}\) 无法被您已经用来预测 \(T\) 的其他变量 \(X\) 预测。非常巧妙的是，\(\widetilde{T}\) 是处理变量的一个版本，它与 \(X\) 中的任何其他变量都不相关。
我知道这有点拗口，但这真的很棒。事实上，这已经是我答应要教给你们的 \(FWL\) 定理的内容了。所以，如果您不太理解多元回归的这部分内容，也不用担心，因为您即将以一种更直观、更形象的方式重新学习它。
FWL定理与正交化
\(FWL\)正交化是你可运用的首个主要去偏技术。它是一种简单却强大的方法，能让非实验数据看起来仿佛处理是随机化的。\(FWL\) 主要与线性回归相关；正如你将在第 \(3\) 部分看到的，\(FWL\) 正交化已被扩展到更一般的情境中。\(FWL\) 定理指出，多元线性回归模型可以一次性全部估计，也可以分 \(3\) 个独立步骤估计。例如，你可以像之前已经做过的那样，对 \(default\) 关于 \(credit\_limit\)、\(wage\)、\(credit\_score1\)、\(credit\_score2\) 进行回归。
formula = 'default ~ credit_limit + wage+credit_score1+credit_score2'
model = smf.ols(formula, data=risk_data).fit()
print(model.summary().tables[1])## =================================================================================
## coef std err t P>|t| [0.025 0.975]
## ---------------------------------------------------------------------------------
## Intercept 0.4037 0.009 46.939 0.000 0.387 0.421
## credit_limit 3.063e-06 1.54e-06 1.987 0.047 4.16e-08 6.08e-06
## wage -8.822e-05 6.07e-06 -14.541 0.000 -0.000 -7.63e-05
## credit_score1 -4.175e-05 1.83e-05 -2.278 0.023 -7.77e-05 -5.82e-06
## credit_score2 -0.0003 1.52e-05 -20.055 0.000 -0.000 -0.000
## =================================================================================但根据\(FWL\)，你也可以将这种估计分解为：
去偏步骤：将处理变量 \(T\) 对混杂因素 \(X\) 进行回归，得到处理残差 \(\widetilde{T} = T - \widehat{T}\)。
去噪步骤：将结果变量 \(Y\) 对混杂变量 \(X\) 进行回归，得到结果残差 \(\widetilde{Y} = Y - \widehat{Y}\)。
结果模型：将结果残差 \(\widetilde{Y}\) 对处理残差 \(\widetilde{T}\) 进行回归，以得到 \(T\) 对 \(Y\) 因果效应的估计。
不出所料，这只是你刚才在《多元线性回归》中看到的公式的重述。\(FWL\) 定理阐述了回归模型在估计程序上的等价性。它还表明，你可以分离出线性回归的去偏部分，也就是前面列表中概述的第一步。
为了更好地理解到底是怎么回事，让我们逐步分析。
去偏步骤
回想一下，最初由于混杂偏差，你的数据看起来是这样的，\(default\) 随着 \(credit\_limit\) 的增加而呈下降趋势。
根据 \(FWL\) 定理，你可以通过拟合一个回归模型，用混杂因素来预测处理变量—信用额度，从而对这些数据进行去偏。然后，你可以从这个模型中获取残差：\(\widetilde{line}_i = line_i - \widehat{line}_i\)。这个残差可以被视为处理变量的一个版本，它与去偏模型中使用的变量不相关。这是因为根据定义，残差与生成预测的变量是正交的。
这个过程会使 \(\widetilde{line}\) 以 \(0\) 为中心。另外，你可以把平均处理量 \(\overline{line}\) 加回来：\(\widetilde{line}_i = line_i - \widehat{line}_i + \overline{line}\)
这对去偏来说不是必需的，但它能让 \(\widetilde{line}\) 处于与原始 \(line\) 相同的范围，这对于可视化来说更好：
debiasing_model = smf.ols('credit_limit ~ wage + credit_score1 + credit_score2', data=risk_data).fit()
risk_data_deb = risk_data.assign(credit_limit_res=(debiasing_model.resid + risk_data["credit_limit"].mean()))
risk_data_deb.head(5)## wage educ exper ... credit_limit default credit_limit_res
## 0 950.0 11 16 ... 3200.0 0 3347.583906
## 1 780.0 11 7 ... 1700.0 0 2331.086567
## 2 1230.0 14 9 ... 4200.0 0 3665.134174
## 3 1040.0 15 8 ... 1500.0 0 1720.933260
## 4 1000.0 16 1 ... 1800.0 0 2084.659827
##
## [5 rows x 9 columns]如果你现在进行简单线性回归，将结果变量对去偏后的或残差化的处理变量进行回归，你就已经能得到在控制了去偏模型中所用混杂变量的情况下，信用额度对风险的影响。此处你得到的参数估计值，与你之前前提到的通过运行完整模型（即同时纳入处理变量和混杂变量的模型）所得到的参数估计值完全相同。
model_w_deb_data = smf.ols('default ~ credit_limit_res', data=risk_data_deb).fit()
print(model_w_deb_data.summary().tables[1])## ====================================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------------
## Intercept 0.1421 0.005 30.001 0.000 0.133 0.151
## credit_limit_res 3.063e-06 1.56e-06 1.957 0.050 -4.29e-09 6.13e-06
## ====================================================================================替代系数公式
只需对处理变量进行残差化这一事实，暗示了一种重写回归系数公式的更简单方法。在单变量情形下，你可以使用：\(\beta_1 = \frac{E\left[(T_i - \overline{T})y_i\right]}{E\left[(T_i - \overline{T})^2\right]}\) 而非用 \(Y\) 和 \(T\) 的协方差除以 \(T\) 的方差。在多元情形下，公式为：\(\beta_1 = \frac{E\left[(T_i - E[T|X])y_i\right]}{E\left[Var(T|X)\right]}\)
不过，还是有区别的。看一下 \(p\) 值，它比你之前得到的略高。那是因为你没有进行去噪步骤而去噪步骤是负责降低方差的。尽管如此，仅通过去偏步骤，鉴于所有混杂因素都已纳入去偏模型，你已经能得到信用额度对风险因果影响的无偏估计。
你也可以通过绘制去偏版本的信用额度与违约率的关系图来直观地了解情况。你会发现，这种关系不再像数据存在偏差时那样呈向下倾斜的状态：
去噪步骤
虽然去偏步骤对于估计正确的因果效应至关重要，但去噪步骤也很不错，尽管没那么重要。它不会改变你的处理效应估计值，但会降低其方差。在这一步，你要将结果变量对非处理的协变量进行回归。然后，你会得到结果的残差 \(\widetilde{default}_i = default_i - \widehat{default}_i\)。
同样，为了更好地可视化，你可以把平均违约率加到去噪后的违约变量上，以便更好地呈现：\(\widetilde{default}_i = default_i - \widehat{default}_i + \overline{default}\)
denoising_model = smf.ols('default ~ wage + credit_score1 + credit_score2', data=risk_data_deb).fit()
risk_data_denoise = risk_data_deb.assign(default_res=denoising_model.resid + risk_data_deb["default"].mean())
risk_data_denoise.head()## wage educ exper ... default credit_limit_res default_res
## 0 950.0 11 16 ... 0 3347.583906 0.000981
## 1 780.0 11 7 ... 0 2331.086567 -0.043175
## 2 1230.0 14 9 ... 0 3665.134174 0.043289
## 3 1040.0 15 8 ... 0 1720.933260 -0.028427
## 4 1000.0 16 1 ... 0 2084.659827 0.000760
##
## [5 rows x 10 columns]回归估计量的标准误
既然我们在谈论噪声，我认为现在是了解如何计算回归标准误的好时机。回归参数估计值的标准误由以下公式给出： \(SE\left( \widehat{\beta} \right) = \frac{\sigma \left( \hat{\epsilon} \right)}{\sigma \left( \widetilde{T} \right) \sqrt{n - df}}\) 其中，\(\hat{\epsilon}\) 是回归模型的残差，\(df\) 是模型的自由度（模型估计的参数数量）。如果你更愿意通过代码来看这个公式的实现，如下所示：
model_se = smf.ols('default ~ wage + credit_score1 + credit_score2',data=risk_data).fit()
print("SE regression:", model_se.bse["wage"])## SE regression: 5.364242347548206e-06model_wage_aux = smf.ols('wage ~ credit_score1 + credit_score2',data=risk_data).fit()
# subtract the degrees of freedom - 4 model parameters - from N.
se_formula = np.std(model_se.resid)/(np.std(model_wage_aux.resid)*np.sqrt(len(risk_data)-4))
print("SE formula: ", se_formula)## SE formula: 5.364242347548193e-06这个公式很不错，因为它能让你进一步直观理解一般意义上的回归，尤其是去噪步骤。首先，分子表明你对结果的预测越好，残差就会越小。因此，估计量的方差也就越低。这正是去噪步骤的关键所在。它还表明如果处理变量能在很大程度上解释结果变量，其参数估计值的标准误也会更小。
有趣的是，误差还与残差化的处理变量的方差成反比。这也很直观。如果处理变量变化很大，就会更容易衡量它的影响。你会在《噪声诱导控制》中了解到更多相关内容。
连续处理的实验
如果您计划设计一个实验且希望将回归的参数估计值作为效应的衡量标准，那么标准误公式也会很有用。如果您想要随机化的处理变量是连续的，这会是个好主意。在这种情况下标准误公式可近似为：\(SE \approx \frac{\sigma(y)}{\sigma(T)\sqrt{n - 2}}\) 。在单变量回归模型的情况下这种近似是保守的，因为\(\sigma(y) \geq \sigma(\hat{\epsilon})\)，处理可能会在一定程度上解释结果。然后，您可以采用这个标准误并代入第 \(2\) 章中的样本量计算公式。重要的是，设计这个检验还有额外的复杂性，即要从 \(T\) 中选择抽样分布，这也会通过 \(\sigma(T)\) 影响标准误。
最终结果模型
有了残差 \(\widetilde{Y}\) 和 \(\widetilde{T}\)，你就可以进行 \(FWL\) 定理所概述的最后一步 —— 只需将 \(\widetilde{Y}\) 对 \(\widetilde{T}\) 进行回归即可。
model_w_orthogonal = smf.ols('default_res ~ credit_limit_res',data=risk_data_denoise).fit()
print(model_w_orthogonal.summary().tables[1])## ====================================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------------
## Intercept 0.1421 0.005 30.458 0.000 0.133 0.151
## credit_limit_res 3.063e-06 1.54e-06 1.987 0.047 4.17e-08 6.08e-06
## ====================================================================================处理变量的参数估计值与你在去偏步骤以及运行包含信用额度和所有其他协变量的回归模型时得到的参数估计值完全相同。此外，标准误和 \(p\) 值现在也和你第一次运行包含所有变量的模型时一样。这就是去噪步骤的效果。
当然，你也可以绘制出去偏后的处理变量与去噪后的结果变量之间的关系，以及最终模型的预测结果，来看看具体情况是怎样的：
FWL 总结
我不知道你是否已经看出来了，但我真的很喜欢说明性的图表。即使它们没有反映任何真实数据，也能非常有用地将一些相当技术性的概念背后的情况可视化。对于 \(FWL\) 来说也不例外。所以，总结一下，假设你想要估计处理变量 \(T\) 和结果变量 \(Y\) 之间的关系，但存在一些混杂因素 \(X\)。你把处理变量放在 \(x\) 轴上，结果变量放在 \(y\) 轴上，混杂因素作为颜色维度。起初，你看到处理变量和结果变量之间呈负斜率，但你有充分的理由（一些领域知识）相信这种关系应该是正的，所以你决定对数据进行去偏。
要做到这一点，你首先使用线性回归估计 \(E[T|X]\)。然后，构建处理变量的去偏版本：\(T - E[T|X]\)（见\(图4-1\)）。有了这个去偏后的处理变量，你已经能看到你希望找到的正相关关系了。但仍然存在很多噪声。
为了处理噪声，你同样使用回归模型来估计 \(E[Y|X]\)。然后，构建结果变量的去噪版本：\(Y - E[Y|X]\)（见\(图4-2\)）。你可以将这个去噪后的结果变量视为在解释了其中所有由 \(X\) 解释的方差之后的结果变量。如果 \(X\) 解释了 \(Y\) 中很大一部分方差，那么去噪后的结果变量的噪声会更少，从而更容易看到你真正关心的关系，即 \(T\) 和 \(Y\) 之间的关系。
最后，在完成去偏和去噪之后，你可以清楚地看到 \(T\) 和 \(Y\) 之间存在正相关关系。剩下要做的就是对这些数据拟合一个最终模型：
这个最终的回归模型将具有与你同时对 \(Y\) 关于 \(T\) 和 \(X\) 进行回归时完全相同的斜率。
去偏与截距
不过，有一点需要注意。在因果推断中，你主要关注的是这个回归的斜率，因为斜率是对连续处理效应 \(\frac{\partial E[Y|t]}{\partial t}\) 的线性近似。但如果你也关心截距—例如，如果你试图进行反事实预测 —— 你应该知道去偏和去噪会使截距等于零。
作为结果模型的回归
在本节中，我一直强调回归主要是通过对处理变量进行正交化来发挥作用的。不过，你也可以将回归视为一种潜在结果插补技术。假设处理变量是二元的。如果在控制组 \(T = 0\) 中，\(Y\) 对 \(X\) 的回归能很好地近似 \(E[Y_0|X]\)，那么你可以使用该模型来插补 \(Y_0\) 并估计 \(ATT\)： \(ATT = \frac{1}{N_1} \sum \mathbb{1}(T_i = 1) (Y_i - \widehat{\mu}_0(X_i))\) 其中 \(N_1\) 是处理组单元的数量。
指示函数
在整本书中，我会用 \(\mathbb{1}(.)\) 来表示指示函数。当括号内的参数为真时，该函数返回 \(1\)，否则返回 \(0\)。
类似地，可以说明：如果对处理组单元的回归能够建模 \(E[Y_1|X]\)，你可以用它来估计未处理组的平均效应。如果将这两个论点放在一起，你可以按如下方式估计 \(ATE\)： \(ATE = \frac{1}{N} \sum (\widehat{\mu}_1(X_i) - \widehat{\mu}_0(X_i))\)。这个估计量会为所有单元插补两个潜在结果。它等同于对 \(Y\) 关于 \(X\) 和 \(T\) 进行回归，并读取 \(T\) 的参数估计值。
或者，你可以只估算缺失的潜在结果：\(ATE=\frac{1}{N} \sum(\mathbb{1}(T_i=1)[Y_i - \hat\mu_0(X_i)] + \mathbb{1}(T_i=0)[\hat{\mu_1}(X_i)-Y_i])\) 。当 \(T\) 是连续的，这有点难以概念化，但你可以理解回归是估算整个处理响应函数，这涉及估算潜在结果 \(Y(t)\)，就好像它是一条线。
回归之所以有效是因为它可以正确估计 \(E[T|X]\) 用于正交化或者可以正确估计潜在结果 \(E[Y_t|X]\)，这赋予了它双重稳健性，你将在第 \(5\) 章中进一步探索。当你学习第 \(4\) 部分中的差异性时，通过这种视角来看待回归也将很重要。
实用例子
公立或私立学校？
在 \(Mastering \space Metrics\) 这本书中，\(Angrist\) 和 \(Pischke\) 展示了如何使用回归来调整在估计上私立学校对个人收入的影响时产生的偏差。私立学校的毕业生通常比公立学校的毕业生挣更多的钱，但是很难说这种关系有多少是因果关系。例如，你父母的收入可能会混淆这种关系，因为更富裕家庭的孩子更有可能上私立学校并赚更多的钱。类似地，由于私立学校非常挑剔，因此它们可能只招收那些无论如何都会过得更好的人。
因此，对收入与私立学校虚拟变量进行的朴素回归几乎肯定会返回一个积极的影响。换句话说，估计以下模型会给你一个积极的和显著的 \(\hat{\beta}_1\)：
然而，\(Angrist\) 和 \(Pischke\) 表明，如果你调整 \(SAT\) 分数和父母的收入，那么测量到的影响会减少。也就是说，如果你用这两个变量来扩充这个模型，你的 \(\hat{\beta}_1\) 会更小，相比于你从简短模型中得到的：\(income_i = \delta_0 + \beta_1 private + \delta_1 SAT_i + \delta_2 ParentInc_i + e_i\)
然而，\(Angrist\) 和 \(Pischke\) 表明，如果你调整 \(SAT\) 分数和父母的收入，那么测量到的影响会减少。也就是说，如果你用这两个变量来扩充这个模型，你的 \(\hat{\beta}_1\) 会更小，相比于你从简短模型中得到的：\(income_i = \delta_0 + \beta_1 private + \delta_1 SAT_i + \delta_2 ParentInc_i + e_i\)
即便在进行了包含父母收入的回归后，私立学校的影响仍然是积极的和显著的，至少在作者使用的数据集中是这样。然而，最后一组控制变量设法使这种关系变得不显著。作者纳入了学生申请的学校的平均\(SAT\)分数（无论他们是否被录取）。这可以被解释为一种对抱负的代理：
\(income_i = \delta_0 + \beta_1 private + \delta_1 SAT_i + \delta_2 ParentInc_i + \delta_3 AvgSATSchool_i + e_i\)
一旦他们加入了抱负代理控制变量，估计的 \(\hat{\beta}_1\) 就变得不显著了。有趣的是，仅保留这些控制变量并删除 \(SAT\) 和父母收入控制变量仍然导致一个不显著的估计。这表明，在给定你的抱负水平的情况下，你去公立学校还是私立学校并不重要，至少就你的收入而言是这样。
正性与外推
由于回归将潜在结果建模为参数函数，因此它允许在外推到您拥有所有处理水平的数据的区域之外。这可能是一种祝福或一种诅咒，这完全取决于外推是否合理。例如，考虑您必须在重叠度较低的数据集中估计处理效应，称其为 \(Dataset1\)。\(Dataset1\) 对于协变量 \(x\) 的高值没有控制单元，对于同一协变量的低值没有处理单元。如果你使用回归来估计该数据的处理效应，它会像第 \(1\) 个图中的直线所示那样，对 \(Y_0\) 和 \(Y_1\) 进行填补：
np.random.seed(1)
n = 1000
x = np.random.normal(0, 1, n)
t = np.random.normal(x, 0.2, n) > 0
y0 = x
y1 = 1 + x
y = np.random.normal((1-t)*y0 + t*y1, 0.2)
df_no_pos = pd.DataFrame(dict(x=x,t=t.astype(int),y=y))
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6), sharey=True)
sns.scatterplot(data=df_no_pos, x="x", y="y", hue="t", ax=ax1)
m0 = smf.ols("y~x", data=df_no_pos.query(f"t==0")).fit()
m1 = smf.ols("y~x", data=df_no_pos.query(f"t==1")).fit()
sns.lineplot(data=df_no_pos.assign(pred=m0.predict(df_no_pos)), x="x", y="pred", color=f"C0", ax=ax1)
sns.lineplot(data=df_no_pos.assign(pred=m1.predict(df_no_pos)), x="x", y="pred", color=f"C1", ax=ax1)
ax1.set_title("Dataset 1")
np.random.seed(1)
n = 1000
x = np.random.normal(0, 1, n)
t = np.random.binomial(1, 0.5, size=n)
y0 = x * (x<0) + (x>0)*0
y1 = 1 + x
y = np.random.normal((1-t)*y0 + t*y1, 0.2)
df_pos = pd.DataFrame(dict(x=x,t=t.astype(int),y=y))
sns.scatterplot(data=df_pos, x="x", hue="t", y="y", ax=ax2)
sns.lineplot(data=df_no_pos.assign(pred=m0.predict(df_pos)), x="x", y="pred", color=f"C0", ax=ax2)
sns.lineplot(data=df_no_pos.assign(pred=m1.predict(df_pos)), x="x", y="pred", color=f"C1", ax=ax2)
ax2.set_title("Dataset 2")
这没问题，只要你在低水平 \(x\) 的对照组中拟合的 \(Y_0\) 与 \(x\) 的关系在高 \(x\) 值下同样成立，并且你在处理组拟合的 \(Y_1\) 也能较好地外推到低 \(x\) 水平。通常，如果在协变量空间中有重叠的结果趋势看起来相似，少量的外推问题就不那么严重。
然而，过多的外推始终是危险的。假设你已经在 \(Dataset1\) 上估计了效果，但随后你收集了更多数据，现在随机分配处理。在这个 \(Dataset2\) 上，你会看到对于正的 \(x\) 值，效果越来越大。因此，如果你在这个新数据上评估之前的拟合，你会意识到你严重低估了处理的真实效果。这表明，在你没有正性的区域，你永远无法真正知道处理效果会发生什么。你可能会选择相信你对该区域的外推，但这并非没有风险。
线性回归中的非线性
到目前为止，处理效应曲线看起来相当线性。似乎信用额度的增加会导致风险持续增加，无论信用额度处于什么水平。从 \(1000\) 的额度增加到 \(2000\)，似乎与从 \(2000\) 的额度增加到 \(3000\)，对风险的增加幅度差不多。然而，你可能会遇到并非如此的情况。
举个例子，考虑和之前相同的数据，但现在你的任务是估计信用额度对信用卡支出的因果效应：
spend_data = pd.read_csv("./data/spend_data.csv")
spend_data.head()## wage educ exper ... credit_score2 credit_limit spend
## 0 950.0 11 16 ... 518.0 3200.0 3848
## 1 780.0 11 7 ... 429.0 1700.0 3144
## 2 1230.0 14 9 ... 571.0 4200.0 4486
## 3 1040.0 15 8 ... 411.0 1500.0 3327
## 4 1000.0 16 1 ... 518.0 1800.0 3508
##
## [5 rows x 8 columns]为简便起见，我们假设此处唯一的混杂变量是工资（假设银行在决定信用额度时仅参考这一信息）。该过程的因果图大致如下：
因此，你必须通过控制工资来识别信用额度对消费支出的影响。如果你想使用正交化方法估计这一效应，可以按以下思路操作：通过将信用额度对工资进行回归并提取残差，来消除信用额度中的偏差。到目前为止，这些都不是新内容。但这里存在一个问题：如果你按不同工资水平分别绘制 “消费支出 - 信用额度” 的关系图，会清晰地发现二者的关系并非线性的。
plt_df = (spend_data
.assign(wage_group = lambda d: pd.IntervalIndex(pd.qcut(d["wage"], 5)).mid)
.groupby(["wage_group", "credit_limit"])
[["spend"]]
.mean()
.reset_index())
plt.figure(figsize=(10,4))
sns.scatterplot(data=plt_df,
x="credit_limit",
y="spend",
hue="wage_group",
palette="gray")
相反，处理响应曲线似乎呈现出某种凹性特征：信用额度越高，该曲线的斜率就越低。或者用因果推断的术语来说 —— 由于斜率与因果效应密切相关 —— 你也可以表述为：随着信用额度的提高，信用额度对消费支出的影响会逐渐减弱。具体而言，信用额度从 \(1000\) 提升至 \(2000\) 时，消费支出的增长幅度会大于信用额度从 \(2000\) 提升至 \(3000\) 时的增长幅度。
处理变量线性化
为解决这一问题（处理响应曲线非线性），首先需要对处理变量进行变换，使其与结果变量呈现线性关系。例如，已知二者关系呈凹性，因此可尝试对信用额度处理变量施加某种凹函数变换。常用的候选变换函数包括对数函数、平方根函数，或任何将信用额度转化为分数次幂的函数。
本案例中，我们尝试使用平方根变换：
plt_df = (spend_data
# apply the sqrt function to the treatment
.assign(credit_limit_sqrt = np.sqrt(spend_data["credit_limit"]))
# create 5 wage binds for better vizualization
.assign(wage_group = pd.IntervalIndex(pd.qcut(spend_data["wage"], 5)).mid)
.groupby(["wage_group", "credit_limit_sqrt"])
[["spend"]]
.mean()
.reset_index())
plt.figure(figsize=(10,4))
sns.scatterplot(data=plt_df,
x="credit_limit_sqrt",
y="spend",
palette="gray",
hue="wage_group")
现在我们有进展了！信用额度的平方根似乎和消费支出存在线性关系。这当然不是完美的。要是你非常仔细地观察，还是能看到一些曲线形态。但就目前而言，它或许已经够用了。
我很遗憾地说，这个过程相当繁琐、依赖人工。你得尝试一大堆不同的方法，然后看看哪种方法能更好地让处理变量（这里指信用额度相关变换后的值）与结果变量呈现线性关系。一旦你找到了一个让自己满意的方法，在运行线性回归模型的时候就可以应用它。在这个例子中，这意味着你要估计的模型是：\(spend_i = \beta_0 + \beta_1 \sqrt{line_i} + e_i\)，而你的因果参数就是\(\beta_1\)。
这个模型可以用 \(statsmodels\) 来估计，方法是在公式中直接使用 \(Numpy\) 的平方根函数。
现在我们有进展了！信用额度的平方根似乎和消费支出存在线性关系。这当然不是完美的。要是你非常仔细地观察，还是能看到一些曲线形态。但就目前而言，它或许已经够用了。
model_spend = smf.ols('spend ~ np.sqrt(credit_limit)',data=spend_data).fit()
print(model_spend.summary().tables[1])## =========================================================================================
## coef std err t P>|t| [0.025 0.975]
## -----------------------------------------------------------------------------------------
## Intercept 493.0044 6.501 75.832 0.000 480.262 505.747
## np.sqrt(credit_limit) 63.2525 0.122 519.268 0.000 63.014 63.491
## =========================================================================================但事情还没结束。回想一下，工资会混淆信用额度和消费支出之间的关系。这一点你可以通过将上述模型的预测结果与原始数据对比作图来观察。你会发现模型的斜率很可能存在向上的偏差，原因是：工资越高，既会导致消费支出增加，也会导致信用额度提高。
plt_df = (spend_data
.assign(wage_group = lambda d: pd.IntervalIndex(pd.qcut(d["wage"], 5)).mid)
.groupby(["wage_group", "credit_limit"])
[["spend"]]
.mean()
.reset_index()
)
x = np.linspace(plt_df["credit_limit"].min(), plt_df["credit_limit"].max())
plt.figure(figsize=(10,4))
plt.plot(x, model_spend.params[0] + model_spend.params[1]*np.sqrt(x), color="C1", label="prediction", lw=4)
plt.legend()
sns.scatterplot(data=plt_df,
x="credit_limit",
y="spend",
palette="gray",
hue="wage_group")
如果将工资纳入模型：\(spend_i = \beta_0 + \beta_1 \sqrt{line_i} + \beta_2 wage_i + e_i\)
然后再次估计 \(\beta_1\)，你会得到信用额度对消费支出影响的无偏估计（当然，假设工资是唯一的混杂因素）。这个估计值比你之前得到的要小，这是因为把工资纳入模型修正了向上的偏差。
非线性 FWL 与去偏
至于 \(\text{FWl}\) 定理在非线性数据上如何发挥作用，其实和之前的情况完全类似，只是现在你得先应用非线性变换。也就是说，你可以将用线性回归估计非线性模型的过程分解如下：
如果将工资纳入模型：\(spend_i = \beta_0 + \beta_1 \sqrt{line_i} + \beta_2 wage_i + e_i\)
然后再次估计 \(\beta_1\)，你会得到信用额度对消费支出影响的无偏估计（当然，假设工资是唯一的混杂因素）。这个估计值比你之前得到的要小，这是因为把工资纳入模型修正了向上的偏差。
找到一个函数 \(F\)，它能使处理变量 \(T\) 和结果变量 \(Y\) 之间的关系线性化。
一个去偏步骤：将经变换的处理变量 \(F(T)\) 对混杂变量 \(X\) 进行回归，得到处理残差 \(\widetilde{F(T)} = F(T) - \widehat{F(T)}\)。
一个去噪步骤：将结果变量 \(Y\) 对混杂变量 \(X\) 进行回归，得到结果残差 \(\widetilde{Y} = Y - \widehat{Y}\)。
一个结果模型：将结果残差 \(\widetilde{Y}\) 对处理残差 \(\widetilde{F(T)}\) 进行回归，以得到 \(F(T)\) 对 \(Y\) 因果效应的估计值。
在这个例子中，\(F\) 是平方根函数，所以下面就是考虑非线性时如何应用 \(\text{FWL}\) 定理的方法（我还分别将 \(\overline{F(lines)}\) 和 \(\overline{spend}\) 加到 \(\widetilde{F(T)}\) 和 \(\widetilde{Y}\) 中。这一步是可选的，但它能让可视化效果更好）：
debias_spend_model = smf.ols(f'np.sqrt(credit_limit) ~ wage',data=spend_data).fit()
denoise_spend_model = smf.ols(f'spend ~ wage', data=spend_data).fit()
credit_limit_sqrt_deb = debias_spend_model.resid + np.sqrt(spend_data["credit_limit"]).mean()
spend_den = denoise_spend_model.resid + spend_data["spend"].mean()
spend_data_deb = (spend_data.assign(credit_limit_sqrt_deb = credit_limit_sqrt_deb, spend_den = spend_den))
final_model = smf.ols(f'spend_den ~ credit_limit_sqrt_deb',data=spend_data_deb).fit()
print(final_model.summary().tables[1])## =========================================================================================
## coef std err t P>|t| [0.025 0.975]
## -----------------------------------------------------------------------------------------
## Intercept 1493.6990 3.435 434.818 0.000 1486.966 1500.432
## credit_limit_sqrt_deb 43.8504 0.065 672.640 0.000 43.723 43.978
## =========================================================================================<<<<<<< HEAD 毫不奇怪，通过运行包含 \(wage\) 和 \(credit\_limit\) 的完整模型，你在这里得到的 \(\beta_1\) 估计值与之前得到的完全相同。而且，如果你将这个模型的预测结果与原始数据作图对比，会发现它不再像之前那样存在向上的偏差。相反，它正好穿过不同工资组的中间位置。
plt_df = (spend_data
.assign(wage_group = lambda d: pd.IntervalIndex(pd.qcut(d["wage"], 5)).mid)
.groupby(["wage_group", "credit_limit"])
[["spend"]]
.mean()
.reset_index())
x = np.linspace(plt_df["credit_limit"].min(), plt_df["credit_limit"].max())
plt.figure(figsize=(10,4))
plt.plot(x, (final_model.params[0] + final_model.params[1]*np.sqrt(x)), color="C1", label="prediction", lw=4)
plt.legend()
sns.scatterplot(data=plt_df, x="credit_limit",y="spend",palette="gray",hue="wage_group")
虚拟变量回归
回归和正交化固然都很不错，但归根结底你得做出一个独立性假设。你必须假定：在考虑了一些协变量之后，处理变量就如同是随机分配的一样。这可能很牵强。要知道模型中是否包含了所有的混杂因素是非常困难的。正因如此，只要有可能大力推行随机实验是很有意义的。比如在银行业的例子中，如果信用额度是随机分配的，那会很棒。因为这样就能相当直接地估计它对违约率和客户消费的影响。问题在于：这样的实验会极其昂贵。你会把随机的信用额度给到风险很高的客户，进而这些客户很可能会违约从而造成巨额损失。
条件随机实验
解决这一难题的方法并非理想的随机对照试验，但却是次优之选：分层随机实验或条件随机实验。你无需设计一个信用额度完全随机且从同一概率分布中抽取的实验，而是创建多个局部实验，根据客户的协变量从不同分布中抽取样本。例如，你知道 \(credit\_score1\) 变量是客户风险的一个代理指标。因此，你可以用它来划分风险程度不同的客户群体，将他们分成 \(credit\_score1\) 相近的组。然后，对于高风险组，你从平均水平较低的分布中抽样来随机分配信用额度；对于低风险客户，你从平均水平较高的分布中抽样来随机分配信用额度：
risk_data_rnd = pd.read_csv("./data/risk_data_rnd.csv")
risk_data_rnd.head()## wage educ exper ... credit_score1_buckets credit_limit default
## 0 890.0 11 16 ... 400 5400.0 0
## 1 670.0 11 7 ... 200 3800.0 0
## 2 1220.0 14 9 ... 400 5800.0 0
## 3 1210.0 15 8 ... 600 6500.0 0
## 4 900.0 16 1 ... 200 2100.0 0
##
## [5 rows x 9 columns]绘制按 \(credit\_score1\_buckets\) 划分的信用额度直方图，你可以看到信用额度是从不同分布中抽样得到的。得分较高的低风险客户其直方图的分布向左偏斜信用额度更高。得分较低风险更高的客户组的信用额度来自向右偏斜的分布其信用额度更低。这类实验所研究的信用额度与可能的最优额度相差不大，这就将测试成本降低到了更易于管理的水平：
plt.figure(figsize=(15,3))
sns.histplot(data=risk_data_rnd, x="credit_limit", hue="credit_score1_buckets", kde=True, palette="gray")
plt.title("Conditional random experiment")
贝塔抽样：在这个实验中，信用额度是从贝塔分布中抽样得到的。贝塔分布可以被看作是均匀分布的推广，这使得当你希望样本被限制在特定范围内时，它特别好用。
这并不意味着条件随机实验比完全随机实验更好。它们确实更省钱，但也增加了大量额外的复杂性。因此，如果你出于某种原因选择进行条件随机实验，要尽量让它尽可能接近完全随机实验。这意味着：
分组数量越少，处理条件随机试验就越容易。在这个例子中，你只有 \(5\) 个组，因为你把 \(credit\_score1\) 按每 \(200\) 分为一组，而分数范围是从 \(0\) 到 \(1000\)。将具有不同处理分布的不同组结合起来会增加复杂性，所以坚持少分组是个好主意。
各组之间处理分布的重叠程度越大，你操作起来就越轻松。这与正性假设有关。在这个例子中，如果高风险组获得高信用额度的概率为零，你就不得不依赖危险的外推来了解如果他们获得那些高信用额度会发生什么。
如果你把这两个经验法则推到极致，就又回到了完全随机实验的情况，这意味着它们两者都存在权衡：分组数量越少、重叠程度越高，实验就越容易解读，但也会更昂贵。反之亦然。
注意：分层实验也可作为一种工具，用于最小化方差，并确保处理组和对照组在分层变量上保持平衡。但在那些应用场景中，处理分布被设计为在所有组或层中是相同的。
虚拟变量
条件随机实验的巧妙之处在于，条件独立性假设要合理得多，因为你知道，在给定你所选的分类变量的情况下，信用额度是随机分配的。缺点是仅对结果变量和处理变量进行简单回归会得到有偏差的估计。例如，当你在模型中不纳入混杂因素时，估计模型会出现如下情况：\(default_i = \beta_0 + \beta_1 lines_i + e_i\)
model = smf.ols("default ~ credit_limit", data=risk_data_rnd).fit()
print(model.summary().tables[1])## ================================================================================
## coef std err t P>|t| [0.025 0.975]
## --------------------------------------------------------------------------------
## Intercept 0.1369 0.009 15.081 0.000 0.119 0.155
## credit_limit -9.344e-06 1.85e-06 -5.048 0.000 -1.3e-05 -5.72e-06
## ================================================================================如你所见，因果参数估计值 \(\widehat{\beta_1}\) 为负，这在此处毫无意义。更高的信用额度很可能不会降低客户的风险。出现这种情况是因为，在这些数据中，由于实验的设计方式，低风险客户 —— 那些 \(credit\_score1\) 高的客户 —— 平均而言获得了更高的信用额度。
为了对此进行调整，你需要在模型中纳入处理变量被随机分配所处的组。在这种情况下，你需要对 \(credit\_score1\_buckets\) 进行控制。尽管这个组是以数字形式呈现的，但它实际上是一个分类变量：它代表一个组。因此，控制组本身的方法是创建虚拟变量。虚拟变量是针对一个组的二元列。如果客户属于该组，其值为 \(1\)，否则为 \(0\)。由于一个客户只能来自一个组，所以最多只有一个虚拟列的值为 \(1\)，其余都为 \(0\)。如果你的背景是机器学习，你可能知道这就是独热编码，它们完全是一回事。
在 \(pandas\) 中，你可以使用 \(pd.get\_dummies\) 函数来创建虚拟变量。在这里，我传入代表组的列 \(credit\_score1\_buckets\)，并指定我想要创建带有后缀 \(\_sb\)（代表 \(score\_bucket\)）的虚拟列。此外，我删除了第一个虚拟变量，即 \(0\) 到 \(200\) 分组的那个。这是因为其中一个虚拟列是冗余的。如果我知道所有其他列都为 \(0\)，那么我删除的那个列的值肯定为 \(1\)。
risk_data_dummies = risk_data_rnd.join(pd.get_dummies(risk_data_rnd["credit_score1_buckets"],prefix="sb",drop_first=True))
risk_data_dummies.head()## wage educ exper married ... sb_400 sb_600 sb_800 sb_1000
## 0 890.0 11 16 1 ... True False False False
## 1 670.0 11 7 1 ... False False False False
## 2 1220.0 14 9 1 ... True False False False
## 3 1210.0 15 8 1 ... False True False False
## 4 900.0 16 1 1 ... False False False False
##
## [5 rows x 14 columns]一旦你有了虚拟变量列，就可以将它们添加到你的模型中，然后再次估计 \(\beta_1\)：
\(default_i = \beta_0 + \beta_1 lines_i + \theta G_i + e_i\)
现在，你会得到一个更合理的估计值，它至少是正的，这表明更高的信用额度会增加违约风险。
model = smf.ols("default ~ credit_limit + sb_200 + sb_400 + sb_600 + sb_800 + sb_1000",data=risk_data_dummies).fit()
print(model.summary().tables[1])## ===================================================================================
## coef std err t P>|t| [0.025 0.975]
## -----------------------------------------------------------------------------------
## Intercept 0.2253 0.056 4.000 0.000 0.115 0.336
## sb_200[T.True] -0.0559 0.057 -0.981 0.327 -0.168 0.056
## sb_400[T.True] -0.1442 0.057 -2.538 0.011 -0.256 -0.033
## sb_600[T.True] -0.2148 0.057 -3.756 0.000 -0.327 -0.103
## sb_800[T.True] -0.2489 0.060 -4.181 0.000 -0.366 -0.132
## sb_1000[T.True] -0.2541 0.094 -2.715 0.007 -0.438 -0.071
## credit_limit 4.652e-06 2.02e-06 2.305 0.021 6.97e-07 8.61e-06
## ===================================================================================我只是向你展示如何手动创建虚拟变量，这样你就能了解底层发生了什么。如果你需要在非 \(Python\) 的其他框架中实现这类回归，这会非常有用。在 \(Python\) 中，如果你使用\(statsmodels\)，公式中的\(C()\)函数可以为你完成所有这些工作。
model = smf.ols("default ~ credit_limit + C(credit_score1_buckets)",data=risk_data_rnd).fit()
print(model.summary().tables[1])## ====================================================================================================
## coef std err t P>|t| [0.025 0.975]
## ----------------------------------------------------------------------------------------------------
## Intercept 0.2253 0.056 4.000 0.000 0.115 0.336
## C(credit_score1_buckets)[T.200] -0.0559 0.057 -0.981 0.327 -0.168 0.056
## C(credit_score1_buckets)[T.400] -0.1442 0.057 -2.538 0.011 -0.256 -0.033
## C(credit_score1_buckets)[T.600] -0.2148 0.057 -3.756 0.000 -0.327 -0.103
## C(credit_score1_buckets)[T.800] -0.2489 0.060 -4.181 0.000 -0.366 -0.132
## C(credit_score1_buckets)[T.1000] -0.2541 0.094 -2.715 0.007 -0.438 -0.071
## credit_limit 4.652e-06 2.02e-06 2.305 0.021 6.97e-07 8.61e-06
## ====================================================================================================最后，这里你只有一个斜率参数。添加虚拟变量来控制混杂因素会为每个组提供一个截距，但所有组的斜率是相同的。我们很快会讨论这个问题，不过这个斜率会是每个组内回归的方差加权平均值。如果你绘制该模型对每个组的预测结果，能清楚看到每个组对应一条直线，但所有直线的斜率都相同。
作为方差加权平均的回归
但是，如果饱和回归和按组计算效应得到的结果完全一致，你可能会问自己一个非常重要的问题。当你运行 \(default \sim credit\_limit + C(credit\_score1\_buckets)\)这个模型（不含交互项）时，会得到一个单一的效应：只有一个斜率参数。重要的是，如果你回头看，这个效应估计值和你通过按组估计效应、并以组规模为权重对结果取平均所得到的结果是不同的。所以，在某种程度上，回归是在整合不同组的效应。而且，它整合的方式并非是样本量加权平均。那它到底是怎样的呢？
要回答这个问题，最好的方法还是使用一些非常具有说明性的模拟数据。在这里，让我们模拟来自两个不同组的数据。组 \(1\) 有 \(1000\) 个样本，平均处理效应为 \(1\)。组 \(2\) 有 \(500\) 个样本，平均处理效应为 \(2\)。此外，组 \(1\) 中处理变量的标准差是 \(1\)，组 \(2\) 中是 \(2\)：
np.random.seed(123)
# std(t)=1
t1 = np.random.normal(0, 1, size=1000)
df1 = pd.DataFrame(dict(t=t1, y=1*t1,g=1))
# std(t)=2
t2 = np.random.normal(0, 2, size=500)
df2 = pd.DataFrame(dict(t=t2, y=2*t2,g=2))
df = pd.concat([df1, df2])
df.head()## t y g
## 0 -1.085631 -1.085631 1
## 1 0.997345 0.997345 1
## 2 0.282978 0.282978 1
## 3 -1.506295 -1.506295 1
## 4 -0.578600 -0.578600 1如果你分别估计每个组的效应，然后以组规模为权重对结果取平均，你会得到大约 \(1.33\) 的 \(ATE\)，即\((1×1000 + 2×500) / 1500\)。
def regress(group, y, t):
# 函数实现，使用group作为当前分组的数据
import statsmodels.api as sm
X = sm.add_constant(group[t])
model = sm.OLS(group[y], X).fit()
return model.params[t] # 返回处理变量t的系数
effect_by_group = df.groupby("g").apply(regress, y="y", t="t")
ate = (effect_by_group * df.groupby("g").size()).sum() / df.groupby("g").size().sum()
ate## 1.333333333333333但如果你在控制组的情况下，对 \(y\) 关于 \(t\) 进行回归，会得到一个非常不同的结果。现在，综合效应更接近组 \(2\) 的效应，尽管组 \(2\) 的样本量只有组 \(1\) 的一半。
model = smf.ols("y ~ t + C(g)", data=df).fit()
model.params## Intercept 0.024758
## C(g)[T.2] 0.019860
## t 1.625775
## dtype: float64原因在于回归并非以样本量作为权重来整合组效应。相反，它使用与每个组内处理变量方差成比例的权重。回归更倾向于处理变量变化大的组。乍一听这似乎有些奇怪，但仔细想想是很有道理的。如果在一个组内处理变量变化不大，你怎么能确定它的效应呢？如果处理变量变化很大，它对结果的影响就会更明显。
总而言之，如果你有多个组且每个组内的处理变量是随机分配的，条件性原则表明：总体效应是每个组内效应的加权平均：\(ATE = E\left\{ \frac{\partial}{\partial t} E[Y_i \vert T = t, Group_i] \ w(Group_i) \right\}\)
根据方法的不同权重会有所不同。对于回归，\(w(Group_i) \propto \sigma^2(T \vert Group)\)，但你也可以选择手动以样本量为权重来对组效应进行加权：\(w(Group_i) = N_{Group}\)。
另见
了解这种差异是理解回归背后运作机制的关键。回归以方差对组效应进行加权，这一点即使是最优秀的研究者也需要不断被提醒。\(2020\)
年，计量经济学领域在双重差分法（你会在第 \(4\)
部分看到更多相关内容）方面经历了一次复兴。问题的核心在于回归并非以样本量对效应进行加权。如果你想了解更多相关内容，我建议你查阅
\(Andrew \space Goodman-Bacon\)
的论文《Treatment Timing Variation in Difference - in -
Differences》，或者只需等到我们讲到第 \(4\) 部分。
去均值与固定效应
你刚刚了解了如何在模型中纳入虚拟变量，以解释不同组间处理分配的差异。而正是在处理虚拟变量时，\(FWL\) 定理真正发挥了作用。要是你有大量的组，为每个组都添加一个虚拟变量，不仅繁琐，而且计算成本很高。你会创建大量大多为 \(0\) 的列。通过应用 \(FWL\) 定理，并理解在涉及虚拟变量时回归是如何对处理变量进行正交化的，你就能轻松解决这个问题。
你已经知道，\(FWL\) 定理中的去偏步骤包括从协变量（在这种情况下就是虚拟变量）预测处理变量：
model_deb = smf.ols("credit_limit ~ C(credit_score1_buckets)",data=risk_data_rnd).fit()
print(model_deb.summary().tables[1])## ====================================================================================================
## coef std err t P>|t| [0.025 0.975]
## ----------------------------------------------------------------------------------------------------
## Intercept 1173.0769 278.994 4.205 0.000 626.193 1719.961
## C(credit_score1_buckets)[T.200] 2195.4337 281.554 7.798 0.000 1643.530 2747.337
## C(credit_score1_buckets)[T.400] 3402.3796 279.642 12.167 0.000 2854.224 3950.535
## C(credit_score1_buckets)[T.600] 4191.3235 280.345 14.951 0.000 3641.790 4740.857
## C(credit_score1_buckets)[T.800] 4639.5105 291.400 15.921 0.000 4068.309 5210.712
## C(credit_score1_buckets)[T.1000] 5006.9231 461.255 10.855 0.000 4102.771 5911.076
## ====================================================================================================由于虚拟变量本质上是组的平均值，你可以看到，在这个模型中正是在做这样的预测：如果 \(credit\_score1\_buckets=0\)，你就在预测 \(credit\_score1\_buckets=0\) 组的平均信用额度；如果 \(credit\_score\_buckets=1\)，你就在预测 \(credit\_score1\_buckets=1\) 组的平均信用额度（通过将截距与该组的系数相加得到的，即\(1173.0769 + 2195.4337 = 3368.510638\)），依此类推。这些正是组的平均值。
risk_data_rnd.groupby("credit_score1_buckets")["credit_limit"].mean()## credit_score1_buckets
## 0 1173.076923
## 200 3368.510638
## 400 4575.456498
## 600 5364.400448
## 800 5812.587413
## 1000 6180.000000
## Name: credit_limit, dtype: float64这意味着，如果你想对处理变量进行残差化，你可以用一种更简单且有效的方式来做。首先，计算每个组的平均处理值：
risk_data_fe = risk_data_rnd.assign(credit_limit_avg = lambda d: (d.groupby("credit_score1_buckets")["credit_limit"].transform("mean")))然后，要得到残差，就从处理变量中减去该组的平均值。由于这种方法减去了处理变量的平均值，它有时被称为对处理变量进行去均值。如果你想在回归公式中进行这样的操作，必须用 \(I(...)\) 把数学运算括起来：
model = smf.ols("default ~ I(credit_limit-credit_limit_avg)",data=risk_data_fe).fit()
print(model.summary().tables[1])## ======================================================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------------------------------
## Intercept 0.0935 0.003 32.121 0.000 0.088 0.099
## I(credit_limit - credit_limit_avg) 4.652e-06 2.05e-06 2.273 0.023 6.4e-07 8.66e-06
## ======================================================================================================你在这里得到的参数估计值和你在模型中添加虚拟变量时得到的完全相同。这是因为，从数学角度而言它们是等价的。这种方法被称为固定效应，因为你在控制组内任何固定不变的因素。它源自具有时间结构的面板数据的因果推断文献，你会在第四部分进一步探讨。
来自同一文献的另一个思路是，在回归模型中纳入按组计算的平均处理（值）（源于蒙德拉克，1978 年）。回归会对纳入的额外变量所对应的处理（变量）进行残差化，所以这里的效应大致相同：
model = smf.ols("default ~ credit_limit + credit_limit_avg",data=risk_data_fe).fit()
print(model.summary().tables[1])## ====================================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------------
## Intercept 0.4325 0.020 21.418 0.000 0.393 0.472
## credit_limit 4.652e-06 2.02e-06 2.305 0.021 6.96e-07 8.61e-06
## credit_limit_avg -7.763e-05 4.75e-06 -16.334 0.000 -8.69e-05 -6.83e-05
## ====================================================================================实际案例
营销组合建模
衡量广告对销售额的影响非常困难，因为通常你无法随机决定谁能看到你的广告。在广告行业，一种流行的替代随机化的方法是称为营销组合建模的技术。尽管名字听起来很花哨，但营销组合模型其实就是将销售额对营销策略指标和一些混杂因素进行回归。例如，假设你想了解电视、社交媒体和搜索广告的预算对你产品销售额的影响。你可以运行一个回归模型，其中每个单元 \(i\) 是一天：
\(\begin{align*} Sales_i &= \beta_0 + \beta_1 TV_i + \beta_2 Social_i + \beta_3 Search_i + \delta_1 CompetitorSales_i + \delta_2 Month_i + \delta_3 Trend_i + e_i \end{align*}\)
为了考虑到你可能在销售情况好的月份增加了营销预算这一情况，你可以通过在回归中纳入额外的控制变量来对这个混杂因素进行调整。例如，你可以纳入竞争对手的销售额、每个月的虚拟变量以及一个趋势变量。
遗漏变量偏误：从回归视角看混杂问题
我希望在第 \(3\)
章中已经把话说得非常清楚了：共同原因即混杂因素会使处理变量与结果变量之间的估计关系产生偏误。这就是为什么你需要对它们加以考虑，比如将其纳入回归模型中。然而，回归对于混杂偏误有着其独特的解读视角。当然，到目前为止所说的所有内容仍然成立，但回归能让你更精确地理解混杂偏误。例如，假设你想估计信用额度对违约行为的影响，且工资是唯一的混杂因素：
在这种情况下，你知道应该估计包含混杂因素的模型：\(default_i = \beta_0 + \beta_1 lines_i + \beta_2 wage_i + e_i\)
但如果你转而估计一个省略了混杂因素的更简单的模型：\(default_i = \beta_0 + \beta_1 lines_i + e_i\)
得到的估计值就会产生偏差：
short_model = smf.ols("default ~ credit_limit", data=risk_data).fit()
short_model.params["credit_limit"]## -2.401961992596881e-05如你所见，看起来更高的信用额度会导致违约率下降，这是毫无意义的。但你已经知道这一点了。你不知道的是，你可以精确计算出这种偏差的大小。通过回归，你可以说，由于遗漏变量导致的偏差等于该变量在包含它的模型中的效应，加上遗漏变量对结果的效应乘以遗漏变量对纳入变量的回归系数。别担心，我知道这听起来很拗口，所以让我们一点一点来理解。首先，这意味着对 \(Y\) 关于 \(T\) 的简单回归会是真实的因果参数 \(\tau\)，加上一个偏差项：
\(\frac{Cov(T,Y)}{Var(T)} = \tau + \beta'_{omitted}\delta_{omitted}\)
这个偏差项是遗漏的混杂因素对结果的系数 \(\beta_{omitted}\)，乘以将遗漏变量对处理变量进行回归得到的系数 \(\delta_{omitted}\)。为了验证这一点，你可以用下面的代码得到你之前得到的有偏差的参数估计值，这段代码重现了遗漏变量偏差公式：
long_model = smf.ols("default ~ credit_limit + wage",data=risk_data).fit()
omitted_model = smf.ols("wage ~ credit_limit", data=risk_data).fit()
(long_model.params["credit_limit"] + long_model.params["wage"]*omitted_model.params["credit_limit"])## -2.401961992596883e-05中性控制变量
到目前为止，你可能已经对回归如何调整混杂变量有了很好的认识。如果你想在调整混杂变量 \(X\) 的同时，了解处理变量 \(T\) 对 \(Y\) 的效应，你所要做的就是将 \(X\) 纳入模型。或者，为了得到完全相同的结果，你可以从 \(X\) 预测 \(T\) 得到残差，并将其用作处理变量的去偏版本。对 \(Y\) 关于这些残差进行回归，会给出在保持 \(X\) 固定的情况下，\(T\) 与 \(Y\) 的关系。
但你应该将哪些类型的变量纳入 \(X\) 呢？同样，并不是因为添加变量能对其进行调整，你就想把所有东西都纳入回归模型。正如前面章节所讲的，你不想纳入对撞因子或中介变量，因为它们会引发选择偏误。但在回归的背景下，还有更多类型的控制变量是你应该了解的。有些控制变量起初看似无害，但实际上相当有害。这些控制变量被称为中性的，因为它们不会影响你回归估计中的偏差。但它们在方差方面会产生严重的影响。正如你将看到的在回归中纳入某些变量时，存在偏差—方差权衡。例如，考虑以下有向无环图：
你应该把 \(credit\_score2\) 纳入你的模型吗？如果不纳入它，你会得到和之前一直看到的相同结果。这个结果是无偏的，因为你在对 \(credit\_score\_buckets\) 进行调整。但尽管你没必要纳入它，还是来看看纳入 \(credit\_score2\) 时会发生什么。将以下结果与你之前没纳入 \(credit\_score2\) 时得到的结果做比较。有什么变化？
formula = "default~credit_limit+C(credit_score1_buckets)+credit_score2"
model = smf.ols(formula, data=risk_data_rnd).fit()
print(model.summary().tables[1])## ====================================================================================================
## coef std err t P>|t| [0.025 0.975]
## ----------------------------------------------------------------------------------------------------
## Intercept 0.5576 0.055 10.132 0.000 0.450 0.665
## C(credit_score1_buckets)[T.200] -0.0387 0.055 -0.710 0.478 -0.146 0.068
## C(credit_score1_buckets)[T.400] -0.1032 0.054 -1.898 0.058 -0.210 0.003
## C(credit_score1_buckets)[T.600] -0.1410 0.055 -2.574 0.010 -0.248 -0.034
## C(credit_score1_buckets)[T.800] -0.1161 0.057 -2.031 0.042 -0.228 -0.004
## C(credit_score1_buckets)[T.1000] -0.0430 0.090 -0.479 0.632 -0.219 0.133
## credit_limit 4.928e-06 1.93e-06 2.551 0.011 1.14e-06 8.71e-06
## credit_score2 -0.0007 2.34e-05 -30.225 0.000 -0.001 -0.001
## ====================================================================================================首先，\(credit\_limit\) 的参数估计值变高了一点。但更重要的是标准误差减小了。这是因为 \(credit\_score2\) 是结果 \(Y\) 的良好预测因子，它会对线性回归的去噪步骤产生作用。在 \(FWL\) 定理的最后一步，由于纳入了 \(credit\_score2\)，\(\tilde{Y}\)的方差会被降低，对其关于 \(\tilde{T}\) 进行回归会得到更精确的结果。
这是线性回归一个非常有趣的性质。它表明线性回归不仅可以用于调整混杂因素，还可以用于降低噪声。例如，如果你有来自设计良好的 \(A/B\) 测试的数据，你不必担心偏差问题。但你仍然可以将回归用作降噪工具，只需纳入对结果有高度预测性且不会引发选择偏差的变量即可。
降噪技术
还有其他的降噪技术。最著名的是 \(CUDPED\)（控制单元后验差异估计），它由微软的研究人员开发，在科技公司中被广泛使用。\(CUDPED\) 与仅执行 \(FWL\) 定理的去噪部分非常相似。
引发噪声的控制变量
就像控制变量能降低噪声一样，它们也能增加噪声。例如，再考虑一个条件随机实验的情况。但这次，你感兴趣的是信用额度对支出的影响而非对风险的影响。和之前的例子一样，在给定 \(credit\_score1\) 的情况下，信用额度是随机分配的。但这次，假设 \(credit\_score1\) 不是一个混杂因素。它会引发处理变量，但不会引发结果变量。这个数据生成过程的因果图如下：
这意味着，要得到信用额度对支出的因果效应，你不需要对 \(credit\_score1\) 进行调整。一个单变量回归模型就足够了。在这里，我保留了平方根函数，以解释处理响应函数中的凹性：
spend_data_rnd = pd.read_csv("data/spend_data_rnd.csv")
model = smf.ols("spend ~ np.sqrt(credit_limit)",data=spend_data_rnd).fit()
print(model.summary().tables[1])## =========================================================================================
## coef std err t P>|t| [0.025 0.975]
## -----------------------------------------------------------------------------------------
## Intercept 2153.2154 218.600 9.850 0.000 1723.723 2582.708
## np.sqrt(credit_limit) 16.2915 2.988 5.452 0.000 10.420 22.163
## =========================================================================================但是，如果你确实纳入了\(credit\_score1\_buckets\)，会发生什么呢？
model = smf.ols("spend~np.sqrt(credit_limit)+C(credit_score1_buckets)",data=spend_data_rnd).fit()
print(model.summary().tables[1])## ====================================================================================================
## coef std err t P>|t| [0.025 0.975]
## ----------------------------------------------------------------------------------------------------
## Intercept 2367.4867 556.273 4.256 0.000 1274.528 3460.446
## C(credit_score1_buckets)[T.200] -144.7921 591.613 -0.245 0.807 -1307.185 1017.601
## C(credit_score1_buckets)[T.400] -118.3923 565.364 -0.209 0.834 -1229.211 992.427
## C(credit_score1_buckets)[T.600] -111.5738 570.471 -0.196 0.845 -1232.429 1009.281
## C(credit_score1_buckets)[T.800] -89.7366 574.645 -0.156 0.876 -1218.791 1039.318
## C(credit_score1_buckets)[T.1000] 363.8990 608.014 0.599 0.550 -830.720 1558.518
## np.sqrt(credit_limit) 14.5953 3.523 4.142 0.000 7.673 21.518
## ====================================================================================================你可以看到，它纳入该变量增大了标准误差，拓宽了因果参数的置信区间。这是因为，正如你在《回归作为方差加权平均》中看到的，\(OLS\)偏好处理变量具有高方差的情况。但如果你控制了一个能解释处理变量的协变量，实际上就是在降低其方差。
特征选择：偏差—方差权衡
实际上，很难出现一个协变量仅引发处理变量而不引发结果变量的情况。大多数情况下，你会有一堆混杂因素，它们对 \(T\) 和 \(Y\) 都有引发作用，只是程度不同。在 \(图4-3\) 中，\(X_1\)是\(T\) 的强引发因素，但却是 \(Y\) 的弱引发因素；\(X_3\)是 \(Y\) 的强引发因素，但却是 \(T\) 的弱引发因素；而\(X_2\)则处于中间位置，这由每个箭头的粗细来表示。
在这些情况下，你很快就会陷入两难的境地。一方面，如果你想消除所有偏差，就必须纳入所有协变量；毕竟，它们是需要进行调整的混杂因素。另一方面，对处理变量的引发因素进行调整会增加估计量的方差。
为了看清这一点，让我们根据 \(图4-3\) 中的因果图模拟数据。这里，真实的 \(ATE\) 是 \(0.5\)。如果你在控制所有混杂因素的情况下尝试估计这个效应，估计值的标准误差会过高，以至于无法得出任何结论：
np.random.seed(123)
n = 100
(x1, x2, x3) = (np.random.normal(0, 1, n) for _ in range(3))
t = np.random.normal(10*x1 + 5*x2 + x3)
# ate = 0.05
y = np.random.normal(0.05*t + x1 + 5*x2 + 10*x3, 5)
df = pd.DataFrame(dict(y=y, t=t, x1=x1, x2=x2, x3=x3))
print(smf.ols("y ~ t + x1 + x2 + x3", data=df).fit().summary().tables[1])## ==============================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------
## Intercept 0.2707 0.527 0.514 0.608 -0.775 1.316
## t 0.8664 0.607 1.427 0.157 -0.339 2.072
## x1 -7.0628 6.038 -1.170 0.245 -19.049 4.923
## x2 0.0143 3.128 0.005 0.996 -6.195 6.224
## x3 9.6292 0.887 10.861 0.000 7.869 11.389
## ==============================================================================如果你知道某个混杂因素是处理变量的强预测因子，却是结果变量的弱预测因子，你可以选择将其从模型中剔除。在这个例子中，那会是 \(X_1\)。现在，注意了！这会使你的估计产生偏差。但如果它也能显著降低方差，或许这是值得付出的代价：
print(smf.ols("y ~ t + x2 + x3", data=df).fit().summary().tables[1])## ==============================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------
## Intercept 0.1889 0.523 0.361 0.719 -0.849 1.227
## t 0.1585 0.046 3.410 0.001 0.066 0.251
## x2 3.6095 0.582 6.197 0.000 2.453 4.766
## x3 10.4549 0.537 19.453 0.000 9.388 11.522
## ==============================================================================归根结底，你在模型中纳入调整的混杂因素越多，因果估计中的偏差就越小。然而，如果你纳入的变量是结果变量的弱预测因子，却是处理变量的强预测因子，这种偏差的减少会以方差大幅增加为沉重代价。换种说法，有时候为了降低方差，接受一点偏差是值得的。此外，你应该清楚并非所有混杂因素都是等同的。当然，它们都是\(T\) 和 \(Y\) 的共同原因。但如果它们对处理变量的解释程度过高，而对结果变量的解释几乎为零，你真的应该考虑在调整时将其剔除。这对回归是有效的，对其他调整策略（如倾向得分加权见第\(5\)章）也同样适用。
不幸的是，就解释处理变量而言，混杂因素要弱到什么程度才值得将其剔除在因果推断中仍是一个未解决的问题。不过，了解这种偏差—方差权衡的存在是很有价值的，因为它有助于你理解和解释线性回归中出现的情况。
核心观点
本章内容围绕回归展开，但视角与你在机器学习书籍中常见的大不相同。这里的回归并非预测工具，注意到我甚至一次都没提及 \(R^2\)！相反，回归在这里主要被用作调整混杂因素的方法，有时也作为一种方差降低技术。
本章的核心是正交化，它是在条件独立性成立时让处理变量尽可能接近随机分配的一种手段。形式上如果\(Y_t \perp T|X\)，你可以通过对\(T\)关于\(X\)进行回归并获取残差来调整由\(X\)导致的混杂偏差。这些残差可被视为处理变量的去偏版本。
这种方法借助 \(FWL\) 定理得到了进一步发展，该定理指出多元回归可分解为以下步骤：
去偏步骤：将处理变量 \(T\) 对混杂因素 \(X\) 进行回归，得到处理残差 \(\tilde{T} = T - \hat{T}\)。
去噪步骤：将结果变量 \(Y\) 对混杂变量 \(X\) 进行回归，得到结果残差 \(\tilde{Y} = Y - \hat{Y}\)。
结果模型：将结果残差 \(\tilde{Y}\) 对处理残差 \(\tilde{T}\) 进行回归，以获取 \(T\) 对 \(Y\) 的因果效应估计值。
本章的其他所有内容都源于这个定理 —— 无论是非线性处理响应函数、理解含分类变量的回归如何实现加权平均，还是回归中好的和不好的控制变量的作用。