AB实验与因果推断：原理剖析与工程实践指南

1次阅读

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

在统计学中，辛普森悖论是一个经典的例子，它展示了即使数据在整体层面呈现某种趋势，在子群体中可能完全相反。比如，某个医院可能整体上 A 疗法的治愈率高于 B 疗法，但在每个年龄段的子群体中，B 疗法的治愈率反而更高。这种现象之所以发生，是因为忽略了潜在的混杂变量（如年龄）。

这引出了因果推断的核心问题：我们如何确定一个变量（如疗法）对另一个变量（如治愈率）的因果影响？随机对照实验（RCT）被认为是因果推断的黄金标准，因为它通过随机分配处理组和对照组，确保了所有潜在混杂变量的分布在不同组间是平衡的。

AB 实验本质上是一种随机对照实验，它依赖于以下三大核心假设：

可忽略性（Ignorability）：处理分配与潜在结果独立，即 $Y(1), Y(0) \perp T$。在 AB 实验中，这通过随机化实现。
正值性（Positivity）：每个单元都有机会被分配到处理组或对照组，即 $0 < P(T=1|X) < 1$。
一致性（Consistency）：观察到的结果等于潜在结果，即 $Y = T \cdot Y(1) + (1-T) \cdot Y(0)$。

尽管 AB 实验在理论上很完美，但在实际工程中会遇到许多问题：

非随机分流 ：由于技术或业务限制，分流可能不完全随机，导致选择偏差。
样本污染 ：用户可能在实验期间切换到另一组，破坏随机性。
新奇效应 ：用户对新产品或功能的初始反应可能不反映长期行为。

双重稳健估计（Doubly Robust Estimation）是一种结合了倾向得分和结果模型的估计方法，即使其中一个模型错误，仍能得到一致的估计。其数学形式为：

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

以下是一个基于倾向得分的因果效应估计示例，对比了 t -test 和贝叶斯方法：

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from scipy import stats
import pymc3 as pm

# 生成模拟数据
np.random.seed(42)
n = 1000
X = np.random.normal(0, 1, (n, 2))
T = np.random.binomial(1, 0.5, n)
Y = 0.5 * T + 0.3 * X[:, 0] + 0.2 * X[:, 1] + np.random.normal(0, 0.1, n)

# 倾向得分估计
ps_model = LogisticRegression().fit(X, T)
ps = ps_model.predict_proba(X)[:, 1]

# 双重稳健估计
def doubly_robust(X, T, Y, ps):
    mu1 = np.mean(Y[T == 1])
    mu0 = np.mean(Y[T == 0])
    return mu1 - mu0

effect_dr = doubly_robust(X, T, Y, ps)
print(f"Doubly Robust Estimate: {effect_dr:.4f}")

# t-test
t_stat, p_val = stats.ttest_ind(Y[T == 1], Y[T == 0])
print(f"t-test p-value: {p_val:.4f}")

# 贝叶斯方法
with pm.Model() as model:
    alpha = pm.Normal('alpha', mu=0, sd=1)
    beta = pm.Normal('beta', mu=0, sd=1, shape=2)
    tau = pm.Normal('tau', mu=0, sd=1)

    mu = alpha + pm.math.dot(X, beta) + tau * T
    sigma = pm.HalfNormal('sigma', sd=1)

    y_obs = pm.Normal('y_obs', mu=mu, sd=sigma, observed=Y)

    trace = pm.sample(1000, tune=1000)

print(pm.summary(trace))