16 重抽样与多重比较

抽样分布有时可以从模型直接推导，有时却难以写成便于计算的形式。计算机重复生成样本，为后一类问题提供了近似途径。Bootstrap 用经验分布替代未知总体，研究统计量随抽样怎样变化；置换检验则在零假设允许的交换中，研究观测结果有多极端。两种方法生成数据的依据不同。

同时研究多个假设还带来另一种概率变化：关注的事件从单项错误，转为一组判断中的错误。单项检验的控制水平无法直接代表成组控制水平。模拟次数、可交换性、检验家族和依赖结构，分别影响这些近似与控制的有效性。

这些计算仍以抽样分布为对象：置信区间利用估计量的波动，假设检验与错误率则利用零假设下的参照分布。下面的 Python 代码将相应概率机制写成重复计算；即使不熟悉代码，也可以沿生成样本、计算统计量和比较结果的过程理解。

非参数 Bootstrap

设已有 $n$ 个独立同分布观测 $x_1,\dots ,x_n$，目标是估计统计量 $T$ 的波动，例如均值、中位数或两组风险差的波动。若能从真实总体反复抽样，就能直接观察 $T$ 的抽样分布；但当前只有一份样本。非参数 bootstrap 将每个观测值赋予 $1/n$ 的概率，构成经验分布，再从这个近似总体中有放回地抽 $n$ 个观测。

有放回抽样使一份新样本可以重复包含某个观测，同时遗漏另一个。组成的变化产生统计量的变化，这就是模拟抽样波动的来源。若无放回抽完全部 $n$ 人，得到的只是原样本的重新排序，均值等统计量不会变化。

重复 $B$ 次，便得到 $T_1^{\ast },\dots ,T_B^{\ast }$。这组数值的样本标准差估计原统计量的标准误，第 $2.5\%$ 与第 $97.5\%$ 分位点则给出简单的 percentile bootstrap 区间。由于经验分布只是总体的近似，这种区间不具有自动精确的覆盖率。偏差较大、分布偏斜或参数位于边界时，BCa、学生化区间等方法可作不同校正，但也各有适用条件。

两独立组风险差的 Bootstrap 计算

虚构实验中，对照组 $40/200$，干预组 $70/200$，原始风险差为 $0.15$。两组独立，因此各有自己的经验分布。一次重抽在对照组内生成 $200$ 人，在干预组内生成 $200$ 人，再计算比例差；这样保留了原设计的组别及样本量，也允许两组成功率不同。

代码中的虚构 $0/1$ 数组逐项对应独立受试者。随机种子固定模拟序列，使这一次数值近似可以复现；它不会改变原研究所包含的受试者信息量。

import numpy as np
 
rng = np.random.default_rng(20261012)
y1 = np.r_[np.ones(70), np.zeros(130)]
y0 = np.r_[np.ones(40), np.zeros(160)]
 
rd_star = np.array([
    rng.choice(y1, size=200, replace=True).mean()
    - rng.choice(y0, size=200, replace=True).mean()
    for _ in range(5000)
])
print(np.quantile(rd_star, [0.025, 0.975]))

所得区间约为 $[0.065,0.235]$，与本例的大样本区间接近。这一接近只说明两种近似在当前样本中的数值关系，其他样本和统计量仍有各自的适用条件。增加 $B$ 减少的是计算中的模拟误差；原研究仍只有 $400$ 人，不会产生 $400\times B$ 位独立受试者的信息。

重抽样单位与依赖结构

经验分布的构造单位决定重抽是否保留原有联合结构。同一患者的前后测量组成一对，整对重抽才能保留配对关系；同一医院内的患者若相关，合适的聚类单位重抽可以保留这种依赖；时间序列的区块方法则保留局部相关。逐条独立重抽相关记录可能低估标准误。复杂抽样中的分层、权重和设计信息同样属于原抽样结构。

经验分布的另一个限制是只能生成已出现过的值。$10$ 人零次成功时，从十个零中重抽仍然全零，percentile 区间于是退化成 $[0,0]$。有限样本未出现成功，并没有证明总体成功率为零；经验分布在这里缺少稀有事件的信息。合适的二项区间或模型方法能够利用不同的概率结构表达这份不确定性。

许多平滑统计量具有良好的 bootstrap 近似性质，但极值、极稀疏类别、分离回归和经过强选择的统计量需要额外的理论依据与诊断。重抽样也没有引入原样本之外的人群信息，因此无法修复样本的代表性偏差。

置换检验与交换性

