구간 추정·가설 검정·부트스트랩

신뢰구간

추정값 하나만 보고하면 독자는 그것이 얼마나 정확한지 알 수 없다. 표준오차를 함께 쓰는 것도 방법이지만, 더 직접적인 답은 "참값이 있을 법한 범위"다.

정의. 모수 $\theta$에 대한 신뢰수준 $1-\alpha$의 신뢰구간(confidence interval)은 표본으로 계산되는 구간 $[L(X), U(X)]$로서, 모든 $\theta$에 대해

를 만족하는 것이다. 확률의 주어는 구간이다. $\theta$는 고정된 미지량이고, 무작위인 것은 표본에서 계산된 $L$과 $U$다. 그래서 "95% 신뢰구간"은 이렇게 읽는다. 같은 절차로 실험과 구간 계산을 반복하면, 그렇게 만든 구간들의 95%가 참값을 덮는다. 이미 계산된 구간 하나는 참값을 덮거나 덮지 않거나 둘 중 하나다.

구성. 정규분포 표본의 평균에 대해서는 다음 통계량이 $\mu$와 $\sigma$에 상관없이 자유도 $n-1$인 $t$ 분포를 따른다.

분포가 미지량에 의존하지 않는 이런 양을 pivot이라 한다. $t_{n-1}$의 $97.5\%$ 분위수 $t^\star$로 $\Pr[-t^\star \le T \le t^\star] = 0.95$를 $\mu$에 대해 풀면 구간 $\bar{X} \pm t^\star S/\sqrt{n}$을 얻는다. 데이터가 정규분포가 아니어도 $n$이 크면 중심극한정리(11.5장)로 근사적으로 성립한다. $n$이 클 때 $t^\star \to 1.96$이므로 "추정값 ± 1.96 × 표준오차"가 된다.

"근사적으로"가 얼마나 근사적인지는 시뮬레이션으로 확인할 수 있다. 같은 실험을 20,000번 반복해 구간이 참 평균을 덮는 비율(포함 확률, coverage)을 잰다. 정규 데이터와 오른쪽 꼬리가 긴 로그정규 데이터를 비교한다.

import numpy as np
from math import erf, sqrt

rng = np.random.default_rng(0)
t975 = {9: 2.262, 29: 2.045}           # t 분포의 97.5% 분위수 (자유도 n-1)

def coverage(sample, mean_true, n, R=20000):
    X = sample((R, n))                                    # 같은 실험을 R번 반복, (R, n)
    m, se = X.mean(axis=1), X.std(axis=1, ddof=1) / np.sqrt(n)
    half = t975[n - 1] * se                               # 구간 반폭, (R,)
    return np.mean((m - half <= mean_true) & (mean_true <= m + half))

for n in (10, 30):
    c_norm = coverage(lambda s: rng.normal(5.0, 2.0, s), 5.0, n)
    c_skew = coverage(lambda s: rng.lognormal(0.0, 1.0, s), np.exp(0.5), n)
    print(f"n = {n}:  정규 데이터 {c_norm:.3f}   로그정규 데이터 {c_skew:.3f}   (목표 0.950)")

n = 10:  정규 데이터 0.951   로그정규 데이터 0.844   (목표 0.950)
n = 30:  정규 데이터 0.951   로그정규 데이터 0.880   (목표 0.950)

정규 데이터에서는 $n = 10$에서도 포함 확률이 명목값 0.95와 맞는다. 로그정규 데이터에서는 $n = 30$에서도 0.88에 그친다. "95% 구간"이라고 써 놓고 실제로는 8번 중 1번꼴로 참값을 놓치는 것이다. 꼬리가 긴 데이터에서는 큰 값이 드물게 나와 표본 평균과 표본 표준편차가 함께 작게 나오는 경우가 많기 때문이다. 이때는 로그 변환 뒤 구간을 만들거나, 뒤에서 다룰 부트스트랩 BCa 구간을 쓴다.

