Bonferroni와 Holm 보정¶

Bonferroni 보정 문턱은 \(\alpha/m = 0.05/5 = 0.01\)이다.

각 p-값을 0.01과 비교한다:

검정 1만 기각된다. Bonferroni는 보수적이다: 보정하지 않은 수준에서는 유의한 \(p = 0.025\)의 검정 2가 기각되지 않는다.

p-값을 정렬한다: \(p_{(1)} = 0.008, p_{(2)} = 0.025, p_{(3)} = 0.040, p_{(4)} = 0.060, p_{(5)} = 0.120\).

\(j\)번째 정렬된 검정의 Holm 문턱은 \(\alpha/(m - j + 1)\)이다:

2단계에서 기각하지 못하면 거기서 멈추고 남은 검정은 기각하지 않는다.

결과: 이 경우에는 Bonferroni와 마찬가지로 검정 1만 기각된다. 그러나 Holm은 Bonferroni보다 균일하게 더 강력하며(결코 더 적게 기각하지 않으며), p-값 배열이 다르면 더 많이 기각할 수 있다.

\(\mathcal{H}_0\)을 참인 귀무가설의 집합이라 하고 \(|\mathcal{H}_0| = m_0 \leq m\)이라 하자. FWER은:

첫 부등식은 Boole 부등식(합집합 한계)이다. 두 번째 등식은 참인 귀무가설 아래에서 p-값이 \((0,1)\) 위의 균등분포를 따르므로 \(P(p_i < \alpha/m) = \alpha/m\)이라는 사실을 쓴다. 마지막 부등식은 \(m_0 \leq m\)을 쓴다.

이 증명은 검정의 종속 구조에 대해 아무 가정도 하지 않는다. 그래서 Bonferroni가 임의의 종속성에서 타당하지만 보수적일 수 있는 것이다. \(\square\)

\(m_0 \ll m\)일 때: 증명에서 \(m_0 \leq m\)을 썼으므로 실제 FWER은 최대 \(m_0 \alpha/m\)이다. 1000개 귀무가설 중 10개만 참이면(\(m_0 = 10\)) 참 FWER은 최대 \(10 \times 0.05/1000 = 0.0005\)로 명목 \(\alpha = 0.05\)보다 훨씬 낮다. Bonferroni 문턱은 최악의 경우(\(m_0 = m\))에 맞추어져 있다.

양의 상관이 있을 때: 합집합 한계 \(P(\cup A_i) \leq \sum P(A_i)\)는 사건들이 서로 배반일 때만 등호에 가깝다. 검정이 양의 상관을 가지면(검정통계량이 자료를 공유할 때 흔하다) 사건 \(\{p_i < \alpha/m\}\)이 함께 일어나는 경향이 있어 합집합이 합보다 훨씬 작아진다. 실제 FWER이 \(\alpha\)보다 훨씬 작을 수 있다.

두 보수성의 원천 모두 검정력을 떨어뜨려 참 효과를 탐지하기 어렵게 만든다. Holm 절차가 첫 번째 문제를 부분적으로 다루고, 재표본추출 기반 방법(예: 순열 FWER)이 두 번째를 다룬다.

보정된 \(p\)-값. "이 값을 \(\alpha\)와 직접 비교하면 해당 절차의 결론이 나오는 수"다. 절차마다 정의가 다르다.

import numpy as np

p = np.array([0.001, 0.008, 0.039, 0.041, 0.23, 0.76])
m = len(p)
ps = np.sort(p)
i = np.arange(1, m + 1)

bonf = np.minimum(ps * m, 1)
sidak = 1 - (1 - ps)**m
holm = np.maximum.accumulate(np.minimum((m - i + 1) * ps, 1))
hoch = np.minimum.accumulate(((m - i + 1) * ps)[::-1])[::-1]
bh = np.minimum.accumulate((m / i * ps)[::-1])[::-1]
c_m = np.sum(1 / i)
by = np.minimum(np.minimum.accumulate((c_m * m / i * ps)[::-1])[::-1], 1)

print(f"{'p':>8s} {'본페로니':>9s} {'시닥':>8s} {'홀름':>8s} "
      f"{'호크버그':>9s} {'BH':>8s} {'BY':>8s}")
