---
title: "经济与商务实证研究方法 - 第 5 周:因果机器学习"
subtitle: "完整讲义:高维控制、LASSO、DML、因果森林、异质性效应与人工智能辅助审计"
author: "陈志远"
institute: "中国人民大学商学院"
date: "2026-05-22"
format:
html:
theme: cosmo
css: ../../lecture-notes.css
html-math-method: mathml
toc: true
toc-depth: 3
number-sections: true
code-fold: false
code-tools: true
highlight-style: github
self-contained: true
embed-resources: true
page-layout: article
execute:
echo: true
warning: false
message: false
eval: true
cache: false
fig-width: 10
fig-height: 6
dpi: 150
lang: zh
jupyter: rmeb-venv
---
```{python}
#| echo: false
#| output: false
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
plt.rcParams["font.sans-serif"] = ["Arial Unicode MS", "SimHei", "DejaVu Sans"]
plt.rcParams["axes.unicode_minus"] = False
plt.rcParams["figure.dpi"] = 150
np.random.seed(42)
```
# 引言 {#sec-intro}
机器学习(Machine Learning, ML)进入经济学和管理学实证研究后,最容易出现的误解是:只要模型预测准确,就可以得到更好的因果推断。这个判断只对了一半。机器学习确实能帮助我们处理高维控制、复杂非线性和异质性效应,但因果推断的核心仍然是识别。也就是说,机器学习负责更好地估计条件期望、倾向得分或局部异质性;研究者仍然要说明为什么处理变量中剩下的变化可以被解释为因果变化。
本讲围绕三个问题展开。第一,预测和因果推断的目标函数为什么不同。第二,最小绝对收缩与选择算子(Least Absolute Shrinkage and Selection Operator, LASSO)、事后 LASSO、双重/去偏机器学习(Double/Debiased Machine Learning, DML)如何帮助我们在高维控制下估计平均处理效应。第三,因果森林如何把问题从平均处理效应推进到异质性处理效应,并如何把机器学习发现的异质性桥接到结构估计中的偏好、成本或技术原语。
# 第一部分:预测与因果推断 {#sec-prediction-causality}
## 两种目标函数
机器学习的标准目标是预测。给定特征 $X$,学习一个函数 $\hat{f}(X)$,让测试集预测误差尽可能小:
$$
\min_{\hat{f}} E[(Y-\hat{f}(X))^2]
$$
或者在分类问题中最大化准确率、AUC 等指标。训练集 / 测试集划分、交叉验证、正则化、套袋法、随机森林、提升法,都是围绕样本外预测组织起来的。
在课堂讨论中,可以把预测问题概括为“尽量猜准 $Y$”。只要某个变量能提高样本外预测表现,预测模型通常愿意使用它。但因果推断不能这样处理变量:我们还必须知道变量是在处理前还是处理后出现,是否是中介变量,是否是碰撞变量,是否会改变我们想估计的总效应。
因果推断关心的是潜在结果:
$$
\tau = E[Y(1)-Y(0)]
$$
或条件处理效应:
$$
\tau(x)=E[Y(1)-Y(0)\mid X=x]
$$
预测 $Y$ 的模型可以使用任何与 $Y$ 相关的变量;因果模型却必须区分变量在因果图中的位置。处理前混杂变量通常需要控制,处理后中介变量可能不能控制,碰撞变量控制后会引入偏误,工具变量也不能简单作为普通控制处理。
## 机器学习的合适角色
机器学习在因果推断中的合适角色是估计干扰函数。例如在部分线性模型中:
$$
Y = \theta_0 D + g_0(X) + U
$$
研究者关心 $\theta_0$,而 $g_0(X)$ 是一个可能复杂的控制函数。若 $X$ 维度很高、存在非线性和交互,用随机森林、提升法、LASSO 或神经网络估计 $g_0$ 可能比线性模型更好。同样,处理方程
$$
D = m_0(X) + V
$$
中的 $m_0(X)$ 也可以用机器学习估计。DML 的思想就是让这些预测任务服务于因果参数估计。这里要特别强调:机器学习改善的是“如何控制”,不是自动解决“为什么可识别”。
# 第二部分:LASSO 与事后 LASSO {#sec-lasso}
## 高维控制与近似稀疏
当控制变量数量 $p$ 接近甚至超过样本量 $n$,普通最小二乘法(Ordinary Least Squares, OLS)会变得不稳定。很多经验研究还会主动构造大量技术变量,例如多项式、交互项、固定效应组合、文本特征、滞后项和分组指标。此时我们需要正则化。
LASSO 的目标函数为
$$
\hat{\beta}_{LASSO}
=
\arg\min_b
\sum_{i=1}^n (Y_i-X_i'b)^2
+
\lambda\sum_{j=1}^p |b_j|\hat{\psi}_j
$$
$\ell_1$ 惩罚会把一部分系数设为 0,因此 LASSO 同时完成收缩和变量选择。LASSO 最适合近似稀疏情形:真实系数未必严格为 0,但重要系数快速衰减,用少数变量即可近似条件期望。它的代价是收缩偏差,收益是高维情况下更稳定的预测和更可控的模型复杂度。
## 事后 LASSO 与变量选择后的推断
LASSO 的收缩会带来系数偏差。事后 LASSO 先用 LASSO 选择变量集合 $\hat{S}$,再在 $\hat{S}$ 上做 OLS:
$$
\hat{\beta}_{post}
=
\arg\min_{\beta_j=0, j\notin \hat{S}}
\sum_{i=1}^n (Y_i-X_i'\beta)^2
$$
事后 LASSO 往往改善预测表现,因为它去掉了 LASSO 惩罚导致的收缩偏差。但在因果推断中,变量选择后的 OLS 不能直接当作普通 OLS。标准误和置信区间没有自动考虑“先选择、再估计”的不确定性。
一个关键问题是:只选择预测 $Y$ 的变量,可能漏掉预测 $D$ 的混杂变量。双重选择的做法是分别用 LASSO 预测 $Y$ 和 $D$,把两个变量集合取并集,再估计处理效应。这个方法的目标不是得到最简模型,而是降低遗漏关键混杂变量的概率。对于课堂讲解,可以用一句话概括:只要一个变量可能影响处理分配,就不能因为它对结果预测力弱而轻易丢掉。
# 第三部分:双重/去偏机器学习 {#sec-dml}
## 残差化直觉
DML 最常见的模型是部分线性模型:
$$
Y = \theta_0 D + g_0(X) + U,\quad E[U\mid X,D]=0
$$
$$
D = m_0(X) + V,\quad E[V\mid X]=0
$$
算法的直觉是 Frisch-Waugh-Lovell 定理的机器学习版本。先用机器学习估计
$$
\hat{g}(X)\approx E[Y\mid X]
$$
和
$$
\hat{m}(X)\approx E[D\mid X]
$$
然后构造残差:
$$
\tilde{Y}=Y-\hat{g}(X)
$$
$$
\tilde{D}=D-\hat{m}(X)
$$
最后回归
$$
\tilde{Y}=\theta_0\tilde{D}+\text{error}
$$
处理效应来自 $D$ 中不能被 $X$ 预测的部分,而不是原始的 $D$。这一步就是去偏的直觉。课堂上可以把它说成“先净化,再比较”:先把 $Y$ 和 $D$ 中都能由控制变量解释的部分拿掉,再看剩下的共同变化。
下面用模拟数据展示残差化后的斜率。
::: {style="font-size: 1em;"}
```{python}
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
n = 600
X = np.random.normal(size=(n, 5))
D = 0.8 * X[:, 0] - 0.4 * X[:, 1] + np.random.normal(size=n)
Y = 2.0 * D + 1.5 * np.sin(X[:, 0]) + X[:, 2] ** 2 + np.random.normal(size=n)
g_hat = RandomForestRegressor(
n_estimators=80, min_samples_leaf=10, random_state=1
).fit(X, Y).predict(X)
m_hat = RandomForestRegressor(
n_estimators=80, min_samples_leaf=10, random_state=2
).fit(X, D).predict(X)
Y_res = Y - g_hat
D_res = D - m_hat
theta_res = LinearRegression().fit(D_res.reshape(-1, 1), Y_res).coef_[0]
theta_res
```
:::
```{python}
#| echo: false
fig, ax = plt.subplots(figsize=(8, 4))
ax.scatter(D_res, Y_res, alpha=0.35, s=18, color="#2E6E9E")
grid = np.linspace(D_res.min(), D_res.max(), 100)
ax.plot(grid, theta_res * grid, color="#B45F06", linewidth=2)
ax.set_xlabel("处理残差:D - E[D|X]")
ax.set_ylabel("结果残差:Y - E[Y|X]")
ax.set_title(f"残差化后的因果斜率约为 {theta_res:.2f}")
plt.show()
```
## 正交性与交叉拟合
DML 的理论核心是 Neyman 正交性。一个典型正交分数为
$$
\psi(W;\theta,\eta)
=
\left(Y-g(X)-\theta[D-m(X)]\right)
\left(D-m(X)\right)
$$
其中 $\eta=(g,m)$。正交性要求干扰函数的小误差不会一阶影响 $\theta$:
$$
\partial_{\eta}E[\psi(W;\theta_0,\eta)]\big|_{\eta=\eta_0}=0
$$
这让我们可以使用相对灵活的机器学习方法估计 $g$ 和 $m$。但如果同一批数据既训练机器学习模型又估计因果参数,过拟合仍会污染残差。交叉拟合的做法是把样本分成 $K$ 折,每次用 $K-1$ 折训练干扰模型,在留出折上生成折外残差,最后汇总所有残差估计 $\theta$。
# 第四部分:因果森林与异质性效应 {#sec-causal-forest}
## 从树模型到因果森林
回归树把特征空间切成多个区域,在每个叶节点用局部均值预测 $Y$。随机森林通过自助抽样和随机特征子集构造许多树,再平均预测:
$$
\hat{f}_{RF}(x)=\frac{1}{B}\sum_{b=1}^B \hat{f}^{(b)}(x)
$$
因果森林的目标不是预测 $Y$,而是估计
$$
\tau(x)=E[Y(1)-Y(0)\mid X=x]
$$
它会寻找处理效应差异大的区域,并在局部比较处理组与控制组。现代因果森林通常结合倾向得分调整、正交化和诚实分裂。
## 诚实推断
诚实推断中的“诚实”含义是:一部分样本用于决定树如何分裂,另一部分样本用于估计叶节点中的处理效应。这样避免同一数据既“寻找异质性”又“估计异质性”,降低数据挖掘后的过度乐观。一个适合课堂的类比是:先用一部分样本“出题”,再用另一部分样本“阅卷”,避免模型对自己刚刚发现的模式过度自信。
解释异质性处理效应时必须克制。变量重要性和异质性曲线是诊断工具,不是机制解释本身。研究者应检查异质性是否稳定、是否由处理前变量定义、分组样本是否足够、是否存在选择性展示和多重检验。换句话说,异质性处理效应可以告诉我们“哪里不同”,但不能自动告诉我们“为什么不同”。
# 第五部分:课堂实验:用 EconML 实现 DML {#sec-lab}
本实验使用 scikit-learn 内置糖尿病数据集。结果变量是疾病进展指标。为了演示 DML,我们把身体质量指数是否高于中位数定义为二元处理变量,其他健康指标作为控制变量。这个例子用于教学,不能自动解释为身体质量指数对疾病进展的因果效应;因果解释仍依赖条件独立、共同支撑和没有坏控制等假设。
## 数据准备
::: {style="font-size: 1em;"}
```{python}
from sklearn.datasets import load_diabetes
data = load_diabetes(as_frame=True)
df = data.frame.copy()
Y = df["target"].to_numpy()
D = (df["bmi"] > df["bmi"].median()).astype(int).to_numpy()
X = df.drop(columns=["target", "bmi"])
pd.DataFrame({
"n_obs": [len(df)],
"n_controls": [X.shape[1]],
"treated_share": [D.mean()]
}).round(3)
```
:::
## 朴素 OLS 与控制变量 OLS
::: {style="font-size: 1em;"}
```{python}
import statsmodels.api as sm
ols_naive = sm.OLS(Y, sm.add_constant(D)).fit(cov_type="HC3")
X_ols = pd.concat(
[pd.Series(D, name="D").reset_index(drop=True), X.reset_index(drop=True)],
axis=1
)
ols_controls = sm.OLS(Y, sm.add_constant(X_ols)).fit(cov_type="HC3")
pd.DataFrame({
"model": ["朴素 OLS", "控制变量 OLS"],
"effect": [ols_naive.params[1], ols_controls.params["D"]],
"se": [ols_naive.bse[1], ols_controls.bse["D"]]
}).round(3)
```
:::
## 使用 EconML 估计 DML
课程环境已经配置 `econml`。下面先给出一个只依赖 `sklearn` 的手写交叉拟合 DML 版本,随后用 `EconML` 的 `LinearDML` 复现同一思路。两段代码的目的不是比较软件包优劣,而是让学生看清楚 DML 的核心步骤:训练干扰模型、生成折外残差、用残差估计处理效应。
::: {style="font-size: 1em;"}
```{python}
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
from sklearn.model_selection import KFold
from sklearn.base import clone
from sklearn.linear_model import LinearRegression
y_res = np.zeros_like(Y, dtype=float)
d_res = np.zeros_like(D, dtype=float)
kf = KFold(n_splits=5, shuffle=True, random_state=42)
model_y = RandomForestRegressor(
n_estimators=60,
min_samples_leaf=10,
random_state=42
)
model_d = RandomForestClassifier(
n_estimators=60,
min_samples_leaf=10,
random_state=42
)
for train_idx, test_idx in kf.split(X):
gy = clone(model_y)
md = clone(model_d)
gy.fit(X.iloc[train_idx], Y[train_idx])
md.fit(X.iloc[train_idx], D[train_idx])
y_res[test_idx] = Y[test_idx] - gy.predict(X.iloc[test_idx])
d_res[test_idx] = D[test_idx] - md.predict_proba(X.iloc[test_idx])[:, 1]
manual_dml = sm.OLS(y_res, sm.add_constant(d_res)).fit(cov_type="HC3")
pd.DataFrame({
"model": ["手写 DML"],
"effect": [manual_dml.params[1]],
"se": [manual_dml.bse[1]]
}).round(3)
```
:::
::: {style="font-size: 1em;"}
```{python}
from econml.dml import LinearDML
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
dml = LinearDML(
model_y=RandomForestRegressor(
n_estimators=60,
min_samples_leaf=10,
random_state=42
),
model_t=RandomForestClassifier(
n_estimators=60,
min_samples_leaf=10,
random_state=42
),
discrete_treatment=True,
cv=5,
random_state=42
)
dml.fit(Y, D, X=X)
ate = dml.ate(X=X)
lb, ub = dml.ate_interval(X=X)
pd.DataFrame({
"model": ["LinearDML"],
"effect": [ate],
"ci_low": [lb],
"ci_high": [ub]
}).round(3)
```
:::
## 解释去偏步骤
朴素 OLS 比较高身体质量指数组与低身体质量指数组的平均差异。控制变量 OLS 加入线性控制变量。DML 则用随机森林分别预测 $Y$ 和 $D$,并通过交叉拟合在留出折上构造残差。最终估计的不是原始 $D$ 对原始 $Y$ 的关系,而是
$$
Y-\hat{E}[Y\mid X]
$$
对
$$
D-\hat{E}[D\mid X]
$$
的关系。这就是去偏的核心。需要提醒学生的是,如果身体质量指数高低本身不是在控制变量条件下近似随机的,DML 也不能自动把这个教学例子变成真正的因果证据。
# 第六部分:人工智能辅助特征工程与审计 {#sec-ai}
大语言模型(Large Language Model, LLM)可以帮助生成候选控制变量、非线性项、交互项和 scikit-learn 管线。例如,在企业数据中,可以让大语言模型根据制度背景列出规模、年龄、行业、地区、历史绩效、融资约束、治理结构等候选混杂变量,并生成变量字典。
但人工智能生成变量之后,必须进行因果审计。每个变量都应标注为处理前混杂变量、处理后中介变量、碰撞变量、工具变量、结果代理或无关变量。预测有用不等于因果可控。审计记录至少包括:变量来源、生成规则、是否由人工智能提议、是否进入最终模型、删除原因、随机种子、分折划分和调参网格。
# 第七部分:桥接结构估计 {#sec-structural}
因果机器学习估计的异质性效应可以为结构模型提供线索。例如因果森林发现高收入消费者对补贴更敏感,这提示结构需求模型中可能需要收入异质的价格敏感度。DML 发现政策对高生产率企业效果更大,这可能对应生产函数或投资模型中的技术原语差异。
但机器学习估计的 $\hat{\tau}(x)$ 仍是局部因果效应描述。结构估计要进一步说明异质性来自偏好、成本、技术、约束还是信息,并用这些原语模拟新的政策情景。换句话说,机器学习告诉我们“哪里不同”,结构模型解释“为什么不同”以及“换一个政策会怎样”。
# 小结 {#sec-summary}
本讲的核心是把机器学习放回因果推断框架。预测能力本身不保证因果识别;LASSO 和事后 LASSO 帮助处理高维控制;DML 通过正交化和交叉拟合降低高维预测误差对因果参数的影响;因果森林把平均效应推进到异质性效应。人工智能可以加快特征工程和变量选择,但必须留下可审计的管线。下一讲进入结构估计,把异质性、选择和反事实放入更明确的经济决策系统。