自助法和置换测试
介绍
大多数与数据工作的人都会观察数据，然后想知道这些观察结果是否具有统计学意义。除非有些统计推断的正式培训和运行显著性测试的经验，否则人们首先会想到找一个能提供如何进行测试的建议，或者至少能确认测试已经正确执行并且结果是有效的统计学家。
这种情况有很多原因。首先，通常不会立即明确需要哪种测试，哪些公式支持测试原理，如何使用这些公式，以及测试是否可以使用，例如数据是否满足正态性等。有许多全面的R和Python包用于估计丰富的统计模型和进行统计测试，比如statsmodels。
然而，如果没有对统计理论的充分理解，只是通过复制用户指南中的示例来使用包，往往会让人感到不安，因为预计会受到经验丰富的统计学家的严厉批评。就我个人而言，我是一个随着时间变成数据分析师的工程师。在我的本科和研究生学习期间，我上过统计课程，但是我并没有广泛使用统计学，因为这通常不是工程师的职业。我相信对于许多其他数据分析师和数据科学家也是一样的，特别是如果他们的正式培训是工程、计算机科学或化学等领域。
我决定写这篇文章是因为我最近意识到，模拟可以很容易地代替更经典的基于公式的统计方法。大多数人可能会立即想到用自助法估计均值的不确定性。但问题不仅仅是自助法。在随机置换测试中使用重抽样可以回答许多统计推断问题。这些测试通常不难编写和执行。它们适用于连续或二元数据，不受样本大小限制，并且不对数据分布做假设。从这个意义上说，置换测试是非参数的，唯一的要求是可交换性，即观察到某个值序列的概率对于序列的任何置换都是相同的。这实际上要求不多。
过去，计算资源的不可用可能是基于公式的统计推断测试取得令人印象深刻进展的原因之一。在那个时候，对具有几十或几千个记录的数据样本进行数千次重抽样是禁止的，但现在不再是禁止的了。这是否意味着不再需要经典的统计推断方法？当然不是。但是，能够运行置换测试并确认结果在相似性上可以提供一些安全感，或者帮助理解当我们观察到差异时哪些假设不成立。能够从头开始运行统计测试而不依赖于一个包，也给人一些掌控感。
当然，置换测试并不是什么新鲜事物，但我认为提供一些示例和相应的代码是一个好主意。这可能会减轻一些数据专家的恐惧，并将使用模拟进行统计推断更加贴近他们的日常实践。本文使用置换测试来回答两个问题。在许多其他情况下也可以使用置换测试，并且对于更复杂的问题，设计置换测试可能不太明显。从这个意义上说，本文并不全面。然而，原理是相同的。通过理解基本原理，更容易查找权威来源，了解如何设计置换测试来回答其他更微妙的业务问题。我的意图是激发一种思考方式，即将模拟人群分布放在中心位置，使用理论抽样来估计观察到的效果发生的随机概率。这就是假设检验的含义。
统计推断从一个假设开始，例如，新药对某种疾病的疗效比传统治疗更好。疗效可以通过检查给定血液指标的减少（连续变量）或者在使用新药和传统治疗（对照）时无法检测到疾病的动物数量（离散变量）来衡量。这种两组比较，也被称为A/B测试，在所有经典统计学教材和流行的技术博客中都有详细讨论，比如这篇文章。利用药物设计的例子，我们将测试新药是否比传统治疗更有效（A/B测试）。在此基础上，我们将估计需要多少动物才能证明新药比传统治疗实际上更有效，假设新药比传统治疗效果更好1%（或其他效果大小）。虽然这两个问题似乎没有关联，但实际上它们是有关联的。我们将从第一个问题中重用代码来回答第二个问题。所有的代码都可以在我的博客代码库中找到。
欢迎留言，但请给予建设性的评论。我并不自称为统计学家，我的意图是帮助他人在进行排列测试时经历类似的学习过程。
A/B测试
让我们回到第一个问题，即新药是否比传统治疗更有效。当我们进行实验时，患病的动物被随机分配到两个组，取决于它们接受的治疗方式。动物被随机分配到组别，因此治疗效果的任何观察差异都是由于药物的有效性，或者只是偶然发生的强免疫系统的动物被分配到新药组。这是我们需要解开的两种情况。换句话说，我们想要检查随机机会是否能够解释使用新药的任何观察到的益处。
让我们假设一些虚构的数字来进行说明：
响应变量是二进制的，即治疗成功与否。排列测试的工作方式与响应变量是连续的情况相同（这不适用于经典统计测试！），但上面的表格将包含平均值和标准差，而不是计数。
我们有意不使用相同大小的治疗组，因为这不是排列测试的要求。这个假想的A/B测试涉及大量的动物，而且新药看起来很有希望。新药比传统治疗更有效率1.5%。鉴于样本较大，这看起来是显著的。我们将回到这个问题。作为人类，我们倾向于认为那些可能并不显著的事情是显著的。这就是为什么标准化假设测试如此重要。
“把零假设想象成什么都没有发生，即机会可以解释一切。”
在A/B测试中，我们使用一个基准假设，即没有观察到特殊的情况。这也被称为零假设。进行测试的人通常希望证明零假设不成立，即已经做出了一个发现。换句话说，备择假设是真实的。证明这一点的一种方法是显示随机机会导致与观察到的差异一样极端的差异的概率非常低。我们已经开始看到与排列测试的联系。
想象一个过程，将所有接受治疗的动物汇集成一个单一的组（由2487 + 1785只动物组成），然后再次随机分成两个与原始治疗组大小相同的组。对于每只动物，我们知道治疗是否成功，因此我们可以计算每个组治愈动物的百分比。使用观察到的数据，我们确定新药将治愈动物的百分比从80.34%提高到81.79%，即增加了近1.5%。如果我们多次对两组进行重新抽样，我们会看到新药导致更高百分比的动物被治愈相对于传统治疗的情况有多常见？这个“多常见”就是统计推断中无处不在的p值。如果这种情况经常发生，即p值大于我们可以接受的阈值（通常是5%的显著水平），那么我们在实验中观察到的情况可能是由于偶然性，因此零假设不被拒绝。如果这种情况很少见，那么仅仅靠偶然性就无法解释观察到的差异，因此零假设被拒绝（如果你的团队发现了新药，你可以举办一场派对！）。如果你仔细观察我们实际上使用排列所做的是模拟零假设，即两个治疗组是等效的。
再次考虑一下零假设是如何被制定的，因为这将决定如何进行排列测试。在上面的例子中，我们想要看到多少次机会会让我们相信备择假设是真实的，即新药更有效。这意味着零假设，作为备择假设的补充，说明新药的效果不如传统治疗，或者与传统治疗一样有效。这也被称为单侧检验（与双侧检验相反，也被称为双向检验）。从另一个角度来看。我们不希望被随机机会欺骗，认为新药更有效。被另一方面欺骗并不重要，因为我们本来就不打算替代传统治疗。双侧检验会导致较高的p值，因此更保守，因为它拒绝零假设的机会更大。然而，这并不意味着如果不是正确的测试就应该使用它。
排列测试可以在最一般的情况下表述如下。假设有Gᵢ，i=1,..,Nᴳ个具有基数∣ Gᵢ ∣，i=1,..,Nᴳ的组：
- 将所有组的数据点汇集在一起；这本质上模拟了零假设，即假设没有发生任何事情。
- 无替换地将∣ G₁ ∣个点随机分配给组G₁，无替换地将∣ G₂ ∣个点分配给组G₂，..，直到所有点都被分配完。
- 计算感兴趣的统计量，就像在原始样本中计算的那样，并记录结果。
- 重复上述过程很多次，并记录每次的感兴趣统计量。
基本上，上述过程建立了一个具有感兴趣统计量的分布。如果观察到的差异至少与观察到的差异一样极端的值的概率是p值。如果p值很大，那么机会很容易产生观察到的差异，我们还没有取得任何发现（尚未）。
“将p值视为如果零假设成立，则观察到与我们的观察一样极端结果的概率。”
上述公式非常通用。回到我们的例子，我们只有两组，一个是新药组，一个是传统治疗组。执行排列检验的代码如下。
我们进行了10,000次排列，大约需要30秒。关键问题是：机会有多少次使新药的效果比传统治疗提高1.5%或更多？我们可以将模拟的效果差异的直方图可视化，并计算p值，如下所示。
这给出了以下直方图：
红色柱表示通过机会发现新药比传统治疗更有效。这似乎并不罕见。p值为0.1084。假设我们希望以0.05的显著性水平运行测试，这意味着无法拒绝零假设。此时没有什么好庆祝的。如果你组织了一个派对，它需要被取消。或者可能推迟。
“将a视为虚警率，即假设零假设为真，我们会在重复实验时有5%的时间得出统计上显著的差异。”
有一些理由可以持乐观态度。我们刚刚进行的A/B测试有两种可能的结果：要么存在效果（在我们的例子中，新药比传统治疗更有效），要么没有足够的证据来得出结论存在效果。测试并未得出不存在效果的结论。新药可能确实更有效。但是到目前为止，根据已有的数据，我们无法证明这一点。该测试本质上保护我们免受错误阳性（也称为第一类错误）的影响；但是可能存在错误阴性（也称为第二类错误）。这是团队所希望的。
我们还可以提出另一个问题。要得出新药比传统治疗更有效的结论，观察到的差异需要多大？显然1.5%是不够的，但是需要多少才足够？可以从生成的直方图中很容易地获得答案。我们可以将对应于观察到的差异的垂直线向右移动，直到红色柱子的尾部占据总面积的5%；换句话说，使用95%百分位数np.percentile(differences, 95)，得到0.0203或2.03%。比我们观察到的1.5%多一点，但是并不太离谱。
使用0.05的显著性水平，如果新药在治疗效果上的增加在区间（-∞，0.0203]内，我们将不会拒绝零假设。这也被称为置信区间：观察到的统计量值的一组值，不会拒绝零假设。因为我们使用了5%的显著性水平，所以这是一个95%的置信区间。假设新药不更有效，那么多次运行实验会在置信区间内95%的时间内给出效果差异。这就是置信区间告诉我们的。当且仅当置信区间包含观察到的效果增加时，p值将超过a，这意味着无法拒绝零假设。这两种检查零假设是否可以被拒绝的方法当然是等效的。
根据目前测试的动物数量，我们无法拒绝零假设，但离置信区间边界也不远。团队持乐观态度，但我们需要收集更有说服力的证据来证明新药的有效性更高。但需要多少更多的证据呢？在下一节中，我们将重新讨论这个问题，因为使用重采样进行模拟可以帮助我们回答这个问题！
在结束本节之前，重要的是要注意，我们也可以使用经典的统计检验来近似计算p值。上面呈现的表格也称为列联表，它提供了两个变量之间的相互关系，并可用于确定它们之间是否存在相互作用。可以使用卡方检验从列联矩阵开始检查两个变量的独立性，但需要注意不要运行双侧检验（没有进行过广泛尝试，但scipy似乎默认使用双侧检验；这将导致较高的p值）。在深入研究统计库的用户指南之前，了解如何运行置换检验是不是很好呢？
功效估计
当然，如果我们无法证明新药的增效是具有统计显著性的，人们会感到失望。毕竟，新药可能真的更好。我们愿意通过治疗更多的动物来做更多的工作，但我们需要多少动物呢？这就是功效的作用。
功效是在给定样本大小和显著性水平的情况下检测到给定效应大小的概率。假设我们预计新药相较于传统治疗能够将治疗效果提高1.5%。假设我们已经对每种治疗方法分别治疗了3000只动物，并将显著性水平固定为0.05，那么检验的功效为80%。这意味着如果我们重复进行实验很多次，我们会发现在5次实验中有4次我们会得出结论：新药比传统治疗更有效。换句话说，错误否定（第二类错误）的概率为20%。以上数字当然是假设的。重要的是这四个量：样本大小、效应大小、显著性水平和功效是相关的，并且设置其中的三个量可以计算出第四个量。最典型的情况是从其他三个量中计算样本大小。这就是我们在本节中要研究的内容。为了简化，我们假设在每次实验中，我们对新药和传统治疗的动物数量是相同的。
下面的过程试图构建一个以样本大小为函数的功效曲线：
- 创建一个合成数据集，其中的动物假设已经接受了传统治疗，以便治疗效果与我们所知道的差不多（下面，我将其设置为0.8034，与上面的列联矩阵相对应）。
- 创建一个合成数据集，其中的动物假设已经接受了新药治疗，通过添加我们想要研究的效应大小（下面，我将其设置为0.015和0.020，以观察其对结果的影响）。
- 从每个合成数据集中绘制大小为n_sample的自助样本（下面我将其设置为3000、4000、5000、6000和7000）。
- 使用我们在前一节中确定的方法进行置换检验以检验治疗效果的差异是否具有统计显著性，并记录结果。
- 不断生成自助样本并计算治疗效果的差异具有统计显著性的频率；这就是检验的功效。
当然，这是一个更长的模拟过程，因此我们将自助样本的数量限制在200个，而显著性检验中的置换次数也与前一节相比减少到500次。
在一台普通的机器上运行这个自助采样/置换模拟需要大约一个小时，并且可以从多处理中获益，但这超出了本文的范围。我们可以使用matplotlib直观地显示结果：
这将生成以下图表：
我们从中学到了什么？如果我们预计新药的效果提高了1.5%，那么为了以80%的功效证明这一点，我们需要治疗超过7000只动物。如果效应大小更大，即2%，那么我们需要做的工作就少了，大约4500只动物就足够了。这是直观的。检测到大效应比小效应更容易。决定进行如此大规模的实验需要进行成本/效益分析，但至少现在我们知道证明新药更有效需要什么。
我们还可以使用statsmodels来计算所需的样本量：
这将打印：
效应大小：0.015，样本量：8426.09
效应大小：0.020，样本量：4690.38模拟结果似乎一致。在模拟中，我们将样本量增加到7000，但当效应大小为1.5％时，无法达到0.8的功效，这也可以通过使用proportion_effectsize函数看到。
总结
希望您喜欢这篇文章。个人而言，我觉得使用简单的自助法和排列法从头开始研究所有这些统计概念非常有意义。
在结束之前，需要注意一点。本文非常强调p值，但p值越来越受到批评。事实上，p值的重要性在历史上被夸大了。p值表示数据与代表零假设的统计模型或排列测试的不兼容程度。p值并不是替代假设为真的概率。此外，一个能够拒绝零值的p值并不意味着效应大小重要。一个很小的效应大小可能在统计上是显著的，但它是如此之小，以至于并不重要。
参考文献
- 《入门统计与分析：重新采样视角》（Peter Bruce著，Wiley，2014年）
- 《数据科学家的实用统计学：使用R和Python的50个+基本概念》（Peter Bruce、Andrew Bruce、Peter Gedeck著，O’Reilly，2020年）
- 《解读A/B测试结果：误报和统计显著性》（Netflix技术博客）