因果推断实战：基于Propensity Score匹配后如何预测实验组与对照组效果

1次阅读

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

在因果推断中，我们经常遇到观察性数据（非随机实验数据），这时处理组（实验组）和对照组之间可能存在混杂变量（confounders）。Propensity Score（倾向性得分）匹配通过模拟随机实验环境，帮助我们减少选择偏差。

但匹配只是第一步，真正的挑战在于匹配后：

匹配后的样本量通常会减少
协变量平衡性需要验证
因果效应的估计方法直接影响结论可靠性

匹配后直接比较两组的结果差异：

# 匹配后的 DataFrame 包含匹配对信息
matched_df = matching_algorithm(pscores, treatment_indicator)

# 简单均值比较（不建议单独使用）ate_naive = matched_df[matched_df['treatment']==1]['outcome'].mean() - \
             matched_df[matched_df['treatment']==0]['outcome'].mean()

问题：未考虑残余偏差和测量误差

结合倾向得分和结果模型的双重保护：
1. 先建立结果预测模型（如线性回归）
2. 再结合倾向得分进行修正

from causalml.inference.meta import LRSRegressor

dr = LRSRegressor()
# 输入特征、处理标识、结果变量
ate = dr.estimate_ate(X, treatment, y)

优势：只要倾向得分模型或结果模型有一个正确，估计就是一致的

from econml.grf import CausalForest

cf = CausalForest(n_estimators=100)
cf.fit(X, treatment, y)
# 获取个体处理效应 (ITE)
ite = cf.predict(X)
# 平均处理效应 (ATE)
ate = ite.mean()

S-Learner：单一模型
T-Learner：双模型
X-Learner：交叉纠正

# 使用 causalml 实现
from causalml.inference.meta import XGBoostRegressor

learner = XGBoostRegressor()
ate = learner.estimate_ate(X, treatment, y)

# 环境准备
!pip install causalml econml
import pandas as pd
import numpy as np
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.model_selection import train_test_split

# 生成模拟数据
np.random.seed(42)
n = 2000
X = np.random.normal(size=(n, 5))
# 生成倾向得分
true_pscore = 1/(1+np.exp(-X[:,0]-X[:,1]))
treatment = np.random.binomial(1, true_pscore)
# 生成结果变量
y = treatment * (X[:,0]+X[:,1]) + (1-treatment)*X[:,1] + np.random.normal(size=n)

# 步骤 1：Propensity Score 匹配（使用最邻近匹配）from causalml.match import NearestNeighborMatch

matcher = NearestNeighborMatch(replace=False)
matched_df = matcher.match(data=pd.DataFrame(X), 
                          treatment=treatment,
                          score=pscore)

# 步骤 2：双重稳健估计
from causalml.inference.dr import DRLearner

dr = DRLearner()
# 拟合模型
dr.fit(X, treatment, y)
# 获取 ATE 和置信区间
ate, ate_lb, ate_ub = dr.estimate_ate(interval=True)
print(f"ATE: {ate:.2f}, 95% CI: [{ate_lb:.2f}, {ate_ub:.2f}]")

# 步骤 3：Causal Forest
from econml.grf import CausalForest

cf = CausalForest(n_estimators=500)
cf.fit(X, treatment, y)
# 可视化变量重要性
import matplotlib.pyplot as plt
plt.bar(range(X.shape[1]), cf.feature_importances_)
plt.title('Feature Importance')
plt.show()

标准化均值差 (SMD)：应 <0.1
QQ 图 ：检查变量分布平衡性
t 检验 ：连续变量组间差异

from causalml.match import create_table_one

# 生成平衡性报告
table1 = create_table_one(data=matched_df, 
                        treatment_col='treatment',
                        features=['feature_0','feature_1'])
print(table1)