发帖
雷达卡
towardsdatascience.com/causal-machine-learning-what-can-we-accomplish-with-a-single-theorem-ea81432f1fac因果推断基础：从ATE到CATE
因果推断的核心任务是量化处理变量（T）对结果变量（Y）的影响。典型问题包括：运动是否有助于减重？营销活动能否提高转化率？价格调整如何影响销量？某种疗法是否改善了患者的健康状况？ 当处理变量T是通过随机化方式分配时，如在随机对照试验（RCT）中，我们可以直接通过比较不同处理组的结果差异来估计因果效应。此时，无需额外控制其他变量，也能在总体层面得到无偏估计。常见的做法是使用线性回归模型预测Y，其中T的回归系数θ即代表处理每增加一个单位时Y的平均变化量：
在这种设定下，若T为随机分配，则θ被解释为**平均处理效应**（Average Treatment Effect, ATE），其一般形式定义如下：
然而，在大多数实际场景中，数据来源于自然发生的观测行为，而非实验设计。这意味着处理的分配可能受到多种潜在因素的影响，从而引入**混杂变量**——那些同时影响处理选择和结果输出的变量。如果不加以控制，这些变量会导致因果效应估计出现偏差。
以饮食习惯与预期寿命的关系为例：数据显示素食者通常寿命更长。但这种关联未必说明素食本身延长寿命，因为素食者往往整体生活方式更为健康——他们更注重锻炼、定期体检、避免高风险行为等。这些“健康意识”相关的特征就是典型的混杂变量，它们同时驱动了饮食选择和长寿结果。
除了混杂变量外，还有一类重要变量称为**背景变量**，它们不引起处理分配（尤其在随机实验中），但却调节处理对结果的作用强度。例如，在一次随机开展的营销活动中，虽然人群是否接收到广告是随机决定的，但年龄、收入水平、性别或过往购买记录等因素仍会影响个体对广告的响应程度。这类变量虽非混杂因子，却决定了因果效应的异质性。
我们将上述两类变量统称为**干扰参数**，并用协变量矩阵X表示。尽管我们关注的核心仍是T→Y的因果路径，但为了准确估计这一路径的强度，必须考虑X的影响。
因此，在现代因果机器学习中，研究者更多地聚焦于估计**条件平均处理效应**（Conditional Average Treatment Effect, CATE），即在给定一组特定干扰参数x的条件下，处理带来的期望结果变化：
当假设Y与T及X之间存在线性关系时，可以通过在回归模型中引入T与X的交互项来估计CATE。这种方法允许因果效应随个体特征的变化而变化，从而实现个性化因果推断。
towardsdatascience.com/causal-machine-learning-what-can-we-accomplish-with-a-single-theorem-ea81432f1fac一个可用于估计条件平均处理效应（CATE）的线性模型表达式。
towardsdatascience.com/causal-machine-learning-what-can-we-accomplish-with-a-single-theorem-ea81432f1fac利用上述线性形式，我们可以通过对结果变量 Y 关于处理变量 T 求导，进而推导出 CATE 的数学表达式：
linear_example1.csv该公式展示了在线性框架下如何计算 CATE。
statsmodels尽管线性回归是一种基础方法，但近年来已发展出多种更复杂的机器学习算法用于 CATE 估计。这类方法通常被称为“元学习器”，其中双重机器学习（也称 R-学习器）是典型代表之一。
在因果分析中，若变量集合 X 能够在被控制后消除所有混杂偏差，从而准确识别处理变量 T 对结果 Y 的因果影响，则称 X 为“充分调整集”。为了实现可靠的因果推断，全面收集并纳入潜在的混杂因素至关重要。忽略关键混杂变量可能导致严重的估计偏误，进而得出错误结论。
虽然因果推断与机器学习领域内容广泛，但本文所涉及的基础知识已足以支撑后续讨论。接下来，我们将深入探讨计量经济学中的一个核心定理——弗里施-沃格-洛夫（FWL）定理，以揭示线性回归背后的本质机制。这不仅有助于理解为何线性回归在因果建模中如此有效，也为后续复杂方法的学习奠定基础。
弗里施-沃格-洛夫（FWL）定理
FWL 定理是计量经济学中的一个重要理论工具，它深刻揭示了线性回归的工作原理，并解释了其在因果推断中的强大表现。尽管本文将从因果视角解读该定理，但它实际上反映了线性回归的一般性质。
假设我们希望评估处理变量 T 对结果变量 Y 的影响，同时控制一组协变量 X。我们假定这种关系可通过如下线性模型刻画：
t_resid该方程表示将因变量 Y 同时对处理变量 T 和协变量矩阵 X 进行回归。
y_resid按照标准做法，我们采用普通最小二乘法（OLS）来估计系数向量 θ 和 θ。在因果推断中，我们重点关注 θ，因为它反映了在控制 X 后 T 对 Y 的净效应。但“控制 X”究竟意味着什么？FWL 定理为此提供了清晰的几何与代数解释。
根据 FWL 定理，θ 可通过以下步骤等价获得：首先，分别使用 X 对 Y 和 T 进行线性回归，得到各自的预测值：
y_resid其中，Ylr 和 Tlr 分别表示基于 X 预测出的 Y 和 T 的拟合值。随后，在第二步中，我们对残差进行回归：
t_resid即，将 Y 在去除 X 影响后的残差（Y – Ylr）对 T 去除 X 影响后的残差（T – Tlr）做简单回归，所得斜率即为原始多元回归中的 θ。
初看之下，这种方法似乎多此一举——既然一次回归即可得到结果，为何要拆分为三步？然而，FWL 定理的价值在于其提供的直观理解：它揭示了线性回归本质上是在“去偏”之后衡量变量间的关系。
进一步解析各残差项的意义：
- Y – Ylr：表示在剔除了所有能由 X 线性解释的部分后，Y 中剩余的不可预测成分。如果 X 能完美预测 Y，则该残差为零，说明 T 对 Y 没有额外解释力；从因果角度看，这意味着混杂因素已充分解释结果变异，治疗的作用可能不显著。
- T – Tlr：表示在排除了 X 对 T 的影响后，处理变量中“纯净”的变动部分。这一残差帮助我们剥离由混杂因素引起的非随机分配偏差。若该值为零，表明 T 完全由 X 决定，无法进行有效的因果推断。反之，当 X 构成充分调整集且满足一定假设时，T – Tlr 可视为一种“类随机化”的处理指标。
当我们对这两个残差进行回归时，实际上是在一个已被“清洗”过的空间中评估 T 对 Y 的独立影响。由于 X 的影响已被分别从 Y 和 T 中移除，最终回归结果有效地隔离了处理变量的真实因果效应。
解释方差的原理正是我们能够将系数 θ 理解为：当自变量 T 增加一个单位，而其他变量 X 保持不变时，因变量 Y 的期望变化量的原因所在。
为了更深入理解理想化的弗里希-沃-洛夫（FWL）定理所表达的核心思想，接下来我们将通过两个在 Python 中实现的具体示例进行说明。
示例一：存在混杂因素的处理效应分析
本示例使用一个基于以下生成机制构建的合成数据集：
示例1的数据生成过程。图片由作者提供。
我们的目标是评估变量 t 对 y 的因果影响。然而，t 与 y 之间的关系受到 x1 和 x2 这两个混杂变量的影响。从第一个方程可以看出，t 可以表示为 x1 和 x2 的线性组合加上高斯噪声项，这表明 t 与这些协变量相关联。
此外，相较于 x1 和 x2，t 对 y 的预测能力较弱，这一点可以从其较小的回归系数中看出。尽管如此，我们关注的重点仅在于 t 对 y 的净效应，而不考虑其效应强度。
以下是本部分所需的程序依赖项：
DATA_PATH = "data/"
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
from lightgbm.sklearn import LGBMRegressor
from sklearn.model_selection import cross_val_predict
接下来加载数据集：
data_linear_1 = pd.read_csv(DATA_PATH + "linear_example1.csv")
print(data_linear_1.shape)
"(10000, 4)"
print(data_linear_1.head())
"""
t x1 x2 y
0 31.139704 1.101262 0.682916 874.357591
1 -46.437301 0.338431 0.064071 -181.341797
2 -94.685877 -0.539972 1.060375 2190.136391
3 -30.234819 -1.260242 1.994309 4525.823182
4 43.309347 -1.894621 -0.086365 870.204667
"""
linear_example1.csv该数据集包含 10,000 条观测记录，涵盖变量 t、x1、x2 和 y。下图展示了 t 与 y 之间的散点关系：
t 与 y 的散点图。图片由作者提供。
从图中可以观察到，t 与 y 之间似乎没有明显的线性趋势或关联模式。我们可以通过对 t 和 y 进行简单线性回归来验证这一直观判断：
t_ols_model = smf.ols("y ~ t", data=data_linear_1)
print(t_ols_model.fit().summary().tables[1])
"""
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept -11.8412 20.760 -0.570 0.568 -52.535 28.853
t 0.2665 0.418 0.637 0.524 -0.553 1.086
==============================================================================
"""
statsmodels该模型将 y 对 t 进行回归。结果显示，t 的估计系数为 0.2665，远低于真实值 2。同时，其 p 值为 0.524，大于常用显著性水平（如 0.05），说明无法拒绝“t 对 y 无影响”的原假设。若忽略潜在混杂因素的存在，我们可能错误地推断 t 对 y 没有显著作用。
当我们将 x1 和 x2 同时纳入回归模型中，对 t 的系数估计将发生显著变化。此时，t 的系数估计值为 1.9998，相较于真实值 2 已非常接近，说明模型的拟合效果更优。此外，该系数对应的 p 值为 0，表明其在统计上高度显著。
从实际意义上看，在控制了 x1 和 x2 的影响后，t 每增加一个单位，y 的平均值预计会上升约 1.9998 个单位。若 x1 与 x2 构成了充分的调整集，则可推断出 t 对 y 的平均处理效应（ATE）约为 1.9998。
这种改进源于 x1 和 x2 是 t 与 y 关系中的混杂变量。为了准确识别 t 对 y 的因果效应，必须排除这些变量带来的干扰。线性回归通过控制协变量的方式实现了这一点，从而提供了一个更为精确的因果估计。
根据弗里施-瓦乌-洛夫利（FWL）定理，我们可以通过一种等价的分解方式获得相同的 t 系数。具体操作如下：首先建立两个辅助回归模型——一个用 x1 和 x2 预测 t，另一个用 x1 和 x2 预测 y：
t_resid_model = smf.ols("t ~ x1 + x2", data=data_linear_1).fit()
y_resid_model = smf.ols("y ~ x1 + x2", data=data_linear_1).fit()
接着，利用这两个模型计算残差，得到去除 x1 和 x2 线性影响后的 t 和 y 的残差版本：
data_linear_1["t_resid"] = data_linear_1["t"].values - t_resid_model.predict(data_linear_1) data_linear_1["y_resid"] = data_linear_1["y"].values - y_resid_model.predict(data_linear_1)
t_resid其中，
y_resid表示的是剔除了 x1 和 x2 可解释部分之后的 y 值。如果 x1 和 x2 是一个充分的调整集，那么
y_resid接下来，我们对残差化的变量进行简单线性回归：
resid_linear_model = smf.ols("y_resid ~ t_resid", data=data_linear_1)
print(resid_linear_model.fit().summary().tables[1])
这种方法所得到的 t_resid 系数与完整模型中 t 的系数完全一致，验证了 FWL 定理的有效性：无论采用全模型回归还是分步残差回归，只要调整集充分，就能获得相同的因果效应估计结果。
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 0.1830 0.504 0.363 0.716 -0.805 1.171
t 1.9998 0.010 196.444 0.000 1.980 2.020
x1 -499.9225 0.505 -990.130 0.000 -500.912 -498.933
x2 2000.2940 0.500 4000.091 0.000 1999.314 2001.274
==============================================================================
Intercept 2.526e-13 0.504 5.02e-13 1.000 -0.987 0.987 t_resid 1.9998 0.010 196.464 0.000 1.980 2.020
根据 FWL 定理，我们可以观察到，在对变量进行残差化处理后，t 的回归系数为 1.9998，这与在完整模型中估计出的系数完全一致。同时，p 值和置信区间也保持不变。通过从 t 和 y 中分别剔除 x1 和 x2 的影响，我们利用简单线性回归准确地还原了 t 对 y 的真实效应。
为了更直观地理解这一过程，我们可以绘制以下两个残差之间的关系图：
t_residy_resid上图展示了 t 与 y 残差项之间的关联情况（图像来源：作者提供）。相比原始数据中 t 与 y 的散点分布，这种残差化后的图形呈现出明显的线性趋势，而原始数据之间的关系则显得较为嘈杂。通过去除 x1 和 x2 这些混杂因素的影响，我们成功揭示了 t 与 y 之间潜在的真实线性关系。
示例 2：虚假相关性的识别
接下来考虑另一种情形，我们使用如下生成机制构造数据集（见下图）：
https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/05be33be814db304721bfc26a215ebb5.png该图为示例 2 的数据生成方程。（图片由作者提供）
注意，在此设定中，t 是 x1 的线性函数加上高斯噪声，但在 y 的生成方程中并未包含 t，这意味着 t 实际上并不影响 y。因此，任何观测到的 t 与 y 之间的关联都属于虚假相关。
加载该示例的数据集：
data_linear_2 = pd.read_csv(DATA_PATH + "linear_example2.csv") print(data_linear_2.shape) "(10000, 4)" print(data_linear_2.head())
t x1 x2 y
0 44.758318 1.101262 0.682916 -44.527435
1 12.625887 0.338431 0.064071 -63.033848
2 -23.450555 -0.539972 1.060375 36.979009
3 -50.945122 -1.260242 1.994309 70.451553
4 -74.999631 -1.894621 -0.086365 102.008320
首先绘制原始的 t 对 y 的散点图以观察其表观关系：
这是 t 与 y 的原始散点图。（图片由作者提供）
从图中可见，t 与 y 似乎存在较强的负向线性关系。为进一步验证，我们对 y 关于 t 进行简单线性回归分析：
t_ols_model = smf.ols("y ~ t", data=data_linear_2)
print(t_ols_model.fit().summary().tables[1])
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 0.0916 0.544 0.169 0.866 -0.974 1.158
t -1.2452 0.014 -91.498 0.000 -1.272 -1.219
==============================================================================
结果显示，t 的回归系数显著不为零（约为 -1.245），且 p 值极小，表明统计上高度显著。然而，这一结果具有误导性——因为实际上 t 并未参与 y 的生成过程。这种“显著”关系源于共同驱动 t 和 y 的潜变量（如 x1），从而造成了虚假相关现象。
在进行回归分析时，我们在变量 t 上得到了一个统计显著的负系数（-1.2452）。如果忽略潜在混杂因素的影响，我们可能会误认为 t 对 y 存在因果关系。这正是“相关性不等于因果性”的经典案例。
为了更准确地理解这种关系，我们可以跳过对 t、x1 和 x2 与 y 的整体回归，转而构建 t 和 y 的残差模型：
t_resid_model = smf.ols("t ~ x1 + x2", data=data_linear_2).fit()
y_resid_model = smf.ols("y ~ x1 + x2", data=data_linear_2).fit()
data_linear_2["t_resid"] = data_linear_2["t"].values - t_resid_model.predict(data_linear_2)
data_linear_2["y_resid"] = data_linear_2["y"].values - y_resid_model.predict(data_linear_2)
经过上述处理后，我们得到的是去除 x1 和 x2 线性影响后的 t 与 y 的残差值。下图展示了这些残差的分布情况：
从图中可以看出，在剔除了 x1 和 x2 所带来的方差之后，t 与 y 残差之间几乎不存在任何明显的关系。此前观察到的相关性实际上是由于 t 与 x1 高度相关，而 x1 又直接影响 y 所致。一旦控制了这两个协变量，t 与 y 之间的虚假关联便消失了，揭示出它们之间并无直接联系。
为验证这一点，我们进一步对两个残差变量进行回归分析：
resid_linear_model = smf.ols("y_resid ~ t_resid", data=data_linear_2)
print(resid_linear_model.fit().summary().tables[1])
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 8.59e-15 0.504 1.71e-14 1.000 -0.987 0.987
t_resid -0.0122 0.509 -0.024 0.981 -1.010 0.985
==============================================================================
回归结果表明，在调整了 x1 和 x2 的影响后，t_resid 对 y_resid 几乎没有可检测到的作用。系数的 p 值接近 1，且置信区间包含零，说明该效应不显著。
t_resid这一结果支持了我们的判断：在排除混杂因素后，t 对 y 并无实际影响。即：
ty通过这一过程，我们实际应用了弗里施-瓦乌法-洛夫利（FWL）定理，并加深了对线性回归机制的理解。核心在于：通过对 t 和 y 分别关于混杂变量 x1 和 x2 进行残差化处理，可以有效剥离混杂路径，从而揭示二者之间的真实关系。
接下来的问题是：当数据呈现非线性特征时，是否仍可应用 FWL 定理的思想？我们能否利用类似方法来识别因果效应？以下将探讨这一问题。
将 FWL 定理推广至非线性数据场景
尽管 FWL 定理严格适用于线性回归框架，但现实中的数据往往具有复杂的非线性结构。幸运的是，我们可以借鉴 FWL 定理背后的逻辑直觉，并将其扩展到更广泛的建模情境中，以帮助识别因果关系。
部分线性模型的应用
虽然理想情况下我们希望避免对结果变量与处理变量及协变量之间的函数形式做出强假设，但一个自然的延伸方向是探索：在存在非线性混杂效应的情况下，是否仍能估计出线性的处理效应？
为此，我们考虑如下模型设定：
该公式表示因变量 Y 是处理变量 T 的线性函数，加上协变量 X 的某个任意函数（可能为非线性），同时允许 T 本身也受 X 的非线性影响。这类模型被称为部分线性模型，它结合了灵活性与可解释性，成为连接传统线性回归与现代因果推断方法的重要桥梁。
部分线性模型因其同时包含线性和非线性成分而具有较强的表达能力。尽管其对部分结构仍存在一定的假设限制，但相比纯线性模型，它在建模灵活性上有了显著提升，使我们能够更逼近真实的数据生成机制。
要拟合一个部分线性模型，可以借助现代机器学习方法的强大非线性拟合能力。例如，基于树的模型或其他任意非线性算法可用于估计条件期望。我们可以用机器学习模型来预测结果变量 Y 关于协变量 X 的期望值，记为 Yml = E[Y | X]；同样地，也可以用另一个模型预测处理变量 T 关于 X 的期望值，即 Tml = E[T | X]。为了消除 X 对 T 和 Y 所带来的干扰，可借鉴弗里施-沃夫鲁-洛夫（FWL）定理的思想，将残差重新定义如下：
towardsdatascience.com/causal-machine-learning-what-can-we-accomplish-with-a-single-theorem-ea81432f1fac该图展示了利用任意机器学习模型构建的残差形式。图片由作者提供。
回顾经典的 FWL 定理，其中 E[T | X] 与 E[Y | X] 通常由线性回归估计。但在当前设定中，我们可以使用任意机器学习模型来近似这两个条件期望，从而允许捕捉协变量影响中的非线性模式。
在假设 Tml 和 Yml 能够准确估计真实条件期望的前提下，参数 θ 可通过以下回归方式恢复：
linear_example1.csv此图为部分线性模型中残差化方法的示意图。图片由作者提供。
这意味着，我们首先使用机器学习模型分别从 T 和 Y 中去除 X 的影响，得到对应的残差项；然后将 Y 的残差对 T 的残差进行简单线性回归，所得斜率即为因果效应或目标参数 θ 的一致估计量。
为了验证这一方法的有效性，我们采用如下数据生成机制构造模拟数据集：
statsmodels上述公式描述了部分线性模型下数据的生成过程。其中，t 是 x 和 x 的非线性函数，而 y 是 t 的线性函数加上 x 与 x 的某个未知非线性函数。
接下来加载该数据集并查看基本信息：
data_partial_linear = pd.read_csv(DATA_PATH + "partial_linear_example.csv")
print(data_partial_linear.shape)
"(10000, 4)"
print(data_partial_linear.head())
"""
t x1 x2 y
0 2.052069 1.101262 0.682916 145.112924
1 0.363920 0.338431 0.064071 81.352175
2 -1.508082 -0.539972 1.060375 -37.170164
3 -1.182025 -1.260242 1.994309 -119.741583
4 -0.481000 -1.894621 -0.086365 11.796844
"""
为进一步理解变量间关系，绘制 t 与 y 的散点图：
t_resid该散点图反映了在部分线性数据集中 t 与 y 的关系。图片由作者提供。
可以看出，两者之间呈现出一种模糊且带有噪声的关联趋势。为量化这种关系，我们先对 y 关于 t 进行普通最小二乘回归：
t_ols_model = smf.ols("y ~ t", data=data_partial_linear)
print(t_ols_model.fit().summary().tables[1])
"""
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 27.1079 0.598 45.330 0.000 25.936 28.280
"""
============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ Intercept 34.2972 0.442 77.676 0.000 33.432 35.163 t 4.9527 0.349 14.188 0.000 4.268 5.637 x1 34.6709 0.443 78.290 0.000 33.803 35.539 x2 -0.0551 0.388 -0.142 0.887 -0.816 0.706 x1:x2 21.1510 0.386 54.814 0.000 20.395 21.907 ==============================================================================完整的线性回归模型相较于简单模型更接近真实的 t 系数值，但依然存在一定程度的高估。值得注意的是，在模型构建过程中引入了 x1 与 x2 的交互项（x1:x2），这是因为我们并未完全掌握 y 的真实函数形式，因此需要考虑变量间可能存在的相互作用对结果的影响。 回顾此前的结果：
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
t 17.8144 0.428 41.616 0.000 16.975 18.653
可以看出，简单线性模型严重高估了 t 的影响系数。这提示我们仅依赖基础线性假设可能会导致误导性的结论。
接下来，假设我们并不知道 y 实际上是关于 x1、x2 和 t 的部分线性函数，而是误以为其整体呈线性关系。在此前提下，通过完整线性回归得到 t 的估计系数约为 4.95。为了检验该线性假设是否成立，可以借助残差分析方法，具体做法如下：
首先拟合两个辅助回归模型：
- 使用 x1、x2 及其交互项预测 t：
`t_resid_model = smf.ols("t ~ x1 + x2 + x1*x2", data=data_partial_linear).fit()`
- 同样结构预测 y：
`y_resid_model = smf.ols("y ~ x1 + x2 + x1*x2", data=data_partial_linear).fit()`
随后计算残差：
```python
data_partial_linear["t_resid"] = data_partial_linear["t"].values - t_resid_model.predict(data_partial_linear)
data_partial_linear["y_resid"] = data_partial_linear["y"].values - y_resid_model.predict(data_partial_linear)
```
这些残差反映了在控制 x1 和 x2 影响后，t 与 y 中无法被线性部分解释的成分。我们利用扰动后的协变量（x1 和 x2）作为特征来建模 t 和 y，进而提取出残差序列，用于后续分析。
towardsdatascience.com/causal-machine-learning-what-can-we-accomplish-with-a-single-theorem-ea81432f1fac我们首先定义扰动变量列表：
nuisance = ["x1", "x2"]
接着，构建两个 LightGBM 回归模型，分别用于预测处理变量 t 和结果变量 y：
model_t = LGBMRegressor(n_estimators=800, max_depth=6)
model_y = LGBMRegressor(n_estimators=800, max_depth=6)
利用交叉验证方法，对 t 和 y 基于扰动变量进行预测，并将预测值存入数据框中：
data_partial_linear["t_hat"] = cross_val_predict(model_t, data_partial_linear[nuisance], data_partial_linear["t"], cv=5)
data_partial_linear["y_hat"] = cross_val_predict(model_y, data_partial_linear[nuisance], data_partial_linear["y"], cv=5)
随后计算残差，即实际值与预测值之间的差异：
data_partial_linear["t_resid"] = data_partial_linear["t"] - data_partial_linear["t_hat"]
data_partial_linear["y_resid"] = data_partial_linear["y"] - data_partial_linear["y_hat"]
LGBMRegressor该表达式用于估计 t 与 y 在扰动参数影响下的表现形式。接下来，结合以下公式
cross_val_predict并采用 5 折交叉验证策略，将 t 和 y 的预测结果添加至
data_partial_linear确保每个样本的预测来源于未参与训练的外部数据，从而避免过拟合。最终，在
data_partial_linear中完成残差的计算与存储。
以下是残差分布的可视化结果：
t 与 y 的线性模型残差图。图片由作者提供。
从图中可以看出，尽管仍存在一定噪声和异常点，但在使用 LightGBM 模型有效去除由 x1 和 x2 非线性关系带来的方差后，t 与 y 之间的线性趋势变得更加清晰。
为了进一步验证这一观察，我们对残差变量进行简单线性回归分析：
resid_linear_model = smf.ols("y_resid ~ t_resid", data=data_partial_linear).fit()
print(resid_linear_model.summary().tables[1])
"""
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept -0.0246 0.039 -0.631 0.528 -0.101 0.052
t_resid 1.9824 0.037 53.344 0.000 1.910 2.055
==============================================================================
"""
回归结果显示了高度显著的关系。其中 t_resid 的系数约为 1.98，相较于多重线性回归的估计更接近真实值 2。这表明我们通过应用 Frisch-Waugh-Lovell（FWL）定理的思想，成功剔除了扰动变量非线性效应所带来的干扰，进而更准确地揭示了 t 与 y 之间潜在的线性关联。
部分线性模型的优势在于，它允许我们使用普通最小二乘法或贝叶斯回归等传统统计方法，对 T 与 Y 之间的关系做出具有统计意义的推断。通过选用灵活的机器学习模型来消除混杂变量的影响，我们将原本复杂的因果建模问题转化为一个标准的线性回归任务，从而能够方便地获取 T 系数的 p 值、置信区间或可信区间。
然而，部分线性模型的假设在某些场景下可能过于理想化。例如，若无法保证 T 与 Y 之间存在线性关系，或怀疑 T 与协变量 X 存在非线性交互作用，则该模型可能不足以准确捕捉真实的因果机制。
我们的最终目标是：当 Y 是 T 和 X 的完全未知函数（即 Y = f(T, X)）时，仍能有效识别 T 对 Y 的因果效应。
信仰的飞跃——非线性结果函数
在本节末尾，我们将彻底放弃对结果函数形式的任何先验假设（即不再假设 Y 与 T, X 之间具有特定函数关系），尝试直接估计 T 与 Y 之间的因果联系。这项任务极具挑战性。
因为 T 与 Y 的因果路径可能极为复杂，并且很可能受到 X 的调节影响，尤其是在存在非线性交互或多层嵌套依赖的情况下。如何在这种高维、非线性的环境中稳健地识别因果效应，将是接下来探讨的核心问题。
尽管函数的近似评估存在难度，尤其是在涉及异质性关系时，但从经验角度来看，这种近似仍然是可行的。我们接下来将探讨如何在不依赖同质线性假设的前提下，对部分线性模型进行扩展与优化。
首先回顾残差化部分线性模型的基本形式：
towardsdatascience.com/causal-machine-learning-what-can-we-accomplish-with-a-single-theorem-ea81432f1fac残差化部分线性模型方程。图片由作者提供。
在传统的部分线性框架中，我们利用机器学习方法分别估计条件期望 E[Y | X] 和 E[T | X]，从而提取出 Y 与 T 的残差，以消除混杂变量 X 带来的方差影响。然而，该模型仍假设处理变量 T 与结果变量 Y 之间存在全局一致的线性关系。我们可以进一步放宽这一限制。
为了实现这一点，我们引入一个更为灵活的建模策略：使用从上述步骤中获得的 T 残差和协变量 X，来预测 Y 的残差。具体而言，我们不再假设残差间的关系是线性的或同质的，而是通过一个非参数化的机器学习模型来捕捉其潜在的复杂结构：
linear_example1.csv使用 T 残差和 X 预测 Y 残差的模型。图片由作者提供。
在这个设定中，我们将 Y 残差视为 T 残差与特征向量 X 的某种函数 f 的输出，即：
Y_residual = f(T_residual, X)
其中，包含 X 的目的是允许处理效应与协变量之间存在交互作用，从而体现异质性因果效应。
已有模型用于生成 T 和 Y 的残差，而第三个独立的机器学习模型（如梯度提升树、神经网络等）可用于拟合函数 f。一旦完成训练，该模型即可用于反事实推理——例如，在保持 X 不变的情况下，改变 T 的取值，观察预测的 Y 残差如何变化。这为我们提供了对干预效果进行精细化分析的能力。
正如 Matheus Facure 在其著作《Python 因果推断手册》中所强调的，这类外推虽具有强大潜力，但也伴随风险：若真实数据中的处理机制复杂且未被充分建模，则模型可能产生误导性预测。然而，若对函数 f 的估计足够稳健，此方法可构建出高度实用的因果推断系统，适用于多种处理方案下的影响评估。
为验证该方法的有效性，我们考虑如下数据生成过程：
statsmodels图例 3 的数据生成方程。图由作者提供。
在此模拟场景中，处理变量 t 是协变量 x1 和 x2 的函数，而结果变量 y 同时依赖于 t、x1、x2 及 x3。因此，x1 和 x2 被视为混杂因素，而 x3 并非混杂变量，但可能通过与 t 的交互作用间接影响 y。
加载并查看数据集：
data_non_linear = pd.read_csv(DATA_PATH + "non_linear_example.csv")
print(data_non_linear.shape)
"(15000, 5)"
print(data_non_linear.head())
"""
t x1 x2 x3 y
0 0.815358 1.101262 -0.378001 7 -0.448740
1 -4.789706 0.338431 -0.585686 7 -0.756166
2 -0.201015 -0.539972 0.964218 8 -0.861377
3 -4.267054 -1.260242 0.105409 2 1.447061
4 -2.871802 -1.894621 -1.388249 5 -0.199239
"""
print(data_non_linear["x3"].nunique())
"9"
注意到 x3 是一个具有 9 个不同整数值的离散变量，这一特性将在后续建模中发挥重要作用。
为进一步理解变量间关系，我们绘制原始 t 与 y 的散点图：
t_residt 对 y。图由作者提供。
从图中可见，t 与 y 之间的关系并不明显，可能存在非线性或异质性效应。为剥离混杂因素的影响，我们采用 LightGBM 分别对 T 和 Y 进行残差化处理，所选协变量包括 x1、x2 和 x3：
nuisance = ["x1", "x2", "x3"]
model_t = LGBMRegressor(n_estimators=300, max_depth=6)
model_y = LGBMRegressor(n_estimators=300, max_depth=6)
data_non_linear["t_hat"] = cross_val_predict(
model_t, data_non_linear[nuisance], data_non_linear["t"], cv=5
)
data_non_linear["y_hat"] = cross_val_predict(
model_y, data_non_linear[nuisance], data_non_linear["y"], cv=5
)
data_non_linear["t_resid"] = data_non_linear["t"] - data_non_linear["t_hat"]
data_non_linear["y_resid"] = data_non_linear["y"] - data_non_linear["y_hat"]
让我们观察一下这些残差的分布情况：
towardsdatascience.com/causal-machine-learning-what-can-we-accomplish-with-a-single-theorem-ea81432f1fac
data_non_linear["t_resid_squared"] = data_non_linear["t_resid"] ** 2
resid_quadratic_model = smf.ols("y_resid ~ t_resid_squared", data=data_non_linear).fit()
print(resid_quadratic_model.summary())
回归结果如下所示：
OLS Regression Results
==============================================================================
Dep. Variable: y_resid R-squared: 0.547
Model: OLS Adj. R-squared: 0.547
Method: Least Squares F-statistic: 1.808e+04
Date: Tue, 19 Mar 2024 Prob (F-statistic): 0.00
Time: 09:54:06 Log-Likelihood: 13503.
No. Observations: 15000 AIC: -2.700e+04
Df Residuals: 14998 BIC: -2.699e+04
Df Model: 1
Covariance Type: nonrobust
===================================================================================
coef std err t P>|t| [0.025 0.975]
-----------------------------------------------------------------------------------
Intercept 0.0758 0.001 77.241 0.000 0.074 0.078
t_resid_squared -0.0028 2.12e-05 -134.461 0.000 -0.003 -0.003
==============================================================================
Omnibus: 3204.283 Durbin-Watson: 2.004
从结果可以看出，t_resid_squared 的系数显著不为零（p < 0.001），且模型解释了约 54.7% 的残差变异，说明二次项在拟合残差关系中起到了重要作用。这进一步支持了原始变量间存在非线性关联的假设。
towardsdatascience.com/causal-machine-learning-what-can-we-accomplish-with-a-single-theorem-ea81432f1fac我们观察到，残差平方与结果残差之间呈现出较为显著的R-squared值。这表明，在排除了干扰参数的影响后，处理变量对结果可能存在明显的二次效应。这一发现为我们进一步分析提供了有力支持。
为了更直观地展示这种关系，可以将预测的二次多项式与残差进行可视化绘图：
towardsdatascience.com/causal-machine-learning-what-can-we-accomplish-with-a-single-theorem-ea81432f1fac上图展示了将 t 和 y 的残差拟合至一个二次多项式的结果，图像由作者提供。可以看出，二次模型在当前数据下表现出良好的拟合效果。然而，这仅是分析的第一步。
需要注意的是，处理变量可能与多个干扰参数（如 x1、x2 和 x3）存在交互作用。这意味着 t 与 y 之间的关系并非固定不变，而是依赖于这些干扰变量的具体取值。此前我们构建的模型仅反映了 y 残差作为 t 残差的函数关系；而我们的最终目标是建立 y 残差同时依赖于 t 残差和各干扰参数的联合函数模型。
为此，我们可以利用已计算出的 t 残差和 y 残差，并引入原始干扰参数，来训练一个新的非线性模型：
data_non_linear["y_resid_2"] = cross_val_predict(
LGBMRegressor(n_estimators=300, max_depth=6),
data_non_linear[nuisance + ["t_resid"]],
data_non_linear["y_resid"],
cv=5,
)
接下来，我们再次检验该二次项对新生成的 y 残差的拟合情况：
resid_quadratic_model_2 = smf.ols(
"y_resid_2 ~ t_resid_squared", data=data_non_linear
).fit()
print(resid_quadratic_model_2.summary())
输出结果如下：
OLS Regression Results
==============================================================================
Dep. Variable: y_resid_2 R-squared: 0.776
Model: OLS Adj. R-squared: 0.776
Method: Least Squares F-statistic: 5.190e+04
Date: Tue, 26 Mar 2024 Prob (F-statistic): 0.00
Time: 08:30:40 Log-Likelihood: 21629.
No. Observations: 15000 AIC: -4.325e+04
Df Residuals: 14998 BIC: -4.324e+04
Df Model: 1
Covariance Type: nonrobust
===================================================================================
统计结果显示，调整后的R方仍维持在较高水平，说明即使在考虑干扰因素后，二次项依然具有较强的解释能力。此外，F检验的p值趋近于零，进一步验证了模型整体的显著性。
配套信息提示：
- [1] 标准误差的计算基于误差协方差矩阵被正确设定的前提。
相关统计指标汇总：
| 左侧指标 | 数值 | 右侧指标 | 数值 | |
|---|---|---|---|---|
| Prob(Omnibus) | 0.000 | Jarque-Bera (JB) | 72693.392 | |
| Skew | -0.461 | Prob(JB) | 0.00 | |
| Kurtosis | 13.745 | Cond. No. | 56.7 |
综上所述，通过逐步控制干扰变量并引入残差建模策略，我们增强了对处理变量非线性影响的理解。后续分析可在此基础上扩展更高阶多项式或采用其他非参数方法以提升模型表达力。
coef std err t P>|t| [0.025 0.975] ----------------------------------------------------------------------------------- Intercept 0.0743 0.001 130.244 0.000 0.073 0.075 t_resid_squared -0.0028 1.23e-05 -227.816 0.000 -0.003 -0.003 ============================================================================== Omnibus: 2487.996 Durbin-Watson: 1.992 Prob(Omnibus): 0.000 Jarque-Bera (JB): 35598.209 Skew: -0.349 Prob(JB): 0.00 Kurtosis: 10.515 Cond. No. 56.7在对数据进行建模时，将二次多项式应用于
t_resid_squaredy_resid_2t_resid_squaredy_residcross_val_predict()t_modely_modelfinal_model
x1 x2 x3 t
0 0.003099 -0.014613 1 -10.000000
1 0.003099 -0.014613 1 -9.889447
2 0.003099 -0.014613 1 -9.778894
3 0.003099 -0.014613 1 -9.668342
4 0.003099 -0.014613 1 -9.557789
这些测试样本可用于后续的残差计算与反事实预测流程。注意，在此测试数据集中，变量 x1 和 x2 均仅取单一固定值，即它们分别等于训练数据中的均值，分别为 0.00309858 和 -0.01461324。而变量 x3 则覆盖了训练集中的全部 9 个不同取值（从 1 到 9）。这样的构造方式使得我们能够在保持干扰变量（nuisance parameters）不变的前提下，观察预测响应如何随着变量 t 的变化而变化。 首先，利用已训练好的模型对测试数据进行初步预测： test_data["y_hat"] = y_model.predict(test_data[nuisance]) test_data["t_hat"] = t_model.predict(test_data[nuisance]) 这一步分别生成了基于干扰协变量对 y 和 t 的预测值 y_hat 与 t_hat。 接下来，计算 t 的残差，即实际 t 值与预测 t_hat 之间的差异： test_data["t_resid"] = test_data["t"] - test_data["t_hat"] 随后，将该残差项作为关键特征之一，输入到最终模型中以预测 y 的残差部分： test_data["y_resid"] = final_model.predict(test_data[nuisance + ["t_resid"]])这一过程的核心思想是剥离干扰因素影响后，捕捉处理变量 t 对结果 y 的净效应。 最后，将预测的 y_hat 与由 t 残差驱动的 y_resid 相加，得到反事实情境下的 y 预测值： test_data["y_counterfactual"] = test_data["y_hat"] + test_data["y_resid"] 完整输出如下所示： x1 x2 x3 t y_hat t_hat t_resid y_resid y_counterfactual 0 0.003099 -0.014613 1 -10.000000 1.55532 0.201618 -10.201618 -0.240874 1.314446 1 0.003099 -0.014613 1 -9.889447 1.55532 0.201618 -10.091065 -0.240874 1.314446 2 0.003099 -0.014613 1 -9.778894 1.55532 0.201618 -9.980512 -0.240874 1.314446 3 0.003099 -0.014613 1 -9.668342 1.55532 0.201618 -9.869960 -0.233701 1.321619 4 0.003099 -0.014613 1 -9.557789 1.55532 0.201618 -9.759407 -0.233701 1.321619 ... ... ... .. ... ... ... ... ... ... 1795 0.003099 -0.014613 9 11.557789 -1.34144 1.441577 10.116212 -0.224248 -1.565688 1796 0.003099 -0.014613 9 11.668342 -1.34144 1.441577 10.226765 -0.230682 -1.572122 1797 0.003099 -0.014613 9 11.778894 -1.34144 1.441577 10.337317 -0.230682 -1.572122 1798 0.003099 -0.014613 9 11.889447 -1.34144 1.441577 10.447870 -0.230682 -1.572122 通过上述步骤构建的 y_counterfactual 反映了在控制其他变量不变时，目标变量 t 对应的因果效应估计。final_model
现在我们已知，变量 x1 和 x2 被固定在其均值水平上，而 x3 则取 9 个不同的离散值。在这些 x1 与 x2 的条件下，我们可以观察并可视化预测的 t 与 y 之间的关系如何随着 x3 的不同取值发生变化。
下图展示了在保持混杂因素不变的前提下，对不同 t 值进行反事实 y 预测的结果：
在保持干扰参数不变的情况下，对不同的 t 值进行反事实 y 预测。图片由作者提供。
从图中可以看出，在 x3 取不同值时，尽管 x1 和 x2 保持恒定，模型仍然预测出 t 与 y 之间存在差异化的关联模式。这表明模型成功捕捉到了处理效应中的非线性异质性特征。此外，残差中呈现出的二次型关系也依赖于干扰变量的变化。
接下来是其余 x3 取值下，预测的 t 与 y 关系的表现形式：
在保持干扰参数不变的情况下，对不同的 t 值进行反事实 y 预测。图片由作者提供。
这一系列结果的核心意义在于：
final_model尽管上述方法可用于生成反事实预测，但在实际应用中其评估难度较高。因此，实践中更常采用双机器学习框架中的 R-learner 方法。直观来说，R-learner 可被视为函数 f(T, X) 对 T 的偏导数在点 (x, t) 处的一种近似，而这正是条件平均处理效应（CATE）的定义。我们通过最小化如下损失函数来训练该模型：
R 损失函数。图片由作者提供。
其中，tml 和 yml 分别代表用于估计 T 和 Y 作为协变量 X 函数的机器学习模型的预测输出；而 τ(X, T) 是我们正在优化的目标模型，用于预测在样本点 (xi, ti) 上的 CATE。虽然该模型结构相对复杂，但具备更强的可评估性，并延续了自部分线性模型以来广泛使用的基于机器学习生成的残差策略。
总体而言，无论是旨在生成反事实预测，还是估计个体层面的处理效应（CATE），我们都可以借鉴弗里施-瓦乌权（FWL）定理的思想，构造相应的残差项，以揭示非线性混杂因素对因果关系的影响机制。
总结与注意事项
本文通过结合 FWL 定理的理论洞见，深入探讨了双重机器学习（DML）的基本原理。借助合成数据示例，我们验证了使用机器学习模型替代传统线性回归计算残差，有助于揭示复杂的非线性因果结构。
然而，必须保持警惕：深度机器学习模型及其衍生技术虽然功能强大，但也是一把双刃剑。它们虽能在控制高维混杂变量和挖掘潜在因果路径方面表现出色，但同时也要求研究者严格检验模型假设，并对输出结果进行审慎解读。若忽视内在机制而盲目依赖模型预测，极易导致错误推论。
现实世界的因果结构往往极为复杂，在面对高维干扰参数时，我们很难事先知晓处理变量如何动态影响结果。而因果机器学习的本质，仍在于“获取、理解并有效利用正确的数据”。我们必须确保所选的干扰变量构成一个充分的调整集，同时建立可靠的评估体系来验证因果模型的有效性。唯有如此，才能真正依托由 FWL 定理启发的理论框架，做出稳健且具有实际意义的因果推断。
参考文献
- 《由机器学习和人工智能驱动的应用因果推理》
arxiv.org/abs/2403.02467- 《勇敢者与真知者的因果推理》
matheusfacure.github.io/python-causality-handbook/landing-page.html
京公网安备 11010802022788号