Bootstrap 从估计出的经验分布生成新样本，保留了各组当前的分布差异。置换检验则固定已有观测，在零假设允许的交换中重新分配标签。前者近似估计量的波动，后者构造零假设参照，二者因而不必产生相同的分布。

若两组来自同一分布、观测可交换，将 $400$ 个结果合并后，任意选 $200$ 个作为干预组、其余作为对照组，就得到零假设允许的一种组成。每种组成产生一个风险差，绝对差异至少达到观测值 $0.15$ 的比例近似双侧尾概率。

这种交换之所以给出合法参照，是因为可交换性：在零假设下，允许的交换保持联合分布不变。独立同分布样本且两组完整分布相同，是经典情形。随机试验中，参照还可以直接来自实际随机分配机制；在“每个人的结果在任何分配下都不改变”的尖锐无效应假设下，个体结果固定，重新随机分配便产生相应的随机化分布。

均值相同并不规定完整分布相同。连续数据若均值相同而方差不同，原始观测通常不再满足简单标签交换所需的条件，未标准化均值差的置换不能一般地保证精确性。适当学生化的方法可能获得大样本有效性。观察性分组若与年龄等因素有关，也不能仅凭标签能够交换，就推定其联合分布在交换后保持不变。

配对实验的允许分配由各对内部的标签交换组成，分层随机实验的允许分配则由各层内部的随机化组成。合并全部个体后任意交换标签，会引入原设计不允许的分配，参照分布也随之改变。

有限次随机置换的 $p$ 值

当全部允许分配过多时，可以随机抽取其中一部分。设共进行 $B$ 次随机置换，其中 $b$ 次至少与观测值同样极端，常用的 Monte Carlo 估计为

$${\hat{p}}_{\mathrm{MC}}=\frac{b+1}{B+1}.$$

这个形式把观测分配也纳入排名，因此不会因有限次模拟没有命中尾部就得到 $p=0$。尾部包括与观测相等的统计量；代码中的微小数值容差用于使浮点计算中的相等情况也得到计数。

rng = np.random.default_rng(20261012)
y = np.r_[y1, y0]
observed = abs(y1.mean() - y0.mean())
b = 0
B = 9999
for _ in range(B):
    perm = rng.permutation(y)
    stat = abs(perm[:200].mean() - perm[200:].mean())
    b += stat >= observed - 1e-12
print((b + 1) / (B + 1))

上述模拟有 $b=7$ 次落入尾部，得到 $0.0008$；完整条件分布给出的相应双侧概率约为 $0.001111$。有限次模拟本身会波动，所以两者可以不同。当真实 $p$ 约为 $0.001$ 时，一万次独立随机置换的模拟标准误仍约为

$$\sqrt{\frac{0.001(1-0.001)}{10000}}\approx \mathrm{0.000316.}$$

这个标准误相对于目标概率已经不小，尤其会影响临近决策阈值或很小的 $p$ 值的数值判断。增加置换次数可以降低这种误差，直接使用可计算的精确分布则避免了这项模拟误差；打印更多小数位只改变显示精度。

多重比较与成组错误率

设 $20$ 个检验相互独立，全部零假设均为真，每项实际拒绝概率恰为 $0.05$，例如理想的连续校准检验。一项也不误拒绝的概率是各项不误拒绝概率的乘积，其补事件给出至少一个假阳性的概率：

$$1-(1-0.05{)}^{20}\approx \mathrm{0.642.}$$

单项检验的水平保持有效，成组出现至少一次误判的概率却已约为 $64\%$。若每项实际错误率仅保证“不超过 $0.05$”，例如离散或保守检验，这个数在独立条件下是上界；检验相关时，乘法分解也不再直接成立。

家族错误率 FWER 关注一组中至少一次错误拒绝；错误发现率 FDR 关注全部发现中错误所占比例的期望。设错误拒绝数为 $V$，全部拒绝数为 $R$，则

$$\mathrm{FDR}=E\left[\frac{V}{\max (R,1)}\right].$$

分母的约定使没有发现时错误发现比例为零。FWER 对任何一次错误都敏感，FDR 则容许错误数随发现数变化，因此探索大量候选指标时可能采用 FDR，少数关键确认性结论常更重视 FWER。控制目标的差别来自研究对错误所作的不同定义。

Bonferroni 方法

对于 $m$ 个检验，Bonferroni 方法以 $p_i\le \alpha /m$ 作为每项的拒绝标准。至少一个错误是各个错误事件的并集，其概率不超过各真零假设误判概率之和，所以 FWER 不超过 $\alpha$。这一上界无需独立性，但每个输入 $p$ 值仍须具有有效的零假设校准。

