两阶段实验统计学推断

作者: 老姚记事本 | 来源:发表于2021-06-25 10:53 被阅读0次

为了选择效果最好的方案，在进行A/B实验中，经常会同时设置多个实验组，并且分别与对照组比较。这样做会有两个问题：1. 多重比较问题；2. 衡量效果时存在选择性偏差。
为了解决这个问题，一种推荐的做法是两阶段实验，先进行一轮多方案比较，将最好的方案再单独进行一轮实验。传统做法会将第二轮的评估作为最终效果，但是会有另一个问题：3. 浪费了第一阶段的数据，降低了结果的敏感度。
微软的实验平台采用了一种将两阶段结果合并的方案——Statistical Inference in Two-Stage Online Controlled Experiments with Treatment Selection and Validation。

两阶段实验的背景

1. 多个实验方案的影响

在第一阶段中，多方案对实验评估影响如下：

实验效果的影响
单实验方案时， $\Delta$ 是实验效果的无偏最大似然估计；
多实验方案时， $max(\Delta^{(i)})$ 不再是无偏的；
多重比较影响
会增加发现“显著”方案的概率，可参考我之前写过的——A/B实验设计——如何避免多重检验错误.

2. 弱相关性

先不考虑多方案，综合评估两阶段实验的一个挑战是两次实验具有相关性。
为了规避这个问题，可以将两次实验的流量互斥，但是有两个缺陷：1.每阶段可用流量少了，合计不能超过100%；2.需要额外的工程方面成本。
对线上随机实验来说，两阶段的 $\Delta$ 可以近似认为是独立的，以下是解释：

假设一个用户参与了两阶段的实验，当他在对照组时，观测结果为 $(X_i, Y_i)$ ：
$X_i = \mu^{(1)} + \alpha_i + \epsilon_i,$
$Y_i = \mu^{(2)} + \alpha_i + \zeta_i,$
其中 $\mu^{(1)}, \mu^{(2)}$ 分别代表两阶段的效果期望值， $a_i$ 代表个人的随机影响， $\epsilon_i,\zeta_i$ 是对应的随机噪音。
则对照组的用户可表示为：
$X_i = \mu^{(1)} + \theta_i + \alpha_i + \epsilon_i,$
$Y_i = \mu^{(2)} + \theta_i +\alpha_i + \zeta_i,$
其中 $\theta_i$ 是实验干预影响。

设N为总线上用户数量，m为第一阶段实验组用户数、n为对照组，第二阶段对应用户数为m'、n'。如果两次用户分组都是随机的，而两阶段的效果分别为 $\Delta_1,\Delta_2$ ，则：

定理1 (Almost Uncorrelated Deltas)
假如随机分组，则：
$Cov(\Delta_1,\Delta_2) = Var(\theta)/N$ (1)
另外，如果 $Var(\theta) \leq \rho Var(X)$ ，而且 $Var(\theta) \leq \rho Var(Y)$ ，则：

$Corr(\Delta_1, \Delta_2) \leq \rho$

如何理解定理1？
$Var(\Delta_1) = Var(X^{(t)})/m + Var(X^{(t)})/n$
$Var(\Delta_2) = Var(Y^{(t)})/m' + Var(Y^{(t)})/n'$
我们感兴趣的是相关性：
$Corr(\Delta_1, \Delta_2) = \frac{Cov(\Delta_1, \Delta_2)}{\sqrt{Var(\Delta_1)\times Var(\Delta_2)} } = \frac{Var(\theta)/N}{\sqrt{Var(\Delta_1)\times Var(\Delta_2)} }$ (2)
如果治疗效果与X有比例关系: $\theta / X = d\%$ ，则 $Var(\theta)= (d\%)^2Var(X)$ ，根据 (1)、(2)可推导出：
$Corr(\Delta_1, \Delta_2) \leq (d\%)^2$ 。

由于在线上随机实验中，往往提升都是比较小的，即使按提升10%来看，则两次实验的 $\Delta$ 相关系数会小等于0.01。
通过上述定理可知，两阶段实验的 $\Delta$ 相关性很弱，可以近似认为不相关。

解决方案

检验方面

需要可以将两阶段合并进行检验的方案。

1. 传统方法

多重比较
比如通过Bonferroni correction矫正，可参考A/B实验设计——如何避免多重检验错误
计算统计显著性
荟萃分析，比如通过Fisher’s method，可参考荟萃分析简介

通过这两步可以有效的控制假阳性，但是条件比较严格，会降低敏感度，也就是假阴性较高。

2. Sharp Multiple Comparison Adjustment

传统做法偏保守，其实在本场景下可以有跟敏感的方案。
设第一阶段实验中，对照组指标数值为 $X_0$ , 实验组为 $X_1,...,X_k$ ，可得到 $\Delta_1,...,\Delta_k$

Generalized step-down procedure
可知 $\Delta_1,...,\Delta_k$ 服从一个多元正态分布，而我们要的统计量是 $max(\Delta_i)$
基于零假设下，如果通过蒙特卡洛模拟大量次，每次计算 $max(\Delta_i)$ ，这样就可以得到 $max(\Delta_i)$ 的经验分布，可以带入实验中真实值得到P值。
Generalized weighted average test
结合到我们的场景，如何计算两阶段的合并效果？
根据定理1，利用combined Z-score思路：
$\Delta = wmax(\Delta_i) + (1 - w)\Delta_*$
通过蒙特卡洛得到 $max(\Delta_i)$ 的经验分布，与模拟验证阶段的 $\Delta_*$ 结合，可以得到 $\Delta$ 的经验分布。
于是可以精准得到两阶段的合并效果，很符合预设假阳性水平。

