因果推断实战：基于BART模型的高效因果效应估计与避坑指南

1次阅读

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

在观察性研究中，我们经常遇到非线性关系和隐藏的混杂变量。传统方法如线性回归假设处理效应是恒定的，而倾向得分匹配在高维数据中会面临维度灾难问题。举个例子，当用户行为受数十个交叉特征影响时，逻辑回归可能连倾向得分都估计不准。

和其他现代因果推断方法相比，BART 有三个显著优势：

自动特征交互：通过树结构天然捕捉变量间的复杂交互
内置不确定性：贝叶斯框架直接给出效应估计的置信区间
稳健拟合：对异常值和模型误设相对不敏感

这个表格对比了几种主流方法的表现：

方法	非线性处理	高维协变量	计算效率
线性回归	×	×	✓
因果森林	✓	✓	×
TMLE	✓	✓	×
BART	✓	✓	✓

先安装必要库：

pip install pymc3 bartpy sklearn arviz

我们构造包含隐藏混杂变量的数据集：

import numpy as np

# 模拟 10 个观测变量和 1 个隐藏变量
np.random.seed(42)
n = 2000
X = np.random.normal(size=(n, 10))
hidden = np.random.binomial(1, 0.3, size=n)  # 未观测的混杂变量

# 处理分配机制（与混杂变量相关）ps = 1 / (1 + np.exp(-X[:,0] - 2*hidden))  # 倾向得分
treatment = np.random.binomial(1, ps, size=n)

# 结果变量（包含非线性效应）y = 3*treatment + 2*np.sin(X[:,1]) + 0.5*X[:,2]**2 + hidden + np.random.normal(size=n)

使用 pymc3 实现：

import pymc3 as pm

with pm.Model() as bart_model:
    # 设置先验：200 棵树，深度不超过 3
    mu = pm.BART("mu", X, y, m=200, 
                response="continuous",
                k=2,  # 控制叶节点参数先验
                alpha=0.95,  # 树深度参数
                beta=2.0)  # 叶节点参数

    # 定义似然
    y_obs = pm.Normal("y_obs", mu=mu, sigma=1, observed=y)

    # 采样
    trace = pm.sample(1000, tune=1000, chains=2)

关键参数说明：
– m: 树的数量（通常 200-500）
– alpha: 控制树深度（0.95 对应较浅的树）
– beta: 叶节点参数先验（越大正则化越强）

检查协变量平衡：

from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import roc_auc_score

# 用 GBDT 检测剩余混杂
gb = GradientBoostingClassifier().fit(X, treatment)
print(f"PS 模型 AUC: {roc_auc_score(treatment, gb.predict_proba(X)[:,1])}")
# 如果 >0.7 说明存在明显混杂

可视化后验分布：

import arviz as az

# 绘制处理效应分布
az.plot_posterior(trace["mu"][:,treatment==1] - trace["mu"][:,treatment==0])

当处理组样本极少时（<5%），可以：
1. 使用分层 BART：对不同子群体训练独立模型
2. 调整先验参数：增大 beta 到 3 - 5 增强正则化

# 使用 numba 加速
from bartpy.sklearnmodel import SklearnModel

model = SklearnModel(n_trees=50, n_jobs=-1)  # 并行化
model.fit(X, y)

通过覆盖率检验评估置信区间质量：

coverage = []
for _ in range(100):
    # 重复模拟数据...
    # 计算 95%CI 是否包含真实效应 3
    coverage.append(3 in ci_interval)
print(f"覆盖率: {np.mean(coverage)}")  # 目标≈0.95