例如 $m=6$、家族 $\alpha =0.05$，每项阈值为 $0.00833$。并集上界可能高于实际并集概率，因此检验较多或高度相关时 Bonferroni 往往保守。Holm 步进方法在一般依赖下也控制 FWER，通常更有力；这些方法的保证都对应一个确定的检验家族。

Benjamini–Hochberg 方法

Benjamini–Hochberg 方法让拒绝标准随发现数量变化。将 $m$ 个 $p$ 值排序为 $p_{(1)}\le \cdots \le p_{(m)}$，以 $q$ 表示目标 FDR 水平，最大合格排名 $k$ 满足

$$p_{(k)}\le \frac{k}{m}q.$$

最大合格排名确定后，前 $k$ 项共同构成拒绝集合。因而某项未通过自身排名阈值，并不一定被排除；后面更大的合格排名可能把它一起包含进来。这是作用于整个排序的步进规则。

例如六个 $p$ 值为 $0.002,0.011,0.018,0.041,0.20,0.67$，$q=0.05$：

排名	$p_{(i)}$	$iq/6$
1	0.002	0.00833
2	0.011	0.01667
3	0.018	0.02500
4	0.041	0.03333
5	0.200	0.04167
6	0.670	0.05000

表中最大合格排名为 $3$，所以 BH 拒绝前三项，而 Bonferroni 只拒绝第一项。差别来自控制目标：前者约束错误发现比例的期望，后者约束至少一次错误的概率，允许的错误结构并不相同。

标准 BH 的 FDR 保证建立在独立检验，或真零假设 $p$ 值满足适当正依赖条件（PRDS）等结构上，不能直接扩展到任意依赖。Benjamini–Yekutieli 等较保守方法及针对设计的方法处理不同依赖问题。FDR 为 $5\%$ 描述错误发现比例的期望，既不是每个发现为真的 $95\%$ 后验概率，也不是这次结果必然恰有 $5\%$ 错误。

检验家族对应一组共同的科学选择，例如同一研究的主要候选终点。删去不显著结果再减小 $m$，或为了增加显著结果而拆分家族，都改变了原控制所针对的错误事件；原来的成组保证不会自动保留。

抽样误差、模拟误差与成组判断

重抽样计算有两个误差来源：原研究只观察有限样本，计算又只生成有限次重复。增加模拟次数减少的是后一种 Monte Carlo 误差，观测单位、依赖结构和实际样本量仍决定前一种抽样误差。Bootstrap 与置换虽共享重复计算的形式，却分别依赖经验分布近似和零假设交换机制。

多重比较改变的则是判断事件的范围，FWER 与 FDR 体现了两种成组错误定义。列联表、效应量和回归模型的推断，也遵循相同的概率原则。理论补充见 参考教材与开放课程。

习题

1. 对 $10$ 个全零结果做非参数 bootstrap，会得到什么？为什么不能当作成功率为零的强证据？
2. 两独立组连续数据均值相等，但方差不同。直接交换标签并比较原始均值差，是否一定是精确检验？
3. 六项检验的 $p$ 值如上表。按 $\alpha =0.05$ 的 Bonferroni 和 $q=0.05$ 的 BH，分别拒绝哪些？解释其控制目标。
4. 四项排序 $p$ 值为 $0.009,0.026,0.027,0.20$，BH 的 $q=0.05$。第二项没达到自己的排名阈值，是否一定不能拒绝？

参考解答与推理

1. 每个重抽样仍是全零，成功比例和 percentile 区间均退化为零。经验分布不能生成未观察到的成功，故这反映方法在边界失效，应用二项 Wilson 或合适精确区间评估不确定性。
2. 不一定。均值相等不保证分布相同，方差不同时原始观测未必可交换。需要与弱零假设匹配的方法，如适当学生化且满足条件的渐近置换检验，或 Welch 方法。
3. Bonferroni 阈值 $0.05/6$，只拒绝 $0.002$；BH 最大合格排名为 $3$，拒绝前三项。前者控制至少一次错误的概率，后者在相应依赖条件下控制错误发现比例的期望。
4. 四个阈值为 $0.0125,0.025,0.0375,0.05$。第三项 $0.027\le 0.0375$，最大合格排名为 $3$，所以前三项都拒绝，包括第二项。BH 是向上寻找最大排名的整体规则。

上一篇：15 假设检验 · 下一篇：17 分类数据与列联表 · 课程总览

相关概念：交换性与重抽样 · p值与错误率。