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

1次阅读

没有评论

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

在因果推断分析中，Propensity Score（倾向性得分）匹配是一种常用的消除混杂变量的方法。它通过估计每个样本接受处理的概率，然后将处理组和对照组中得分相似的样本进行匹配，从而模拟随机对照试验的环境。

然而，Propensity Score 匹配只是第一步，匹配后如何正确预测实验组和对照组的因果效应才是真正关键的问题。许多开发者在匹配后直接比较两组的结果，这可能会引入新的偏差。本文将介绍几种更可靠的方法来估计因果效应。

在完成 Propensity Score 匹配后，我们实际上创建了一个平衡的数据集，其中处理组和对照组在可观测的混杂变量上具有相似分布。这时，我们可以采用以下方法来估计因果效应：

直接比较匹配后的处理组和对照组的平均结果差异
使用回归模型调整剩余的微小差异
采用更高级的双重稳健估计方法

双重稳健估计是一种强大的方法，它结合了倾向性得分模型和结果回归模型的优势。其 ” 双重 ” 特性体现在：

只要倾向性得分模型或结果模型其中一个是正确的，就能得到无偏估计
如果两个模型都正确，估计效率更高

近年来，机器学习方法在因果推断领域取得了显著进展。一些值得关注的模型包括：

Causal Forest：基于随机森林的因果效应估计方法
X-Learner：分步估计处理效应的高级框架
S-Learner/T-Learner：简单但实用的处理效应估计方法

以下是一个完整的 Python 示例，展示如何使用 causalml 库实现上述方法：

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestRegressor
from causalml.propensity import ElasticNetPropensityModel
from causalml.match import NearestNeighborMatch
from causalml.inference.meta import BaseXLearner, BaseRLearner

# 生成模拟数据
np.random.seed(42)
n = 1000
X = np.random.normal(size=(n, 5))
propensity = 1 / (1 + np.exp(-X[:, 0] - 0.5*X[:, 1]))
treatment = np.random.binomial(1, propensity)
y0 = X[:, 0] + X[:, 2] + np.random.normal(size=n)
y1 = y0 + 1 + 0.5*X[:, 1]
y = treatment * y1 + (1-treatment) * y0

# 数据准备
data = pd.DataFrame(X, columns=[f'X{i}' for i in range(X.shape[1])])
data['treatment'] = treatment
data['outcome'] = y

# Propensity Score 匹配
ps_model = ElasticNetPropensityModel()
ps_score = ps_model.fit_predict(data.drop(['treatment', 'outcome'], axis=1), 
                               data['treatment'])
matcher = NearestNeighborMatch(replace=False, ratio=1, random_state=42)
matched_data = matcher.match(data=data, treatment_col='treatment', 
                            score_cols='propensity_score')

# 双重稳健估计
# 步骤 1：建立倾向性得分模型
ps_model = LogisticRegression()
ps_model.fit(matched_data[['X0', 'X1', 'X2', 'X3', 'X4']], matched_data['treatment'])
ps = ps_model.predict_proba(matched_data[['X0', 'X1', 'X2', 'X3', 'X4']])[:, 1]

# 步骤 2：建立结果回归模型
outcome_model = RandomForestRegressor()
outcome_model.fit(matched_data[matched_data['treatment']==0][['X0', 'X1', 'X2', 'X3', 'X4']],
                 matched_data[matched_data['treatment']==0]['outcome'])

# 步骤 3：计算双重稳健估计
def doubly_robust_estimate(df, ps, outcome_model):
    treatment = df['treatment']
    outcome = df['outcome']

    # 预测对照结果
    y0_pred = outcome_model.predict(df[['X0', 'X1', 'X2', 'X3', 'X4']])

    # 计算双重稳健估计
    dr = (treatment * (outcome - y0_pred) / ps + y0_pred).mean() - \
         ((1-treatment) * (outcome - y0_pred) / (1-ps) + y0_pred).mean()

    return dr

dr_effect = doubly_robust_estimate(matched_data, ps, outcome_model)
print(f"双重稳健估计的因果效应: {dr_effect:.3f}")

# 使用 Causal Forest
from causalml.inference.forest import CausalForest

cf = CausalForest(n_estimators=100)
tau = cf.fit_predict(X=matched_data[['X0', 'X1', 'X2', 'X3', 'X4']].values,
                    treatment=matched_data['treatment'].values,
                    y=matched_data['outcome'].values)
print(f"Causal Forest 估计的平均因果效应: {np.mean(tau):.3f}")

在应用上述方法前，务必检查匹配质量：