가설 검정

신뢰구간이 "어디쯤인가"를 묻는다면, 가설 검정은 "이 차이가 우연만으로 설명되는가"를 묻는다. 새 촉매가 수율을 높였는지, 새 모델이 오차를 줄였는지처럼 예·아니오로 답해야 하는 판단에 쓴다.

구성 요소.

• 귀무가설(null hypothesis) $H_0$: "효과가 없다"는 기준 가설. 예: 두 그룹의 평균이 같다.
• 대립가설 $H_1$: $H_0$이 아닐 때. 예: 평균이 다르다.
• 검정통계량 $T$: $H_0$에서 벗어날수록 커지도록 만든 통계량. 예: 평균 차이를 표준오차로 나눈 값.
• p-value: $H_0$이 참이라고 가정할 때, 관측된 값만큼 또는 그보다 극단적인 $T$가 나올 확률.
• 유의수준 $\alpha$: p-value가 이보다 작으면 $H_0$을 기각한다고 미리 정한 기준. 관례적으로 0.05를 쓰지만 근거가 있는 값은 아니다.

p-value는 $H_0$이 참일 확률이 아니다. $\Pr[\text{데이터} \mid H_0]$를 $\Pr[H_0 \mid \text{데이터}]$로 뒤집으려면 $H_0$의 사전확률이 필요하다. 미국통계학회(ASA)는 p-value의 오용이 널리 퍼진 상황에서, p-value가 가설이 참일 확률을 재지 않고 효과의 크기나 중요성을 재지도 않는다는 원칙을 공식 성명으로 정리했다(Wasserstein & Lazar, 2016).

두 가지 오류. 판단이 틀리는 방식은 두 가지다.

$H_0$이 참이면 p-value는 $[0, 1]$에서 균등분포를 따르므로, $\alpha = 0.05$로 검정하면 정확히 5%의 확률로 1종 오류가 난다. 검정력은 참 효과의 크기, 표본 크기, 잡음 크기에 따라 정해진다. 표준편차 $\sigma$를 아는 두 그룹 비교(그룹당 $n$개)에서 효과 크기를 $d = (\mu_1 - \mu_0)/\sigma$라 하면

이다. $\Phi$는 표준정규 누적분포함수, $z_q$는 그 $q$ 분위수다. $d = 0.5$, $\alpha = 0.05$에서 검정력 80%($z_{0.8} \approx 0.84$)를 원하면 $n \approx 2(1.96 + 0.84)^2/0.25 \approx 63$, 즉 그룹당 60여 개가 필요하다. 시뮬레이션으로 확인한다.

Phi = np.vectorize(lambda z: 0.5 * (1 + erf(z / sqrt(2))))   # 표준정규 CDF

def z_test_p(a, b):                      # σ = 1을 아는 두 표본 z 검정, 양측 p-value
    z = (b.mean(axis=1) - a.mean(axis=1)) / np.sqrt(2 / a.shape[1])
    return 2 * (1 - Phi(np.abs(z)))

R = 20000
for d, n in ((0.0, 20), (0.5, 20), (0.5, 64)):           # d: 참 효과 크기 (σ 단위)
    p = z_test_p(rng.normal(0, 1, (R, n)), rng.normal(d, 1, (R, n)))
    theory = Phi(d * np.sqrt(n / 2) - 1.96) + Phi(-d * np.sqrt(n / 2) - 1.96)
    print(f"d = {d}, n = {n:2d}:  P(p < 0.05) = {np.mean(p < 0.05):.3f}  (이론 {theory:.3f})"
          f"   P(p < 0.01) = {np.mean(p < 0.01):.3f}")

d = 0.0, n = 20:  P(p < 0.05) = 0.052  (이론 0.050)   P(p < 0.01) = 0.010
d = 0.5, n = 20:  P(p < 0.05) = 0.354  (이론 0.353)   P(p < 0.01) = 0.160
d = 0.5, n = 64:  P(p < 0.05) = 0.810  (이론 0.807)   P(p < 0.01) = 0.599

