多重比較與次族群分析

0.05 承諾的是一次，不是二十次

α = 0.05 的意思是：在沒有效果的情況下，這一次檢定有 5% 的機率會宣稱有效果。 它對「這一次」以外的任何事都沒有承諾。

如果二十個結果彼此獨立、而且全部都真的沒有效果，至少出現一個「顯著」的機率是

這不是模擬，是算術。但下面那個是模擬：5,000 次虛構試驗， 每次每組 60 人、20 個結果指標， 而且每一個指標的兩組在程式碼裡都是從同一個分布抽出來的——真實效果全部是零， 所以每一個「顯著」都是假陽性，不需要推論。

• 平均每次試驗跑出 0.995 個「顯著」結果（理論值正好是 20 × 0.05 = 1）
• 64.3% 的試驗至少有一個
• 最多的那一次跑出 5 個

三種調整方法，各自保證的東西不一樣

上表三種調整的族系錯誤率都落在 5.06% 附近， 看起來一模一樣——那是因為這個模擬裡真實效果全部是零。 在「全部虛無」的情況下偽發現率控制與族系錯誤率控制會重合； 一旦有些效果是真的，Benjamini–Hochberg 會拒絕得比 Bonferroni 多，代價是容許被宣稱的陽性裡有一定比例是錯的。

體學論文（基因體、蛋白質體、代謝體、PheWAS）裡那張火山圖（volcano plot）就是上面這一節換一個畫法：橫軸放效果量、縱軸放 $-\log_{10}(p)$，於是 p 值愈小的點愈高、效果量愈極端的點愈往兩側，兩條決策線就是兩條水平線。下面這張用的不是新資料，是本頁同一個模擬裡的第 2 次試驗（挑選規則事先定好：第一個假陽性個數剛好等於期望值 20 × 0.05 = 1 的那一次，不是挑最好看的一次），20 個結果指標的真實效果全部是零。未調整的 0.05 那條線之上有 1 個點，Benjamini–Hochberg 的線之上有 0 個。站在紅線之上那一點不是訊號，是這頁從頭到尾在講的那種假陽性——它的原始 p 是 0.029，BH 調整後變成 0.57。換句話說，火山圖上「有沒有點浮在線上面」是一個結論，不是資料的性質：同一批點換一條線就換一個結論，而換線不需要重跑任何分析。看這種圖的第一件事因此是問那條線是哪一種門檻。

那條綠線落在哪裡還需要多說一句，因為它不是一個固定高度。 Benjamini–Hochberg 是逐級上行（step-up）的程序：把 20 個 p 值由小排到大，第 i 小的那個要跟自己那一階比，階高是 α × i / 20。能通過自己那一階的最大 p 值就是那條線，它以下的全部宣告為陽性——所以線的高度跟著資料走，真的有訊號時它可以遠高於 Bonferroni 的門檻。Bonferroni 則相反，它讓每一個 p 值都只跟最低那一階比。而這次試驗裡沒有任何 p 值通過任何一階，這樣的 p 值並不存在，能誠實畫出來的只有最低那一階本身，α / 20 = 0.0025——數值上就等於 Bonferroni 的門檻。這不與上表矛盾：沒有人爬上梯子的時候，兩個程序在梯腳相遇。

次族群分析：問題不在有沒有調整，在問錯了問題

臨床試驗最常見的多重比較不是二十個結果指標，是一張次族群森林圖。 下面的模擬只有一個結果、一個真實效果，而且每一個次族群的真實效果完全相同—— 沒有任何交互作用存在。試驗總人數固定，切成 g 個等分。

中間那一欄是「這個藥對某一群人有效、對另一群無效」這種句子的來源，而它的可信度取決於切了幾組。 切成 6 組時，有 96.7% 的試驗 會出現至少一組顯著、至少一組不顯著的分裂結果——而真相是每一組的效果一模一樣。 單純的男女二分是另一回事：只切 2 組時分裂率是 11.3%， 下面的 Callout 會回到這個對比。

右欄是對照組：正確的問法不是「哪一組顯著」，而是「各組的效果是不是真的不同」， 也就是交互作用檢定。它的假陽性率一路維持在 5.2% 到 7.0% 之間——貼著名目的 5%， 上緣略高（組數多時每組只剩幾十人，大樣本近似開始失準，見本節末）。

兩個誠實的但書，因為它們也在同一張表上：

• 中間那一欄在 8 組時是 98.7%， 到 10 組反而略降到 97.5%。 它不是單調的：切得夠細之後每一組都太小，於是「全部都不顯著」開始變常見， 而全部不顯著就不算分裂。分裂率下降不代表問題變小，代表資訊已經稀釋到什麼都看不出來。
• 右欄在 10 組時是 7.0%， 比名目的 5% 高一些。這是大樣本近似在每組只剩幾十人時的誤差， 不是交互作用檢定本身失效——但它提醒了一件事：組數多到一定程度，連正確的檢定也開始不可靠。

什麼時候該調整，什麼時候不該

讀論文時看哪裡

1. 計畫書或註冊資料裡的主要結果是哪一個。 摘要裡的那一句如果不是它，先問為什麼。
2. 總共量了幾個結果。 方法段落列了十五個、結果段落報了三個，另外十二個在哪裡。
3. 次族群是事先指定還是事後的。 CONSORT 要求講清楚；沒講清楚就當作事後的。
4. 有沒有報交互作用的 p 值。 只報各組的效果與各組的 p 值，等於在請讀者做那個會出錯的比較。
5. 調整方法有沒有寫，以及控制的是哪一種錯誤率。 「已用 Bonferroni 校正」與「FDR < 0.05」是兩個不同的宣稱。
6. 期中分析的 alpha 花掉多少。 最終的顯著門檻通常小於 0.05。

動手跑一次

set.seed(20260822)

# One trial: 20 endpoints, no effect anywhere
one_trial <- function(k = 20, n = 60) {
  replicate(k, t.test(rnorm(n), rnorm(n), var.equal = TRUE)$p.value)
}

p <- one_trial()
sum(p < 0.05)                          # how many "significant" findings from nothing

p.adjust(p, "bonferroni") < 0.05       # all FALSE, almost always
p.adjust(p, "holm")       < 0.05
p.adjust(p, "BH")         < 0.05

# The family-wise error rate, measured rather than asserted
mean(replicate(1000, any(one_trial() < 0.05)))                       # about 0.64
mean(replicate(1000, any(p.adjust(one_trial(), "holm") < 0.05)))     # about 0.05

import numpy as np
from scipy import stats
from statsmodels.stats.multitest import multipletests

rng = np.random.default_rng(20260822)

def one_trial(k=20, n=60):
    return np.array([stats.ttest_ind(rng.normal(size=n), rng.normal(size=n)).pvalue
                     for _ in range(k)])

p = one_trial()
print((p < 0.05).sum())

for method in ("bonferroni", "holm", "fdr_bh"):
    print(method, multipletests(p, alpha=0.05, method=method)[0].sum())

print(np.mean([(one_trial() < 0.05).any() for _ in range(1000)]))     # about 0.64

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B8-05-multiplicity.R

回報這一頁的問題