AIPW因果推断入门指南：从理论到Python实战

1次阅读

没有评论

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

在数据分析中，AB 测试是评估干预效果的常见方法，但它有一个根本性局限：只能测量相关性，无法直接推断因果关系。尤其是在存在混淆变量（Confounders）的情况下，简单的组间比较会产生严重偏差。

举个例子，假设我们想评估某药物对血压的影响。如果服药组本身年龄偏大，而年龄又会影响血压，那么直接比较两组血压差异就会包含年龄带来的干扰。这就是为什么我们需要专门的方法来控制混淆变量——因果推断技术应运而生。

$$\hat{\tau}{IPW} = \frac{1}{n}\sum\right]$$}^n\left[\frac{T_iY_i}{e(X_i)} – \frac{(1-T_i)Y_i}{1-e(X_i)

其中 $e(X_i)$ 是倾向得分（Propensity Score）。IPW 完全依赖倾向得分模型，如果模型设定错误，估计就会有偏。

$$\hat{\tau}{Reg} = \frac{1}{n}\sum_0(X_i)]$$}^n[\hat{m}_1(X_i) – \hat{m

这种方法完全依赖结果回归模型。同样，如果回归模型错误，估计也会不准。

$$\hat{\tau}{AIPW} = \frac{1}{n}\sum_0(X_i))\right]$$}^n\left[\hat{m}_1(X_i) – \hat{m}_0(X_i) + \frac{T_i}{e(X_i)}(Y_i-\hat{m}_1(X_i)) – \frac{1-T_i}{1-e(X_i)}(Y_i-\hat{m

AIPW 的关键优势在于双重稳健性(Double Robustness)：只要倾向得分模型或结果回归模型有一个正确，估计就是一致的。

from sklearn.datasets import make_regression
import pandas as pd
import numpy as np

# 生成包含混淆变量的数据
def generate_data(n_samples=1000):
    X, _ = make_regression(n_samples=n_samples, n_features=5, noise=0.1)
    # 倾向得分模型
    ps_score = 1 / (1 + np.exp(-X[:, 0] - 0.5*X[:, 1]))
    T = np.random.binomial(1, ps_score)
    # 结果模型
    y = 2*T + X[:, 0] + 0.5*X[:, 1] + np.random.normal(0, 0.5, size=n_samples)
    return pd.DataFrame(X, columns=[f'X{i}' for i in range(5)]), T, y

X, T, y = generate_data()

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split

# 拆分训练 / 测试集检测过拟合
X_train, X_test, T_train, T_test = train_test_split(X, T, test_size=0.3)

# 拟合逻辑回归模型
ps_model = LogisticRegression(penalty='l2', C=1.0)
ps_model.fit(X_train, T_train)

# 检查测试集 AUC
from sklearn.metrics import roc_auc_score
print(f'Test AUC: {roc_auc_score(T_test, ps_model.predict_proba(X_test)[:,1]):.3f}')

from econml.dml import LinearDML
from sklearn.linear_model import LinearRegression

# 初始化模型
dml_estimate = LinearDML(model_y=LinearRegression(),
                        model_t=LogisticRegression(),
                        discrete_treatment=True)

# 拟合数据
dml_estimate.fit(y, T, X=X)

# 输出平均处理效应(ATE)
print(f'ATE: {dml_estimate.ate_:.3f}')

import matplotlib.pyplot as plt

# 获取各样本的处理效应
te = dml_estimate.effect(X)

plt.figure(figsize=(10,5))
plt.hist(te, bins=30, alpha=0.7)
plt.axvline(x=dml_estimate.ate_, color='r', linestyle='--', 
            label=f'ATE={dml_estimate.ate_:.2f}')
plt.xlabel('Individual Treatment Effect')
plt.ylabel('Count')
plt.legend()
plt.show()