原文：
towardsdatascience.com/clustered-standard-errors-in-ab-tests-a993f29b9225
因果数据科学
封面，图片由作者提供
A/B 测试是因果推断的金标准，因为它们允许我们在最小假设下做出有效的因果陈述，这要归功于随机化。事实上，通过随机分配处理（药物、广告、产品等），我们能够比较结果（疾病、公司收入、客户满意度等）在不同主体（患者、用户、客户等）之间的差异，并将结果平均差异归因于处理的因果效应。
有时会出现处理分配的单位与观察的单位不同的情况。换句话说，我们不是独立地对每个观察值是否进行处理做出决定，而是以组为单位。例如，我们可能会决定对某个地区的所有客户进行处理，同时在客户层面观察结果，或者对某个品牌的所有商品进行处理，同时在商品层面观察结果。通常这种情况是由于实际限制造成的。在第一个例子中，所谓的地理实验，这种情况发生是因为我们无法追踪用户，因为 cookie 的弃用。
当这种情况发生时，处理效应在观察之间就不再是独立的了。实际上，如果一个地区的客户接受了处理，那么同一地区的其他客户也会接受处理。如果一个品牌的商品没有接受处理，那么同一品牌的其他商品也不会接受处理。在进行推断时，我们必须考虑这种依赖性：标准误差、置信区间和 p 值应该进行调整。在这篇文章中，我们将探讨如何使用聚类稳健标准误差来实现这一点。
客户、订单和标准误差
假设你是一个在线平台，并且你想要增加销售额。你刚刚有一个很好的想法：在结账时展示相关文章的轮播图来激励客户将其他文章添加到购物篮中。为了了解轮播图是否增加了销售额，你决定对其进行 A/B 测试。原则上，你可以为每个订单随机决定是否显示轮播图。然而，这会给客户带来不一致的体验，可能损害业务。因此，你决定在客户层面随机分配处理，即轮播图。对于接受处理的客户，你在每个订单上显示轮播图，而对于控制组客户，你永远不会显示轮播图。
我从src.dgp_collection导入数据生成过程（DGP），并从src.theme导入绘图主题。
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import numpy as np
from numpy.linalg import inv
import pandas as pd
from src.theme import *
from src.dgp_collection import dgp_carousel
我们为3000笔订单生成了模拟数据，来自100个客户。对于每笔订单，我们观察到order_id（也是数据集索引），下单客户的customer_id，结账时是否显示了carousel（处理），以及订单revenue（我们感兴趣的成果）。
n = 3000
dgp = dgp_carousel(n=n, w="carousel", y=["revenue"])
df = dgp.generate_data()
df.head()
数据快照，作者提供
由于处理是随机化的，我们可以通过比较处理订单的平均收入与控制（未处理）订单的平均收入来估计平均处理效果。随机化确保了处理订单和未处理订单在平均上是可以比较的，除了处理本身之外，因此我们可以将任何可观察到的差异归因于处理的因果效应。我们可以通过回归revenue到carousel虚拟变量来估计均值差异。报告的 OLS 系数是估计的平均处理效果。
import statsmodels.formula.api as smf
smf.ols("revenue ~ carousel", data=df).fit().summary().tables[1]
线性回归摘要表，作者提供
看起来包含carousel使得每笔订单的revenue减少了-1.63€。治疗效果似乎在*1%*的统计水平上具有显著性，因为报告的 p 值为0.005。然而，这不是一个标准的 AB 测试，因为我们没有对每笔订单进行随机化，而是对客户进行了随机化。同一个客户下的两笔订单不能分配到不同的处理组。因此，我们的治疗效果在观察之间是相关的，而我们在计算标准误差时假设它们是独立的。
报告的标准误差是否正确？我们如何“检查”它，我们能做些什么？
标准误差
哪些标准误差是“正确的”？
这个问题的答案取决于我们认为什么是随机的，什么又是固定的。在频率主义意义上，标准误差衡量的是“世界状态”或“平行宇宙”的不确定性，这些状态或宇宙是在我们看到了数据生成过程的随机部分的不同实现时可能发生的。
在这个特定的情况下，有一个变量是无可争议的随机变量：处理分配，实际上是我们自己随机化的。不仅如此，我们还知道它是如何随机的：每个客户有 50%的机会看到carousel，有 50%的机会看不到它。
因此，我们希望我们的估计标准误差能够衡量在不同治疗方案分配下估计的治疗效果的变异。这通常是一个抽象的概念，因为我们无法重新运行完全相同的实验。然而，在我们的案例中，我们可以做到这一点，因为我们处于一个模拟环境中。
DGP类有一个evaluate_f_redrawing_assignment函数，它允许我们通过重新采样数据并绘制不同的治疗方案来评估数据函数f(X)。我们想要评估的函数是我们用于治疗效应估计的线性回归。
def estimate_treatment_effect(df):
reg = smf.ols("revenue ~ carousel", data=df).fit()
return reg.params["carousel"], reg.bse["carousel"]
我们重复进行1000次不同的随机分配的治疗效应估计。
n_draws = 1_000
ols_coeff = dgp.evaluate_f_redrawing_assignment(f=estimate_treatment_effect, n_draws=n_draws)
在下面的图中，我们绘制了估计的 OLS 系数和标准误差的分布。在估计标准误差的图中（在右侧），我们还报告了模拟中估计系数标准差的一条垂直线（在左侧）。
def plot_coefficients(coeffs):
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))
# Estimates
coeffs_est = list(zip(*coeffs))[0]
sns.histplot(coeffs_est, color="C1", ax=ax1)
ax1.set(title=f"Distribution of estimated coefficients")
ax1.legend([f"Average: {np.mean(coeffs_est):.3f}nStd. Dev: {np.std(coeffs_est):.3f}"]);
# Standsed error
coeffs_std = list(zip(*coeffs))[1]
sns.histplot(coeffs_std, color="C0", ax=ax2)
ax2.set(title=f"Distribution of estimated std. errors")
ax2.legend([f"Average: {np.mean(coeffs_std):.3f}nStd. Dev: {np.std(coeffs_std):.3f}"]);
ax2.axvline(np.std(coeffs_est), lw=2, ls="--", c="C1")
plot_coefficients(ols_coeff)
模拟中回归系数和标准误差的分布，图片由作者提供
平均系数非常接近零。事实上，真正的治疗效果为零。然而，估计系数的标准差为2.536，与从线性回归表中得到的估计标准差0.576相当不同。这个标准差将意味着一个大约为0.5的p 值，与任何统计显著性阈值相差甚远。在右侧面板中，我们看到这并非偶然：标准 OLS 始终低估了系数的变异性，大约低估了 5 倍。
我们是否能够在现实中检测到这个问题，而没有重新随机分配治疗的可能性？是的，通过自助法计算标准误差。如果你从未听说过自助法，我写了一篇文章关于它，以及一种使自助法更快的技术：贝叶斯自助法。
自助法基本上包括有放回地抽取数据样本，并在自助样本中重新计算目标统计量。然后我们可以通过简单地计算样本间统计量的标准差来估计其标准差。
在这种情况下，重要的是要一致地根据数据生成过程来采样数据：按客户而不是按订单。
boot_estimates = []
customers = df.customer_id.unique()
for k in range(1000):
np.random.seed(k)
boot_customers = np.random.choice(customers, size=len(customers), replace=True)
df_boot = pd.concat([df[df.customer_id == id] for id in boot_customers])
reg = smf.ols("revenue ~ carousel", data=df_boot).fit()
boot_estimates += [(reg.params["carousel"], reg.bse["carousel"])]
在下面的图中，我们绘制了估计系数和标准误差的分布。标准误差仍然不正确，但自助估计的分布的标准差为2.465，非常接近模拟值2.536。
plot_coefficients(boot_estimates)
模拟过程中回归系数和标准误差的分布，图片由作者提供
注意，为了使用自助法检测线性回归标准误差中的问题，需要意识到分配是在客户级别进行的。在订单级别进行数据自助法仍然会低估标准误差。
聚类标准误差
线性回归表中报告的标准误差有什么问题？
问题在于，计算线性回归中标准误差的常用公式（稍后详述数学原理）假设观测值是独立同分布的，即 i.i.d.。然而，在我们的案例中，独立性假设不成立。为什么？
在我们的例子中，处理分配的单位，即一个客户，与观测的单位不同，即一个订单。这是一个问题，因为在不同的处理分配下，某个客户的全部订单要么被处理，要么不被处理。它们以块的形式移动，不可能出现同一客户的两个订单被分割在控制组和处理组之间。订单“以块移动”的后果是，我们的估计值比如果订单独立移动时可能存在的更多变异性。直观上，两组之间的平衡平均来说更少。因此，假设独立性计算出的标准误差通常低估了估计量的实际变异性。
有解决方案吗？是的！
解决方案是使用所谓的聚类稳健标准误差，它允许观测值在观测值的聚类内相关，在我们的案例中是一个客户。聚类稳健标准误差通常很容易实现，并且在所有统计软件包中都有提供。在使用statsmodels时，我们必须指定协方差类型为cluster，并且组别是按客户划分的。
smf.ols("revenue ~ carousel", data=df).fit(cov_type='cluster', cov_kwds={'groups': df["customer_id"]}).summary().tables[1]
带有聚类调整标准误差的回归表摘要，图片由作者提供
估计的聚类稳健标准误差等于 2.471，与模拟的标准误差 2.536 非常接近。请注意，估计的系数没有变化（之前为 -1.6292）。
聚类稳健标准误差是如何工作的？我们将在下一节深入探讨数学原理。
一些数学
OLS 估计量方差的通用公式是
OLS 估计量的方差，图片由作者提供
其中 X 是回归协变量，ε 是残差。该公式的关键组成部分是中心矩阵 n × n 矩阵 Ω = ε ε’，其中 n 是观测值的数量。这是关键，因为它是我们为了计算 OLS 估计量方差而需要估计的唯一对象。
起初，可能非常诱人只是使用回归残差e = y – ŷ来估计Ω，其中y是结果向量，ŷ是回归预测。然而，问题在于，根据构造，X和e的乘积为零，因此估计的方差将为零。
X.T @ e
array([[5.72555336e-11],
[6.68904931e-11]])
最简化的假设被称为**齐次性**：回归残差彼此独立，并且它们都具有相同的方差。在实践中，齐次性意味着Ω矩阵是对角线矩阵，每个单元格中的数字相同，σ²。
齐次性下的 Omega 矩阵，图片由作者提供
其中Iₙ是维度为n的单位矩阵。
在齐次性下，OLS 方差简化为
齐次性下的 OLS 系数方差，图片由作者提供
并且通过插入残差的方差来估计
y_hat = X @ inv(X.T @ X) @ (X.T @ y)
e = y - y_hat
np.var(e)
245.20230307247473
因此，OLS 估计量的估计方差由以下公式给出
齐次性下 OLS 系数的估计方差，图片由作者提供
np.var(e) * inv(X.T @ X)
array([[ 0.14595375, -0.14595375],
[-0.14595375, 0.33171307]])
估计的标准误差仅仅是右下角值的平方根。
np.sqrt(np.var(e) * inv(X.T @ X))[1,1]
0.5759453727032665
计算出的标准误差确实与之前线性回归表中报告的相同，0.576。
齐次性是一个非常严格的假设，通常会被放宽，允许观测值之间有不同的残差方差。这个假设被称为异方差性。
异方差性下的 Omega 矩阵，图片由作者提供
在异方差性下，OLS 估计量的方差公式不再简化。
Sigma = np.eye(len(df)) * (e @ e.T)
np.sqrt(inv(X.T @ X) @ X.T @ Sigma @ X @ inv(X.T @ X))[1,1]
0.5757989320663413
在我们的例子中，异方差性下的估计标准误差几乎相同：0.576。
在某些情况下，即使是异方差性的假设也可能过于严格，当观测值相关且回归残差不独立时。在这种情况下，我们可能希望允许Ω矩阵的某些非对角元素与零不同。但哪一些呢？
在许多情况下，包括我们的例子，假设观测值在某个集群内相关，但具有相同的集群内方差是合理的。例如，如果集群是观测值对，Ω矩阵具有以下形式。
集群相关下的 Omega 矩阵，图片由作者提供
我们现在可以计算在簇内相关性的情况下估计的标准误差。
customer_id = df[["customer_id"]].values
Sigma = (customer_id == customer_id.T) * (e @ e.T)
np.sqrt(inv(X.T @ X) @ X.T @ Sigma @ X @ inv(X.T @ X))[1,1]
2.458350214507729
我们确实得到了线性回归表中报告的相同数字。但这个公式的直观理解是什么？
直观理解
为了对估计的标准误差有一个直观的理解，想象一下我们有一个简单的回归模型，只有一个协变量：一个常数。在这种情况下，OLS 估计量的簇稳健方差就是每个簇内残差交叉乘积的总和。
无控制项的 OLS 估计量方差，图片由作者提供
其中 g 指代簇，ng 是簇 g 内的观测数。如果观测值独立，不同观测值残差乘积的期望值为零 𝔼[εᵢεⱼ]=0，估计的方差将接近残差平方和的总和。
无控制项和没有簇内相关性的 OLS 估计量方差，图片由作者提供
在另一个极端情况下，如果观测值高度相关，同一簇中的残差将非常相似，估计的方差将接近每个簇所有交叉乘积的总和。
无控制项和最大簇内相关性的 OLS 估计量方差，图片由作者提供
注意，这正是你在簇级别运行分析时获得的估计方差，在我们的例子中，是在客户级别汇总订单。
这意味着在实践中，如果你的分配单位与观测单位不同，你的有效样本量将介于实际样本量和簇数（在我们的例子中是客户数）之间。换句话说，你的观测数少于数据中的行数。你更接近哪个极端取决于残差的簇间相关性相对于簇内相关性的程度。
我们现在可以使用我们的数据生成过程来验证这个直观理解。首先，我们将残差的簇内方差缩放到零。
dgp = dgp_carousel(n=n, w="carousel", y=["revenue"])
dgp.revenue_per_customer_std = 0
df = dgp.generate_data()
现在未调整的和聚类标准误差应该是相同的，并且它们确实相当接近。
smf.ols("revenue ~ carousel", data=df).fit().summary().tables[1]
回归表总结，包含未调整的标准误差，图片由作者提供
smf.ols("revenue ~ carousel", data=df).fit(cov_type='cluster', cov_kwds={'groups': df["customer_id"]}).summary().tables[1]
回归表总结，包含聚类调整的标准误差，图片由作者提供
现在让我们将跨聚类的方差缩放到零。
dgp = dgp_carousel(n=n, w="carousel", y=["revenue"])
dgp.revenue_std = 0
df = dgp.generate_data()
在这种情况下，所有方差都在聚类层面，因此标准误差和聚类标准误差之间的差异将更加明显。
smf.ols("revenue ~ carousel", data=df).fit().summary().tables[1]
回归表总结，包含未调整的标准误差，图片由作者提供
smf.ols("revenue ~ carousel", data=df).fit(cov_type='cluster', cov_kwds={'groups': df["customer_id"]}).summary().tables[1]
回归表总结，包含聚类调整的标准误差，图片由作者提供
结论
在这篇文章中，我们看到了聚类稳健标准误差的重要性以及它们在随机实验中的相关性。如果你在比你的观察单位更高的层面上分配处理，这将在你的观察的处理效应之间产生相关性，使用假设独立性的常规公式计算标准误差可能会严重低估估计量的真实方差。我们探讨了两种解决方案：聚类稳健标准误差和分配单位的自助标准误差。
第三种保守的解决方案是在观察单位上聚合数据。这将在存在额外的非线性协变量时给出保守的标准误差估计。如果我们每个聚类有不同的观察数量，我们还需要使用回归权重。
聚类稳健标准误差的一个重要优点是，只有当观察之间存在实际依赖性时，它们才比通常的标准误差大。正如我们在上一节中看到的，如果观察在聚类之间只有轻微的相关性，聚类稳健标准误差将非常相似。
参考文献
- A. Abadie, S. Athey, G. W. Imbens, J. M. Wooldridge (2023), 何时应该调整聚类标准误差?.
相关文章
代码
你可以在以下链接找到原始的 Jupyter Notebook：
Blog-Posts/notebooks/clustering.ipynb at main · matteocourthoud/Blog-Posts
感谢阅读！
我真的很感激！ 🤗 *如果你喜欢这篇文章，并想看到更多，请考虑关注我。跟随我。我每周发布一次与因果推断和数据分析相关的话题。我尽量让我的文章简单但精确，总是提供代码、示例和模拟。
此外，一个小小的免责声明：我写文章是为了学习，所以错误是常态，尽管我尽力了。请，当你发现它们时，告诉我。我也很欣赏对新主题的建议！
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