for j in range(m):
    print(f"{ps[j]:8.3f} {bonf[j]:9.4f} {sidak[j]:8.4f} {holm[j]:8.4f} "
          f"{hoch[j]:9.4f} {bh[j]:8.4f} {by[j]:8.4f}")
print("\nq = 0.10 에서 기각 수:",
      f"본페로니 {(bonf <= 0.10).sum()}, 홀름 {(holm <= 0.10).sum()}, "
      f"호크버그 {(hoch <= 0.10).sum()}, BH {(bh <= 0.10).sum()}, "
      f"BY {(by <= 0.10).sum()}")

       p      본페로니       시닥       홀름      호크버그       BH       BY
   0.001    0.0060   0.0060   0.0060    0.0060   0.0060   0.0147
   0.008    0.0480   0.0471   0.0400    0.0400   0.0240   0.0588
   0.039    0.2340   0.2123   0.1560    0.1230   0.0615   0.1507
   0.041    0.2460   0.2221   0.1560    0.1230   0.0615   0.1507
   0.230    1.0000   0.7916   0.4600    0.4600   0.2760   0.6762
   0.760    1.0000   0.9998   0.7600    0.7600   0.7600   1.0000

q = 0.10 에서 기각 수: 본페로니 2, 홀름 2, 호크버그 2, BH 4, BY 2

여섯 방법의 성격.

읽기.

1. 홀름이 본페로니를 지배한다. 모든 \(p\)에서 홀름의 보정값이 본페로니 이하다. 본페로니를 쓸 이유가 없다.

2. 호크버그가 홀름보다 강하다. 셋째·넷째에서 0.156 대 0.123이다. 다만 양의 종속 가정이 필요하다.

3. 시닥의 이득이 미미하다. 본페로니와 소수 셋째 자리에서 갈린다. 독립 가정을 추가로 지불할 가치가 없다.

4. BH가 훨씬 많이 기각한다. 4개 대 2개다. 통제 대상이 FWER이 아니라 FDR이기 때문이다. 공정한 비교가 아니다 — 목표가 다르다.

5. BY는 BH보다 훨씬 보수적이다. \(c(m)=\sum_{i=1}^m1/i\)배 만큼 문턱을 낮추며, \(m=6\)에서 \(c(6)=2.45\)다. \(m\)이 크면 \(c(m)\approx\log m\)으로 커져 대가가 크다.

선택 지침.

• FWER이 필요하면 홀름(가정 없음) 또는 호크버그(양의 종속이 그럴듯할 때).
• FDR로 충분하면 BH. 탐색적 대규모 검정의 표준이다.
• 종속 구조를 전혀 모르면 BY. 다만 대가가 크므로, 순열 기반 방법을 먼저 고려한다.
• 본페로니와 시닥은 교육용으로만.

홀름 절차. \(p_{(1)}\le\cdots\le p_{(m)}\)에 대해 \(k=1,2,\dots\) 순서로

를 확인하고, 처음 실패하는 곳에서 멈춘다. 그 앞까지 기각한다.

지배의 증명. 본페로니는 \(p_{(k)}\le\alpha/m\)를 요구한다. 홀름의 문턱 \(\alpha/(m-k+1)\)은 \(k\ge1\)에서

이므로 모든 단계에서 홀름의 문턱이 더 느슨하다. 따라서 본페로니가 기각하는 것은 홀름도 기각한다. \(\square\)

엄밀하게는 순차성을 확인해야 한다. 본페로니가 \(H_{(j)}\)를 기각하면 \(p_{(j)}\le\alpha/m\)이고, \(k\le j\)인 모든 \(k\)에 대해 \(p_{(k)}\le p_{(j)}\le\alpha/m\le\alpha/(m-k+1)\)이므로 홀름의 단계가 \(j\)까지 끊기지 않는다. 따라서 홀름도 기각한다.