첫 줄은 $H_0$이 참일 때 p-value가 0.05 아래로 떨어지는 비율이 5%, 0.01 아래가 1%라는 것을 보여 준다. p-value가 균등분포라는 뜻이다. 둘째 줄이 더 중요하다. 실제로 중간 크기의 효과($d = 0.5$)가 있어도 그룹당 20개로는 세 번 중 두 번 유의하지 않게 나온다. 검정력이 낮은 실험에서 "유의하지 않다"는 결과는 "효과가 없다"는 증거가 거의 되지 못한다.

신뢰구간과 검정의 관계. 양측 검정에서 $p < \alpha$인 것과 $1-\alpha$ 신뢰구간이 귀무가설의 값(예: 차이 0)을 포함하지 않는 것은 같은 판단이다. 신뢰구간은 같은 정보에 더해 효과가 어느 크기인지까지 보여 주므로, 검정 결과만 쓰기보다 구간을 함께 쓰는 것이 낫다.

효과 크기

p-value는 "차이가 0인가"에 답할 뿐 "차이가 얼마나 큰가"에는 답하지 않는다. 그리고 p-value는 표본 크기의 함수다. 같은 효과도 $n$을 키우면 p-value가 얼마든지 작아진다. 그래서 결과를 보고할 때는 효과 크기(effect size)와 그 신뢰구간을 함께 쓴다.

• 원래 단위의 차이: 평균 차이 $\bar{x}_1 - \bar{x}_0$와 그 신뢰구간. 도메인의 기준(공정 허용오차, 측정 장비의 분해능, 임상적으로 의미 있는 차이)과 바로 비교할 수 있어 가장 유용하다.
• 표준화된 차이: Cohen's $d = (\bar{x}_1 - \bar{x}_0)/s_{\text{pooled}}$. 단위가 없어 서로 다른 연구를 비교할 때 쓴다. 차이가 측정값 퍼짐의 몇 배인지를 말한다.

아주 큰 표본에서의 작은 효과와 작은 표본에서의 중간 효과를 비교한다.

for d, n in ((0.02, 100000), (0.5, 20)):
    a, b = rng.normal(0, 1, n), rng.normal(d, 1, n)        # 그룹당 n개
    diff = b.mean() - a.mean()
    se = np.sqrt(a.var(ddof=1) / n + b.var(ddof=1) / n)
    pooled_sd = np.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
    p = 2 * (1 - Phi(abs(diff / se)))
    print(f"n = {n:6d}:  차이 {diff:+.3f}  95% CI [{diff - 1.96 * se:+.3f}, {diff + 1.96 * se:+.3f}]"
          f"  p = {p:.1e}  Cohen's d = {diff / pooled_sd:.2f}")

n = 100000:  차이 +0.024  95% CI [+0.015, +0.033]  p = 1.0e-07  Cohen's d = 0.02
n =     20:  차이 +0.939  95% CI [+0.296, +1.581]  p = 4.2e-03  Cohen's d = 0.90

첫 줄은 p-value가 매우 작지만 차이는 측정값 퍼짐의 2%에 불과하다. 신뢰구간이 좁으므로 "효과가 작다"는 것까지 정밀하게 알 수 있다. 둘째 줄도 유의하지만 신뢰구간이 0.3에서 1.6까지 걸쳐 있어 효과의 크기를 거의 모른다. 참 효과는 0.5인데 추정값은 0.94로 두 배 가까이 나왔다. 검정력이 낮은 실험에서는 우연히 크게 나온 경우만 유의 기준을 넘으므로, 유의하게 나온 효과 추정값은 참값보다 체계적으로 크기 쉽다. 앞 절의 시뮬레이션에서 $n = 20$일 때 유의한 결과는 35%뿐이었고, 그 35%는 효과가 크게 나온 표본들이다.

부트스트랩

