因果推断实战：基于Propensity得分预测实验组与对照组的效果差异

1次阅读

共计 4302 个字符，预计需要花费 11 分钟才能阅读完成。

在观察性研究或非随机实验中，直接比较实验组（Treatment Group）和对照组（Control Group）的结果会产生偏差，因为两组样本的协变量分布可能不平衡。Propensity 得分（倾向性评分）通过估计个体接受处理的概率来解决这一问题，其定义为：

$$e(X) = P(T=1|X)$$

其中 $T$ 为处理变量，$X$ 为协变量。匹配后理论上满足条件独立假设 $Y(1), Y(0) \perp T | e(X)$，但实际仍存在两个关键问题：

选择性偏差残余 ：即使匹配后，未观测变量仍可能导致偏差
模型误设风险 ：Propensity 得分模型或结果模型的错误设定会传递至效应估计

T-Learner
分别训练实验组模型 $\hat{\mu}_1(x)$ 和对照组模型 $\hat{\mu}_0(x)$
处理效应：$\hat{\tau}(x) = \hat{\mu}_1(x) – \hat{\mu}_0(x)$
复杂度：$O(2n)$，适合处理组样本量均衡场景
X-Learner
阶段一：同 T -Learner 构建基模型
阶段二：对实验组样本估算反事实 $\hat{D}_1 = Y_1 – \hat{\mu}_0(X_1)$，对照组样本估算 $\hat{D}_0 = \hat{\mu}_1(X_0) – Y_0$
阶段三：训练两个效应模型 $\hat{\tau}_1(x), \hat{\tau}_0(x)$ 并加权组合
复杂度：$O(3n)$，对小样本处理组更稳健
R-Learner
通过残差学习直接建模条件平均处理效应（CATE）
目标函数：$\hat{\tau}(\cdot) = \arg\min_{\tau}\left{\frac{1}{n}\sum_{i=1}^n\left[\left(Y_i-\hat{m}^{(-i)}(X_i)\right) – \left(T_i-\hat{e}^{(-i)}(X_i)\right)\tau(X_i)\right]^2\right}$
复杂度：$O(n^2)$，适合高维稀疏特征

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split

# 生成模拟数据
np.random.seed(42)
n_samples = 2000
X = np.random.normal(size=(n_samples, 5))
# 真实 propensity 得分
true_ps = 1 / (1 + np.exp(-X[:, 0] - 0.5*X[:, 1]))
T = np.random.binomial(1, true_ps)
# 潜在结果模型
y0 = X[:, 0] + 0.5*X[:, 1] + np.random.normal(scale=0.1)
y1 = y0 + 0.8 + 0.3*X[:, 2]  # 异质性处理效应
Y = np.where(T == 1, y1, y0)

# 估计 propensity 得分
ps_model = LogisticRegression(penalty='l2', C=1.0)
ps_model.fit(X, T)
propensity_scores = ps_model.predict_proba(X)[:, 1]

import matplotlib.pyplot as plt
from sklearn.neighbors import NearestNeighbors

# 最近邻匹配
nn = NearestNeighbors(n_neighbors=1)
matched_pairs = []
for i in np.where(T == 1)[0]:
    # 在对照组中寻找最近邻
    control_mask = (T == 0) & (propensity_scores >= propensity_scores[i] - 0.1) & (propensity_scores <= propensity_scores[i] + 0.1)
    if sum(control_mask) > 0:
        nn.fit(X[control_mask])
        _, match_idx = nn.kneighbors(X[i].reshape(1, -1))
        matched_pairs.append((i, np.where(control_mask)[0][match_idx[0][0]]))

# 可视化匹配前后协变量差异
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
for j in range(X.shape[1]):
    ax1.scatter([j], np.mean(X[T==1, j]) - np.mean(X[T==0, j]), 
                color='red', label='原始' if j == 0 else None)
    matched_X1 = X[[p[0] for p in matched_pairs]]
    matched_X0 = X[[p[1] for p in matched_pairs]]
    ax1.scatter([j], np.mean(matched_X1[:, j]) - np.mean(matched_X0[:, j]), 
                color='blue', label='匹配后' if j == 0 else None)
ax1.axhline(0, linestyle='--', color='gray')
ax1.set_title('协变量均值差异')
ax1.legend()

# 倾向得分分布检查
ax2.hist(propensity_scores[T==0], bins=30, alpha=0.5, label='Control')
ax2.hist(propensity_scores[T==1], bins=30, alpha=0.5, label='Treatment')
ax2.set_title('Propensity Score Distribution')
ax2.legend()
plt.show()

from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_squared_error

# 划分训练集 / 测试集
X_train, X_test, T_train, T_test, Y_train, Y_test = train_test_split(X, T, Y, test_size=0.3, stratify=T)

# 构建子模型
model_t = GradientBoostingRegressor(n_estimators=100, max_depth=3)
model_c = GradientBoostingRegressor(n_estimators=100, max_depth=3)

# 分别在实验组和对照组上训练
model_t.fit(X_train[T_train == 1], Y_train[T_train == 1])
model_c.fit(X_train[T_train == 0], Y_train[T_train == 0])

# 计算 ATE
ate = model_t.predict(X_test).mean() - model_c.predict(X_test).mean()
print(f"Estimated ATE: {ate:.3f} (True ATE: 0.800)")

# 计算 CATE
cate = model_t.predict(X_test) - model_c.predict(X_test)
print(f"CATE range: [{cate.min():.3f}, {cate.max():.3f}]")

import shap

# SHAP 分析处理效应异质性
explainer = shap.Explainer(model_t, X_test)
shap_values_t = explainer(X_test)
explainer = shap.Explainer(model_c, X_test)
shap_values_c = explainer(X_test)

# 计算特征对处理效应的贡献
shap_diff = np.abs(shap_values_t.values - shap_values_c.values)
mean_shap_diff = np.mean(shap_diff, axis=0)

plt.figure(figsize=(10, 5))
plt.bar(range(X.shape[1]), mean_shap_diff)
plt.xticks(range(X.shape[1]), [f"X{i}" for i in range(X.shape[1])])
plt.title("Feature Importance for Heterogeneous Treatment Effect")
plt.show()

必须确保共同支持域（Common Support）：
$$0 < e(X) < 1$$
实际检查方法：
可视化 Propensity 得分分布直方图
计算标准化均值差异（SMD）：
$$SMD = \frac{|\bar{X}_1 – \bar{X}_0|}{\sqrt{(s_1^2 + s_0^2)/2}}$$
SMD > 0.1 表明存在显著不平衡

广义 Propensity 得分模型：
$$r(t,x) = \frac{f_{T|X}(t|x)}{f_T(t)}$$
可通过核密度估计或深度学习模型实现

双重选择 Lasso（Double-Selection LASSO）：
用 Lasso 选择影响 T 的重要特征
用 Lasso 选择影响 Y 的重要特征
合并两组特征后估计处理效应
目标函数：
$$\min_{\beta, \gamma} \frac{1}{2n}|Y – D\beta – X\gamma|_2^2 + \lambda_1|\beta|_1 + \lambda_2|\gamma|_1$$