FWER 통제의 증명 개요. 참인 귀무가설의 집합을 \(\mathcal I_0\)(\(|\mathcal I_0|=m_0\))라 하고, 그중 가장 작은 \(p\)-값을 \(p_{\min}\)이라 하자. 홀름이 참인 가설을 하나라도 기각했다면, 그 가설이 기각된 단계 \(k^*\)에서 \(p_{\min}\le\alpha/(m-k^*+1)\)이다. 그 단계 이전에 기각된 \(k^*-1\)개는 모두 거짓 귀무가설이므로 \(k^*-1\le m-m_0\), 즉 \(m-k^*+1\ge m_0\)이다. 따라서

합집합 한계로 \(P(p_{\min}\le\alpha/m_0)\le m_0\cdot\alpha/m_0=\alpha\)다. \(\square\)

그런데 왜 본페로니가 여전히 쓰이는가.

1. 단순함. "\(\alpha\)를 \(m\)으로 나눈다"는 한 줄로 설명된다. 홀름은 순차 절차라 설명이 길다.

2. 신뢰구간으로 곧장 확장된다. 본페로니는 \(1-\alpha/m\) 수준의 구간을 만들면 끝이다. 홀름에 대응하는 동시 신뢰구간은 만들기 어렵다 — 이것이 실무적으로 가장 큰 이유다.

3. 계획 단계에서 쓰기 쉽다. 표본크기 계산에 \(\alpha/m\)을 넣으면 된다. 홀름은 기각 순서에 의존해 사전 계산이 어렵다.

4. 차이가 작은 경우가 많다. 가장 작은 \(p\)-값에서는 두 절차가 동일하다. 주 결과 하나만 중요하다면 차이가 없다.

권고. 검정만 보고한다면 홀름을 쓴다. 이유 없이 검정력을 버릴 필요가 없다. 동시 신뢰구간이 필요하면 본페로니가 현실적인 선택이다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(19)
M, m = 40_000, 10

print(f"{'ρ':>5s} {'무보정 FWER':>12s} {'본페로니 FWER':>14s} "
      f"{'홀름 FWER':>11s} {'최대통계량 FWER':>16s}")
for rho in [0.0, 0.3, 0.6, 0.9, 0.99]:
    S = rho * np.ones((m, m)) + (1 - rho) * np.eye(m)
    L = np.linalg.cholesky(S)
    z = rng.standard_normal((M, m)) @ L.T
    p = 2 * stats.norm.sf(np.abs(z))
    ps = np.sort(p, axis=1)

    raw = (p < 0.05).any(1).mean()
    bonf = (p < 0.05 / m).any(1).mean()
    thr = 0.05 / (m - np.arange(m))
    holm = (ps <= thr).any(1).mean()       # 근사(첫 실패까지 보는 대신)

    # 최대통계량의 귀무분포를 몬테카를로로 구해 문턱을 잡는 방법
    zz = rng.standard_normal((M, m)) @ L.T
    cmax = np.percentile(np.abs(zz).max(1), 95)
    step = (np.abs(z).max(1) > cmax).mean()

    print(f"{rho:5.2f} {raw:12.4f} {bonf:14.4f} {holm:11.4f} {step:16.4f}")

    ρ     무보정 FWER      본페로니 FWER     홀름 FWER       최대통계량 FWER
 0.00       0.3985         0.0492      0.0492           0.0487
 0.30       0.3548         0.0451      0.0452           0.0501
 0.60       0.2462         0.0329      0.0330           0.0476
 0.90       0.1268         0.0169      0.0185           0.0467
 0.99       0.0693         0.0077      0.0352           0.0495

상관이 커질수록 본페로니가 보수적이 된다. \(\rho=0.9\)에서 0.017, \(\rho=0.99\)에서 0.008로 명목의 6분의 1이다.

왜 그런가. 본페로니는 합집합 한계 \(P(\cup A_i)\le\sum P(A_i)\)를 쓰는데, 사건들이 겹칠수록 이 한계가 느슨해진다. \(\rho=1\)이면 모든 검정이 같은 사건이라 실제 FWER이 \(\alpha/m=0.005\)가 된다. 위 표의 \(\rho=0.99\) 값 0.0077이 그 극한에 가깝다.

최대통계량 방법(순열/부트스트랩)이 정확하다. 모든 \(\rho\)에서 0.047~0.050이다. 상관 구조를 자료에서 반영하기 때문이다. 위 코드의 홀름 값은 첫 실패에서 멈추지 않은 근사라 \(\rho\)가 매우 클 때 조금 커지는데, 정식 홀름은 본페로니와 마찬가지로 보수적이 된다.