신뢰구간을 만들려면 추정량의 표본분포를 알아야 한다. 평균에는 $t$ 분포가 있지만, 중앙값·분위수·상관계수·두 모델의 AUC 차이·회귀계수의 비 같은 통계량에는 손으로 유도한 공식이 없거나 가정이 너무 강하다. 부트스트랩(bootstrap)은 이 표본분포를 계산으로 근사한다(Efron, 1979).

아이디어. 표본분포는 "모집단에서 표본을 반복해서 뽑으면 통계량이 어떻게 흩어지는가"다(12.1장). 모집단은 모르지만 표본은 모집단을 닮았다. 그래서 표본을 모집단 대신 쓴다. 크기 $n$인 표본에서 $n$개를 복원추출(같은 관측이 여러 번 뽑힐 수 있게)해 새 표본을 만들고, 통계량을 계산하는 일을 $B$번 반복한다. 이 $B$개 값의 분포가 표본분포의 근사다.

• 부트스트랩 표준오차: $B$개 값의 표준편차.
• percentile 구간: $B$개 값의 $2.5\%$, $97.5\%$ 백분위수.

배터리 셀 40개의 사이클 수명(오른쪽으로 치우친 분포)에서 중앙값의 표준오차와 95% 구간을 구한다. 시뮬레이션이므로 모집단에서 새 표본을 반복해 뽑아 "참" 표준오차와도 비교할 수 있다.

life = rng.lognormal(np.log(800), 0.4, size=40)           # 셀 40개의 사이클 수명, (n,)
B = 5000
idx = rng.integers(0, len(life), size=(B, len(life)))     # 복원추출 인덱스, (B, n)
boot = np.median(life[idx], axis=1)                       # 부트스트랩 중앙값, (B,)
lo, hi = np.percentile(boot, [2.5, 97.5])
print(f"표본 중앙값 {np.median(life):.0f} 사이클,  부트스트랩 SE {boot.std():.0f},  95% percentile 구간 [{lo:.0f}, {hi:.0f}]")

pop = np.median(rng.lognormal(np.log(800), 0.4, size=(B, 40)), axis=1)   # 모집단에서 새로 뽑은 표본들
print(f"모집단에서 반복 추출한 중앙값의 표준편차 {pop.std():.0f}  (부트스트랩이 흉내 내려는 값)")

표본 중앙값 815 사이클,  부트스트랩 SE 69,  95% percentile 구간 [740, 971]
모집단에서 반복 추출한 중앙값의 표준편차 63  (부트스트랩이 흉내 내려는 값)

표본 하나만으로 얻은 부트스트랩 표준오차(69)가 모집단을 알아야 계산할 수 있는 값(63)에 가깝다. 구간은 중앙값 815를 중심으로 비대칭이다(아래로 75, 위로 156). 분포의 오른쪽 꼬리를 반영한 것으로, "± 1.96 × SE" 꼴의 대칭 구간에서는 얻을 수 없는 정보다.

BCa 구간. percentile 구간은 추정량에 편향이 있거나 표본분포의 퍼짐이 모수에 따라 달라지면 포함 확률이 명목값에서 벗어난다. Efron(1987)의 BCa(bias-corrected and accelerated) 구간은 부트스트랩 분포의 중앙값이 원래 추정값에서 얼마나 벗어났는지(bias correction)와 퍼짐이 얼마나 빨리 변하는지(acceleration, jackknife로 추정)를 이용해 잘라 낼 백분위수를 조정한다. 많은 문제에서 percentile 구간보다 포함 확률이 명목값에 가깝다. 직접 구현하기보다 라이브러리를 쓰는 것이 안전하다. 집필 시점 예시로, SciPy의 scipy.stats.bootstrap은 method 인자로 percentile·basic·BCa 구간을 고를 수 있다. 사용하는 버전의 문서에서 인자 이름과 기본값을 확인한다.

