计算 Statsmodels 中 A/B 测试的最小样本量：如何计算以及为什么

最新推荐文章于 2025-12-05 12:37:18 发布

原创最新推荐文章于 2025-12-05 12:37:18 发布 · 951 阅读

23 ·

CC 4.0 BY-SA版权

License CC BY-NC-SA 4.0 / 自豪地采用谷歌翻译

文章标签：

#数据库 #前端 #linux

榛樿鍒嗙被专栏收录该内容

791 篇文章

订阅专栏

原文：towardsdatascience.com/computing-minimum-sample-size-for-a-b-tests-in-statsmodels-how-and-why-398e357945d9?source=collection_archive---------3-----------------------#2024-05-31

深入探讨 Statsmodels 如何以及为什么使用数值优化而非封闭形式的公式

https://jasonjiajs.medium.com/?source=post_page---byline--398e357945d9--------------------------------https://towardsdatascience.com/?source=post_page---byline--398e357945d9-------------------------------- Jason Jia

·发表于Towards Data Science ·阅读时间 10 分钟·2024 年 5 月 31 日

–

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/cfc8f9479bd992d2d46961c47bd201c2.png

由 DALL-E 生成的图像

介绍

目前没有关于 Statsmodels 如何计算最小样本量的好资源。

在进行 A/B 测试之前，计算所需的最小样本量至关重要。一种常用的方法是在 Python 的 Statsmodels 包中调用tt_ind_solve_power 函数，但目前在理解它的工作原理时存在两个空白：

有很多优秀的文章（例如 Stan Nsky，TDS 2019）解释了参数的含义并提供了函数调用的示例。然而，它们并没有解释该函数是如何实际计算样本量的，也没有说明为什么该过程是正确的。
也有许多优秀的文章（例如 Mintao Wei，TDS 2023）解释了基于 z 检验的比例统计推导，例如转换率，这是许多在线样本量计算器（例如Evan Miller 的计算器）也常用的方法。然而，这并不是 Statsmodels 使用的方法，因此结果可能会有所不同。

这对数据科学家来说很重要，因为 Statsmodels 常用于 Python 中计算样本量。

数据科学家经常使用 Statsmodels 来获取最小样本量，但可能不知道它采用了与大多数文章描述的方式以及大多数在线计算器使用的方法不同的方式。了解该函数的工作原理是至关重要的，以便我们能够信任其结果。

本文通过解释 Statsmodels 的实际工作原理来填补这一空白。

本文旨在做出创新性贡献，解释 tt_ind_solve_power 如何实际计算样本量，为什么该过程是正确的，并且它相对于封闭形式解法带来的优势。[1]

第一部分： 它将首先解释样本量如何计算以及为什么该程序是正确的，分为两个步骤：

显示样本量计算的统计推导。
编写一个简化版的tt_ind_solve_power，它是统计推导的精确实现，并且产生与原始函数相同的输出。

第二部分： 接下来，它将解释它相对于封闭形式解法带来的两个优势：

泛化性优势
统计直觉的优势

第一部分：Statsmodels 如何计算最小样本量以及为什么它是正确的

1.1. 显示样本量计算的统计推导

核心思想

一般的 A/B 测试是一个无配对的双样本 t 检验。Statsmodels 并不使用封闭形式的解法，而是通过两个步骤来获得最小样本量：

对于给定的样本量，计算该检验的功效。
运行数值优化算法，找到返回目标功效的样本量。

符号和概念

以下是我们将在本文中使用的一些术语：

n：最小所需样本量。n = n_1 + n_2
n_1, n_2：治疗组和对照组的最小所需样本量
比例： n_2 = n_1 * 比例，对于 50:50 的分配，比例 = 1
p：p 值
𝛼: 显著性水平 / I 类错误
𝛽: II 类错误；1-𝛽是检验的功效
μ_1, μ_2：治疗组和对照组的均值
X̄1, X̄2: 治疗组和对照组的样本均值
t_(1-𝛼): 临界值 / t-score，截断标准 t 分布的上 100𝛼(%)部分。
MDE: 最小可检测效应，或给定所有其他参数下可以检测到的统计显著性差异的水平（例如，基础转换率为 10%，预期提升为 50%，因此预期的治疗转换率为 15%，意味着 MDE 为 15-10=5%=0.05）
𝜎: 每组观察值的标准差，假设是相同的
d: Cohen’s d / 标准化效应大小，由 MDE / 𝜎 给出
H_0, H_1: 原假设，备择假设

推导检验功效的公式

定义原假设和备择假设：

2. 推导原假设下检验统计量的分布（H_0）：

我们发现在原假设下，检验统计量 t 服从t 分布，其自由度为(n_1 + n_2 - 2)。

这可以从以下内容得出：

其中 X 的样本方差计算如下：

3. 推导备择假设(H_1)下检验统计量的分布：

我们发现在备择假设下，假设均值差异为最小可检测效应（MDE），检验统计量 t 服从非中心 t 分布，其非中心性参数θ = d * sqrt((n1 * n2) / (n1 + n2))，自由度为**(n_1 + n_2 - 2)**。

