P值校正：Bonferroni与Benjamini-Hochberg方法详解

原创

用户11764306

发布于 2026-01-04 17:37:13

1410

Bonferroni与Benjamini-Hochberg：选择你的P值校正方法

P值可能是一个敏感的话题。或许初次与统计学家接触时最好避免讨论它。对这个话题的态度导致大家默认α = 0.05是黄金标准——实际上，这只是罗纳德·费舍尔本人设定的一个“方便的惯例”，一个经验法则。

他是谁？不认识他？别担心。他是第一个引入最大似然估计、方差分析和费舍尔信息的人。费舍尔是统计学界的核心人物，是统计学之父。他对孟德尔遗传学和进化生物学有深厚的兴趣，并为此做出了几项关键贡献。不幸的是，费舍尔也有棘手的过去。他与优生学学会及其对“智力低下者”的自愿绝育政策有关。

是的，著名的统计学家并非完美无缺。但由统计学之父设定的经验法则有时会被误认为是定律，而它并不是定律。

有一处关键情况，你不仅被允许，而且必须改变这个α水平，这完全取决于多重假设检验。

进行多次检验而不使用Bonferroni校正或Benjamini-Hochberg程序，问题不仅仅是麻烦。如果没有这些校正，我们可以证明任何假设：

H₁: 太阳是蓝色的

只需反复进行实验直到运气降临。但这些校正如何工作？你应该使用哪一个？它们不能互换！

P值和一个问题

要理解为什么，我们需要看看我们的P值究竟告诉了我们什么。要深入理解它，而不仅仅是“小就好，大就坏”。为此，我们需要一个实验，没有什么比发现超重元素更令人兴奋——或者更有争议的了。

这些元素极不稳定，在粒子加速器中一次一个原子地产生。按磅计算，这是有史以来生产的最昂贵的东西。它们只存在于超新星等宇宙事件中，持续时间仅千分之一或百万分之一秒。

但它们的不稳定性反而成为检测的优势，因为新的超重元素会表现出独特的放射性衰变。反应堆中传感器捕获的衰变序列可以告诉我们是否存在新元素。

作为我们的零假设，我们陈述：

H₀ = 该序列是背景噪声衰变。（没有新元素）

现在，如果我们想证明我们已经创造了一个新元素，就需要收集证据证明H₀不成立。这是通过我们的检验统计量T(X)完成的。一般来说，这捕获了传感器观察到的结果与背景辐射预期结果之间的差异。所有检验统计量都是衡量在H₀为真的情况下我们预期观察到的情况与样本数据实际显示情况之间的“差异”程度。T(X)越大，我们就有越多证据表明H₀是假的。

这正是Schmidt检验统计量对放射性衰变时间序列所做的。

σₒ₆ₛ = √1/(n-1) * Σᵢ₌₁ⁿ (ln tᵢ – ln t̄)²

Schmidt检验统计量被用于发现以下元素：钅黑、钅麦、钅鐽、錀、鎶，以及镆、鿬。

为H₀指定一个分布至关重要，这样我们才能计算检验统计量至少与观察数据计算出的检验统计量一样极端的概率。

我们假设噪声衰变遵循指数分布。有一百万个理由说明这是一个很好的假设，但我们不要在这里纠缠。如果我们没有H₀的分布，计算我们的概率值将是不可能的！

H₀(Schmidt): t₁, …, tₙ 独立同分布 ∼ Exp(λ)

那么，P值就是在零模型下，获得至少与从样本数据计算出的检验统计量一样极端的检验统计量的概率。我们的检验统计量越不可能，H₀为假的可能性就越大。

p = Pr_{H₀}(T(X) ≥ T(xₒ₆ₛ))

当然，这引出了一个有趣的问题。如果我们观察到一个罕见的背景衰变率，一个仅仅是类似于未发现的衰变粒子的衰变率呢？如果我们的传感器检测到一个不太可能但有可能的衰变序列，从而产生一个大的检验统计量呢？每次我们运行检验时，都有很小的机会仅凭运气就得到一个异常值。这个异常值将给出一个大的检验统计量，因为它与我们预期在H₀为真时看到的情况大不相同。大的T(x)将位于我们预期的H₀分布的尾部，并将产生一个小的P值。观察到比这个异常值更极端的事件的概率很小。但并没有新元素存在！我们只是玩了一百万次轮盘赌，得到了31次红色。

这看起来不太可能，但当你考虑到质子连续数月被射向目标粒子时，这种可能性是存在的。那么我们如何解释这一点呢？

有两种方法：一种保守的方法和一种不那么保守的方法。你的选择取决于实验。我们可以使用：

族错误率（FWER）和Bonferroni校正
错误发现率（FDR）和Benjamini-Hochberg程序

这些不能互换！你需要仔细考虑你的研究并选择正确的方法。

族错误率（Bonferroni）