그룹 부트스트랩. 부트스트랩은 분포의 모양을 가정하지 않지만 관측이 독립이라는 가정은 그대로 한다. 같은 배터리 셀의 사이클들, 같은 환자의 영상들, 같은 기상 관측소의 날짜들은 서로 상관되어 있다. 이런 데이터에서 개별 관측을 재표본화하면 독립 정보의 양을 과대평가해 구간이 지나치게 좁아진다. 해법은 독립인 단위(셀·환자·관측소)를 재표본화하고, 뽑힌 단위의 관측을 통째로 쓰는 것이다. 이를 그룹(또는 cluster, block) 부트스트랩이라 한다. 무엇이 독립 단위인지는 결론을 어디에 일반화하려는지가 정한다. "새 셀에서도 이 모델이 낫다"고 말하려면 셀이 단위다(실습).

부트스트랩이 실패하는 곳. 표본이 아주 작으면(예: 10개 미만) 표본이 모집단을 잘 닮지 않는다. 최댓값·최솟값처럼 표본의 극단에 의존하는 통계량은 복원추출로 표본 밖의 값을 만들 수 없어 표본분포를 흉내 내지 못한다. 시계열처럼 의존 구조가 있으면 연속된 구간을 묶어 뽑는 block bootstrap이 필요하다.

다중 비교와 선택 편향

검정 하나의 1종 오류 확률이 5%라는 것은 그 검정 하나만 했을 때의 이야기다. 여러 검정을 하고 그중 유의한 것을 보고하면, 아무 효과가 없어도 무엇인가는 유의하게 나온다. 독립인 검정 $m$개를 모두 $\alpha$로 하면 적어도 하나가 거짓 양성일 확률, 즉 가족별 오류율(family-wise error rate, FWER)은

이다. $m = 10$이면 40%, $m = 100$이면 99.4%다.

대응은 무엇을 제어할지에 따라 둘로 나뉜다.

• Bonferroni: 각 검정을 $\alpha/m$으로 한다. FWER을 $\alpha$ 이하로 제어하며, 검정들 사이의 의존성과 무관하게 성립한다. 대신 $m$이 크면 매우 보수적이어서 참 효과도 놓친다.
• Benjamini–Hochberg(BH): p-value를 정렬해 $p_{(k)} \le \frac{k}{m}q$를 만족하는 가장 큰 $k$까지 기각한다. 기각한 것들 중 거짓 발견의 비율의 기댓값, 즉 FDR(false discovery rate)을 $q$ 이하로 제어한다(Benjamini & Hochberg, 1995). 원 논문의 증명은 검정들이 독립이라고 가정한다. 수천 개 유전자나 화합물을 선별하는 연구처럼 "후보 목록에 거짓이 조금 섞여도 되지만 그 비율은 알고 싶다"는 상황에 맞는다.

검정 100개 중 실제 효과가 있는 것이 10개인 상황을 2,000번 반복해 세 방법을 비교한다.

m, m1, R = 100, 10, 2000                       # 검정 100개, 그중 실제 효과 10개, 반복 R번
mu = np.r_[np.full(m1, 3.0), np.zeros(m - m1)]
p = 2 * (1 - Phi(np.abs(rng.normal(mu, 1.0, size=(R, m)))))   # (R, m)

def bh(p, q=0.05):                             # 행마다 Benjamini–Hochberg 기각 여부
    order = np.argsort(p, axis=1)
    ok = np.take_along_axis(p, order, axis=1) <= q * np.arange(1, p.shape[1] + 1) / p.shape[1]
    k = np.where(ok.any(axis=1), p.shape[1] - np.argmax(ok[:, ::-1], axis=1), 0)   # 통과한 최대 순위
    rej = np.zeros_like(p, dtype=bool)
    np.put_along_axis(rej, order, np.arange(p.shape[1]) < k[:, None], axis=1)
    return rej