대안들.

유효 검정 수의 예. 고유값 기반 어림으로

를 쓰면, \(\rho=0.9\), \(m=10\)에서 \(m_{\text{eff}}=1.21\)이다. 사실상 하나의 검정임을 반영한다.

실무 권고. 유전체처럼 상관이 강한 대규모 검정에서는 순열 기반 방법이나 FDR을 쓴다. 본페로니를 그대로 쓰면 검정력을 크게 잃는다.

고정 순서 절차. 가설을 \(H_1,H_2,\dots,H_m\) 순서로 미리 정하고, 순서대로 각각 수준 \(\alpha\)로 검정하되 처음 기각에 실패하면 멈춘다.

왜 보정이 필요 없는가. \(H_1\)이 참이라면 그것을 기각할 확률이 \(\alpha\)이고, 그때만 오류가 시작된다. \(H_1\)이 거짓이고 \(H_2\)가 참이면, 오류는 \(H_1\)을 기각한 뒤에 \(H_2\)를 기각할 때 생기는데 그 확률도 \(\alpha\) 이하다. 일반적으로

핵심은 "멈춤"이다. 앞에서 막히면 뒤를 볼 수 없으므로, 여러 번 기회를 얻지 못한다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(27)
M, m = 40_000, 4

# H1, H2 는 거짓(효과 있음), H3, H4 는 참
shift = np.array([3.0, 2.5, 0.0, 0.0])
z = rng.standard_normal((M, m)) + shift
p = 2 * stats.norm.sf(np.abs(z))

# 고정 순서 절차
err = 0
for row in p:
    for j in range(m):
        if row[j] > 0.05:
            break
        if shift[j] == 0:          # 참인 가설을 기각
            err += 1
            break
print(f"고정 순서 절차의 FWER = {err / M:.4f}")
print(f"보정 없이 모두 검정할 때 FWER = "
      f"{(p[:, 2:] < 0.05).any(1).mean():.4f}")

고정 순서 절차의 FWER = 0.0314
보정 없이 모두 검정할 때 FWER = 0.0984

보정 없이도 FWER이 0.031로 유지된다. 모두 검정하면 0.098이다.

설계 시 주의점 다섯.

1. 순서를 자료를 보기 전에 정한다. 이것이 전부다. 결과를 보고 순서를 바꾸면 보장이 사라진다.

2. 중요한 것을 앞에. 앞에서 막히면 뒤를 아예 못 본다. 주 결과 → 핵심 부수 결과 → 그 밖 순서가 자연스럽다.

3. 앞의 가설이 기각될 가능성이 높아야 한다. 첫 가설의 검정력이 낮으면 전체가 막힌다. 검정력이 높은 것을 앞에 두는 전략과 중요도 순서가 충돌할 수 있다.

4. 분기 구조도 가능하다. 고정 순서가 너무 경직되면 그래프 기반 절차를 쓴다. 가설들을 노드로 두고, 기각 시 \(\alpha\)를 다른 노드로 "물려주는" 규칙을 미리 정한다. 유연하면서 FWER이 보장된다.

5. 부분적 실패의 처리. \(H_1\)은 기각, \(H_2\)는 실패, \(H_3\)는 \(p=0.001\)인 상황에서 \(H_3\)을 보고할 수 없다. 이것이 답답하지만 규칙을 지켜야 보장이 유지된다. 탐색적 결과로 따로 보고한다.

실무 사례. 임상시험에서 "주 평가변수 → 주요 부 평가변수 1 → 2 → 3" 순서를 계획서에 명시하는 것이 표준이 되었다. 규제기관이 이 구조를 권장한다.

핵심 질문 — "무엇이 하나의 가족인가." 보정의 범위를 정하는 것이 가장 어렵고 가장 임의적인 부분이다.

보정해야 하는 경우.

1. 하나의 결론을 여러 검정으로 뒷받침할 때. "어느 하나라도 유의하면 효과가 있다"는 논리를 쓴다면 FWER 통제가 필요하다.