这是我们的保守方法，如果我们不能容忍任何误报，就应该使用这种方法。这种方法将出现至少一个I类错误的概率保持在我们的α水平以下。

Pr（族中至少有一个I类错误） ≤ α

这也是一种更简单的校正。只需将α水平除以实验运行的次数。因此，对于每次检验，当且仅当满足以下条件时拒绝零假设：

pᵢ ≤ α/m

或者，你可以调整你的P值。如果你运行m次检验，取：

p_adjᵢ = min(1, m * pᵢ)

并且在满足以下条件时拒绝零假设：

p(Bonf)ᵢ ≤ α

我们在这里所做的只是将不等式两边都乘以m。

对此的证明也是一句简单的证明。如果我们让Aᵢ表示检验i中出现误报的事件。那么出现至少一个误报的概率将是所有这些事件并集的概率。

Pr（至少一个误报） = Pr(⋃ᵢ₌₁ᵐ Aᵢ) ≤ Σᵢ₌₁ᵐ Pr(Aᵢ) ≤ m * (α/m) = α

这里我们利用了并集界限。这是概率论中的一个基本概念，它指出A₁，或A₂，或Aₖ发生的概率必须小于或等于每个事件发生概率的总和。

Pr(A₁ ∪ A₂ ∪ ⋯ ∪ Aₖ) ≤ Σᵢ₌₁ᵏ Pr(Aᵢ)

错误发现率（Benjamini-Hochberg）

Benjamini-Hochberg程序也并不太复杂。简单地：

对你的P值排序：p₁ ≤ … ≤ pₘ。
接受第一个满足 pₖ > α/(m−k+1) 的k。

在这种方法中，目标是控制错误发现率（FDR）。

FDR = EV / max(R, 1)

其中R是我们拒绝零假设的次数，V是（不幸地）属于误报（I类错误）的拒绝次数。目标是将这个指标保持在特定阈值q = 0.05以下。

BH的阈值是：

(1/m)q, (2/m)q, …, (m/m)q = q

我们拒绝满足以下条件的最小的前k个P值：

P₍ₖ₎ ≤ (k/m)q

当你能够接受一些误报时，使用这种方法。当你主要关注的是最小化II类错误率时，也就是说，你希望确保有更少的漏报，即当H₀实际上是假的时候我们接受了H₀的情况。

把这想象成一项基因组学研究，你的目标是识别出每个人都拥有一个使其对某种特定癌症更易感的基因。治疗一些没有该基因的人，总比让拥有该基因的人得不到治疗而离开的危害要小。

快速对比

Bonferroni:

控制族错误率（FWER）。
保证出现单个错误发现的概率 ≤ α
更高的漏报率 ⇒ 更低的统计功效
零风险容忍度

Benjamini-Hochberg:

控制错误发现率（FDR）
保证在所有发现中，误报 ≤ q
更少的漏报 ⇒ 更高的统计功效
有一定的风险容忍度

超重原子的超小p值

我们的元素周期表中不能有任何不存在的元素，所以当涉及到寻找新元素时，Bonferroni校正是正确的方法。但是，当涉及到由位置敏感硅探测器收集的衰变链数据时，选择一个m并不那么简单。

物理学家倾向于使用在整个数据集上进行整个搜索所预期的随机链的数量：

Pr（≥1 条随机链）≈ 1 – e^{-n_b}

1 – e^{-n_b} ≤ α_family ⇒ n_b ≈ α_family （近似地，对于罕见事件）

随机链的数量来自于在没有进行实验时观察背景数据。从这些数据中，我们可以通过蒙特卡洛模拟建立零分布H₀。

我们通过对背景事件率建模并重新采样观察到的背景事件来估计随机链的数量。在H₀（无重元素衰变链）下，我们使用蒙特卡洛模拟许多零现实，并计算搜索算法产生与观察到的链一样极端的链的频率。

更精确地说：

H₀: 背景事件作为泊松过程到达，速率为λ ⇒ 到达间隔时间是指数分布的。

那么，一个意外链是τ时间内k个连续命中。我们使用我们的检验统计量扫描数据以确定是否存在极端聚类。

lambda_rate = 0.2   # 每秒事件数
T_total = 2_000.0   # 数据采集秒数（平均事件数 ~ 400）
k = 4               # 链长度
tau_obs = 0.20      # "观察到的极端值"：0.10秒内4个事件

Nmc = 20_000
rng = np.random.default_rng(0)

def dmin_and_count(times, k, tau):
    if times.size < k:
        return np.inf, 0
    spans = times[k-1:] - times[:-(k-1)]
    return float(np.min(spans)), int(np.sum(spans <= tau))
...

如果你对数字感兴趣，在发现第117号元素鿬时，使用了5×10⁻¹⁶的P值。我想，如果从未使用过任何校正，我们的元素周期表恐怕无法做成海报大小，化学也会陷入混乱。