for name, rej in (("보정 없음", p < 0.05), ("Bonferroni", p < 0.05 / m), ("BH (q=0.05)", bh(p))):
    false = rej[:, m1:].sum(axis=1)
    fdp = false / np.maximum(rej.sum(axis=1), 1)
    print(f"{name:11s}  발견 {rej.sum(axis=1).mean():5.2f}  참 발견 {rej[:, :m1].sum(axis=1).mean():4.2f}"
          f"  P(거짓 발견 ≥ 1) {np.mean(false > 0):.3f}  평균 FDP {fdp.mean():.3f}")

보정 없음        발견 12.97  참 발견 8.51  P(거짓 발견 ≥ 1) 0.989  평균 FDP 0.330
Bonferroni   발견  3.21  참 발견 3.15  P(거짓 발견 ≥ 1) 0.052  평균 FDP 0.014
BH (q=0.05)  발견  5.37  참 발견 5.08  P(거짓 발견 ≥ 1) 0.238  평균 FDP 0.043

보정하지 않으면 거의 매번 거짓 발견이 섞이고, 발견 목록의 3분의 1이 거짓이다. Bonferroni는 거짓 발견이 하나라도 나올 확률을 5% 근처로 묶지만 참 효과 10개 중 3개 정도만 찾는다. BH는 거짓 발견 비율을 평균 4.3%로 유지하면서 참 효과를 더 많이 찾는다. 대신 거짓 발견이 하나 이상 섞인 반복이 24%다. 어느 쪽이 맞는지는 거짓 발견 하나의 비용이 정한다.

보이지 않는 다중 비교. 더 위험한 것은 검정의 수가 드러나지 않는 경우다. 연구자는 이상치 제거 기준, 공변량, 데이터 부분집합, 변환, metric, 하이퍼파라미터를 데이터를 보면서 고를 수 있다. 각 선택이 합리적이어도, 데이터에 따라 분석이 달라질 수 있었다는 사실 자체가 다중 비교와 같은 효과를 낸다. 연구자가 의도적으로 여러 분석을 돌려 보지 않았더라도 그렇다. Gelman과 Loken(2014)은 이를 갈림길의 정원(garden of forking paths)이라 불렀다. 의도적으로 유의한 결과가 나올 때까지 분석을 바꾸는 것은 p-hacking이다.

ML 연구에서 같은 문제가 다른 이름으로 나타난다. 테스트 세트를 보면서 모델·seed·하이퍼파라미터를 고르면 테스트 점수는 그 선택에 맞춰진다(15.4장의 누수, 16.6장의 모델 선택과 최종 보고의 분리). seed 10개 중 가장 좋은 결과만 보고하는 것은 검정 10개 중 가장 작은 p-value만 보고하는 것과 같다(12.5장).

대응책은 분석 경로를 데이터와 분리하는 것이다.

• 사전 등록: 가설·분석 방법·주 지표를 데이터를 보기 전에 기록한다. 탐색적 분석은 탐색적이라고 밝힌다.
• holdout: 탐색에 쓰지 않은 데이터(테스트 세트, 새 실험)에서 최종 확인을 한 번만 한다.
• 모두 보고: 시도한 분석·seed·설정의 수와 전체 결과를 함께 보고한다.
• 보정: 검정의 수가 명확하면 Bonferroni·BH로 보정한다.

참고 문헌

• Efron (1979). Bootstrap methods: Another look at the jackknife. The Annals of Statistics, 7(1), 1–26. — 부트스트랩의 원 논문
• Efron (1987). Better bootstrap confidence intervals. Journal of the American Statistical Association, 82(397). — BCa 구간
• Benjamini, Hochberg (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B, 57(1), 289–300. — FDR의 정의와 BH 절차
• Gelman, Loken (2014). The statistical crisis in science. American Scientist, 102(6), 460. — 갈림길의 정원
• Wasserstein, Lazar (2016). The ASA's statement on p-values: Context, process, and purpose. The American Statistician, 70(2), 129–133. — p-value 해석에 대한 미국통계학회의 원칙