2. 여러 결과변수·시점·하위집단을 훑고 유의한 것을 보고할 때. 앞서 본 HARKing의 구조다.

3. 규제 제출. 확증적 시험에서 여러 주 평가변수나 용량군이 있으면 보정이 요구된다.

4. 대규모 탐색. 유전체, 뇌영상, 대규모 A/B 시험. 다만 FWER보다 FDR이 적절한 경우가 많다.

보정하지 말아야 하는(또는 불필요한) 경우.

1. 각 검정이 독립적인 연구 질문일 때. 한 논문에 서로 무관한 세 연구가 실렸다고 보정하지 않는다.

2. 모든 결과를 빠짐없이 보고하고, 독자가 스스로 판단할 수 있을 때. 로드먼·펄먼 등이 주장하는 입장이다.

3. 주 결과가 하나로 사전에 정해져 있고 나머지는 탐색적이라고 명시할 때.

4. 추정이 목적일 때. 여러 모수의 크기를 보고하는 것이 목적이면, 각각의 구간을 보고하고 "동시 보장은 없다"고 밝히면 된다.

애매한 회색지대.

널리 받아들여지는 원칙 셋.

1. "하나의 결론"이 기준이다. 여러 검정이 같은 주장을 뒷받침하면 한 가족이다.

2. 투명성이 보정을 대체할 수 있다. 몇 개의 검정을 했고 무엇이 유의했는지 전부 보고하면, 보정하지 않아도 독자가 할인해 읽을 수 있다. 숨기는 것이 진짜 문제다.

3. 사전등록이 가족을 정의한다. 계획서에 "주 결과 하나, 부수 결과 다섯, 보정 방법 X"라고 적으면 논쟁이 사라진다.

베이즈 관점. 계층모형으로 여러 모수를 함께 추정하면 축소(shrinkage)가 다중성을 자동으로 처리한다. 별도의 보정이 필요 없다는 것이 겔먼 등의 주장이다. 다만 모형 가정이 필요하다.

보고할 것 일곱.

1. 검정의 총 개수 \(m\). 보정의 분모다.
2. 가족의 정의. 무엇을 한 가족으로 묶었고 왜인지.
3. 통제 대상. FWER인가 FDR인가. 수준은 얼마인가.
4. 방법 이름. 본페로니/홀름/BH/순열 기반 등. 순열이면 횟수도.
5. 원 \(p\)-값과 보정된 \(p\)-값을 모두. 독자가 다른 보정을 적용해 볼 수 있어야 한다.
6. 효과크기와 구간. 보정된 동시 구간이면 그렇다고 밝힌다.
7. 사전 계획 여부. 보정 방법을 미리 정했는가.

좋은 표의 형태.

이 표가 좋은 이유. 원 \(p\)와 보정 \(p\)를 나란히 두어 보정의 영향이 한눈에 보이고, 효과크기와 구간이 있어 실무적 판단이 가능하다.

흔한 잘못.

• 보정 후 결과만 싣기. 독자가 다른 보정이나 메타분석에 쓸 수 없다.
• \(m\)을 밝히지 않기. 보정을 재현할 수 없다.
• 유의한 것만 표에 넣기. 선택적 보고다.
• "BH 보정했다"만 적고 \(q\) 수준을 안 밝히기.
• 구간은 보정하지 않고 \(p\)만 보정하기. 앞서 본 쌍대성의 불일치다.
• 사후에 보정 방법을 바꾸기. 유의해지는 방법을 고르면 보정의 의미가 사라진다.

FDR을 쓸 때의 추가 사항.

• \(q\)-값을 보고한다. "이 가설을 기각할 때 예상되는 FDR"로 해석된다.
• 기각된 가설 수와 예상 거짓발견 수를 함께 적는다. "\(q=0.05\)에서 200개를 기각했으며, 이 중 약 10개가 거짓발견으로 예상된다."
• 독립성 가정이 그럴듯한지 논한다. 아니면 BY나 순열 기반을 쓴다.

한 문장. 보정은 결론을 바꾸는 절차이므로, 그 절차 자체가 검증 가능해야 한다. 원 \(p\)-값, \(m\), 방법을 모두 적으면 그것이 가능해진다.