一个具有正非中心性参数的非中心 t 分布（nct）可以粗略地看作是一个标准 t 分布向右平移的结果。[2] 直观地说，标准 t 分布发生在原假设下，当我们期望平均效果为 0 时，而非中心 t 分布发生在备择假设下，当我们期望一个正的效果，并且该效果的平均值大致等于 MDE。

定义：具有非中心性参数θ和自由度ν的非中心 t 分布随机变量 T 定义如下：

其中 Z 是标准正态随机变量，V是自由度为ν的卡方分布随机变量。

证明从以下观察开始：在备择假设下，真实的均值差异为 MDE，因此我们可以减去 MDE 并除以总体标准差，从而得到标准正态变量。

4. 计算功效

既然我们知道了原假设和备择假设下检验统计量的分布，并且两种分布的累积分布函数（cdf）都是已知的，我们可以轻松地计算功效，给定显著性水平和检验类型（双尾、大于、小于）。下图可视化了这一过程：

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/f7537362230bff4e8f958aaa8ca1cac5.png

作者提供的图表

在 Python 中，实现如下：

def power(self, effect_size, nobs1, alpha, ratio=1, df=None,
          alternative='two-sided'):
    nobs2 = nobs1*ratio
    if df is None:
        df = (nobs1 + nobs2 - 2)

    # Get non-centrality parameter
    nobs = nobs1 * nobs2 / (nobs1 + nobs2)
    d = effect_size
    nc_param = d * np.sqrt(nobs)

    # Get effective level of signifiance, alpha_
    if alternative in ['two-sided']:
        alpha_ = alpha / 2.
    elif alternative in ['smaller', 'larger']:
        alpha_ = alpha
    else:
        raise ValueError("alternative has to be 'two-sided', 'larger' " +
                        "or 'smaller'")

    # Compute power of a t-test
    power = 0
    if alternative in ['two-sided', 'larger']:
        crit_upp = stats.t.isf(alpha_, df) # isf = inverse survival function = value where Pr(t > value) = alpha
        power += 1 - special.nctdtr(df, nc_param, crit_upp)  # 1 - Pr(t < crit_upp) = Pr(t > crit_upp) for non-central t distribution
    if alternative in ['two-sided', 'smaller']:
        crit_low = stats.t.ppf(alpha_, df) # ppf = percent point function = value where Pr(t < value) = alpha
        power += special.nctdtr(df, nc_param, crit_low) # 
    return power

通过数值优化获取最小样本量

既然我们现在知道了如何为一组给定的参数计算功效，我们就可以运行数值优化方法来找到实现目标功效的最小样本量。由于总样本量是治疗组样本量的函数（n = n_1 + 比例 * n_1），我们将求解 n_1。

之所以有效，是因为功效随着样本量 n_1 的增加而单调增加。直观地说，更多的样本意味着 A/B 测试结果更加可靠，因此，如果备择假设为真，更多的值将拒绝原假设（见下图左侧子图）。

但这也意味着，减去目标功效后，得到的是一个单调递增的函数，起点为负，终点为正。根据介值定理和函数的单调性，存在一个唯一的根，对应于我们的最小样本量（见下图右侧子图）。

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/c66e6dffb8db57745b2cf9d9a2543f10.png

作者提供的图表

一种流行且高效的数值优化方法是Brent 方法。Brent 方法是一种根求解算法，结合了如二分法、割线法和反二次插值等多种技术。有关其在 Statsmodels 中实现的更多细节，可以参见此处。

在 Python 中，实现如下：

def solve_power(self, effect_size=None, nobs1=None, alpha=None, power=None,
                ratio=1., alternative='two-sided'):
    print('--- Arguments: ---')
    print('effect_size:', effect_size, 'nobs1:', nobs1, 'alpha:', alpha, 'power:', power, 'ratio:', ratio, 'alternative:', alternative, '\n')

    # Check that only nobs1 is None
    kwds = dict(effect_size=effect_size, nobs1=nobs1, alpha=alpha,
                power=power, ratio=ratio, alternative=alternative)
    key = [k for k,v in kwds.items() if v is None]
    assert(key == ['nobs1'])

    # Check that the effect_size is not 0
    if kwds['effect_size'] == 0:
        raise ValueError('Cannot detect an effect-size of 0\. Try changing your effect-size.')

    # Initialize the counter
    self._counter = 0

    # Define the function that we want to find the root of
    # We want to find nobs1 s.t. current power = target power, i.e. current power - target power = 0
    # So func = current power - target power
    def func(x):
        kwds['nobs1'] = x
        target_power = kwds.pop('power') # always the same target power specified in keywords, e.g. 0.8
        current_power = self.power(**kwds) # current power given the current nobs1, note that self.power does not have power as an argument
        kwds['power'] = target_power # add back power to kwds

        fval = current_power - target_power
        print(f'Iteration {self._counter}: nobs1 = {x}, current power - target power = {fval}')
        self._counter += 1
        return fval

    # Get the starting values for nobs1, given the brentq_expanding algorithm
    # In the original code, this is the self.start_bqexp dictionary set up in the __init__ method
    bqexp_fit_kwds = {'low': 2., 'start_upp': 50.}

    # Solve for nobs1 using brentq_expanding
    print('--- Solving for optimal nobs1: ---')
    val, _ = brentq_expanding(func, full_output=True, **bqexp_fit_kwds)

    return val