点估计方面

A/B 另一个挑战是要给出效果估计。
在第一个筛选阶段，如果我们有k个实验组，则 $X = (X_1, ... X_k) \sim N(\mu, \Sigma)$ ；对照组为 $X_0 \sim N(\mu_0, \sigma_0)$ 。
不失一般性的，我们假设 $\sigma_0 = \sigma_1 = ...= \sigma_k$ ；另外 $\Delta_i = X_i - X_0$ ，参数 $maxi = argmax_i\Delta_i = argmax_iX_i$ 。
在验证阶段， $X_{maxi}^* \sim N(\mu_{maxi}, \sigma'^2)$ ，可以得出 $\Delta_{maxi}^*$ ，它是真实干预效果 $\theta = \mu_{maxi} - \mu_0$ 的无偏估计，但是 $\Delta_{maxi}$ 是有偏的。如何与 $\Delta_{maxi}$ 结合，并且纠正偏斜呢？
定义bias为 $\lambda(\mu) = E(X_{maxi} - \mu_{maxi})$ ，需要找出它的估计方法。

Naive-correction estimator
将筛选阶段得到的各组效果 $x$ 作为 $\mu$ 的点估计，通过 $N(x, I\sigma^2)$ 来蒙特卡洛得到 $\lambda$ 的估计，则可以用 $\hat{\mu}_{maxi,naive} = X_{maxi} - \hat{\lambda}$ 作为其点估计。
需要注意的是，这是一种“plug-in"估计，因为 $\hat{\lambda} = \lambda(x)$ 。
将各实验组结果排序 $X_1,X_2,X_2...$ ，定义 $H(X) = \frac{X_{(k)} - X_{(k - 1)}}{\sqrt{2\sigma^2}}$ ，通过模拟得到 $H$ 与 $\lambda$ 关系如下：
The λ(µ) = E(xmaxi − µmaxi) and E(H(x)) of 2000 randomly sampled µ’s.
因此可以通过单参数模型拟合，通过 $H$ 计算 $\lambda$
Simple linear-correction estimator
沿袭上一节，令 $y_b \sim N(x, I\sigma^2), b = 1...B$ ，定义 $d(y_b) = max(y_b) - x_{argmaxi(y_b)}$ ，可以通过多项式线性模型拟合： $d(y) = f(H(y))$
这里原文推荐使用natural cubic splines模型，可参考三次样条（cubic spline）插值。
之后可以用它的 $\hat{\lambda}$ 矫正bias，这种做法被微软称为NS方法。
Bayesian posterior driven linear-correction estimator
上面通过MLE矫正bias的做法，会有较大方差，可以结合贝叶斯进行优化。
假设有先验 $µ \sim N (0, \tau^2I)$ ，通过James-Stein estimator得到后验：
$\mu|X \sim N((1 - \frac{(k-1)\sigma^2}{\left \| X \right \|^2_2})X,(1 - \frac{(k-1)\sigma^2}{\left \| X \right \|^2_2})I\sigma^2)$
采用分层结构，先通过上面分布抽样得到 $\mu$ ，在通过 $\mu$ 分布抽样得到 $y_b$ ，后续流程与之前一直。
此方法称为JS-NS估计法，它会有更小的方差。
Beyond simple linear-correction
可以采用更复杂的模型来取代Simple linear-correction，但是微软的研究发现没有明显区别；将蒙特卡洛阶段次数增加到10000次后，再增加模拟次数没有明显区别。

方案验证

以下是微软的数据验证结论。

1. 假设检验

Type I error under the null hypothesis

固定 $\theta$ 下，模拟10000次：
BF: 3.2%假阳性，符合之前设想，过于严格；
WAvg: 5.1%假阳性，接近与期望水平。
$\theta \sim N(0, 0.04)$ ，模拟10000次：
BF: 3.3%；
WAvg: 5.0%.
$\theta \sim N(0, 1)$ ，模拟10000次：
BF: 4.2%；
WAvg: 12.3%。

结论：当效果较固定时(上文假设)，两者都可以较好控制假阳性；当效果较随机时，BF方法更加健壮。

Statistical power under the alternative hypothesis

Power curve of generalized weigthed average method, BF method and validation run t-test.

如图，采用WAvg方法，power会明显优于单独验证二次实验或者BF方法。

2. 点估计

The estimated MSE for candidate estimators on 16 of the 200 randomly generated µ’s. The µ’s are ordered according the value H(µ) and are grouped by H(µ) < 1 and H(µ) > 1.
Naive-Correction在

H(\mu)

较大时，会有较小的bias和较大的方差；
JS-NS在

H(\mu)

较小时，有最好的效果，而且方差始终小于其他方法。因此JS-NS是最好的选择。

后记

现在的实验系统大部分不支持分阶段实验的荟萃分析，但是在样本量受限等场景下，此类方法是非常有用的，但是也会更加复杂。
需要注意，文中没有考虑实验的残留效应，实际使用中需要考虑。