두 집단 비교를 위한 표본크기¶

임계값은 \(z_{0.025} = 1.96\), \(z_{0.20} = 0.842\)이다(검정력 \(= 1 - \beta = 0.80\)이므로 \(\beta = 0.20\)).

올림하면 집단당 \(n = 63\)명(총 126명)이다.

올림하면 집단당 \(n = 63\)명이다. 널리 인용되는 기준값이다: 검정력 80%로 중간 크기의 효과를 탐지하려면 집단당 약 63명이 필요하다.

공식을 정리하면:

집단당 30명이면 검정력 80%로 약 7.23 단위의 차이(표준화 효과로는 \(d = 7.23/10 = 0.72\))를 탐지할 수 있다.

공식 \(n = 2(z_{\alpha/2} + z_\beta)^2 \sigma^2 / \delta^2\)을 \(\delta\)에 대해 풀면 \(\delta \propto 1/\sqrt{n}\)이다. \(n\)을 두 배로 하면 탐지 가능한 최소 차이가 2배가 아니라 \(\sqrt{2} \approx 1.414\)배 작아진다.

직관적으로, 평균 차이의 표준오차는 \(\sigma\sqrt{2/n}\)으로 \(1/\sqrt{n}\)에 따라 줄어든다. \(\delta\)를 탐지하는 검정력은 비 \(\delta / \text{SE}\)에 달려 있으므로, (\(n\)을 두 배로 하여) 정밀도가 \(\sqrt{2}\)배 좋아지면 \(\sqrt{2}\)배 작은 차이를 탐지할 수 있다. 이 수확 체감 관계 때문에 탐지 가능한 효과크기를 절반으로 줄이려면 표본크기가 네 배 필요하다.

필요한 \(n_1\). 분산이 \(\sigma^2(1/n_1+1/n_2)=\sigma^2(1+1/k)/n_1\)이므로

import numpy as np
from scipy import stats

za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
sigma, delta = 10.0, 5.0
base = (za + zb)**2 * sigma**2 / delta**2

print(f"{'k':>5s} {'n1':>6s} {'n2':>6s} {'총':>6s} {'균등 대비':>10s}")
for k in [1, 1.5, 2, 3, 4, 10]:
    n1 = np.ceil((1 + 1 / k) * base)
    n2 = np.ceil(k * (1 + 1 / k) * base)
    print(f"{k:5.1f} {n1:6.0f} {n2:6.0f} {n1 + n2:6.0f} "
          f"{(n1 + n2) / (2 * np.ceil(2 * base)):10.2f}")

    k     n1     n2      총      균등 대비
  1.0     63     63    126       1.00
  1.5     53     79    132       1.05
  2.0     48     95    143       1.13
  3.0     42    126    168       1.33
  4.0     40    157    197       1.56
 10.0     35    346    381       3.02

균등 배분(\(k=1\))이 총 표본을 최소화한다. 총수 \(n_1(1+k)=(1+1/k)(1+k)\cdot\text{base}=(2+k+1/k)\cdot\text{base}\)이고, \(k+1/k\)가 \(k=1\)에서 최소이기 때문이다.

불균등의 대가는 처음에 완만하다. \(k=2\)면 총 표본이 13%만 늘고, \(k=3\)이어야 33%다. \(k\le2\) 정도의 불균형은 실무적으로 큰 손해가 아니다.

그래도 불균등을 쓰는 이유.

1. 한쪽 집단이 부족할 때. 희귀질환 환자를 다 넣고 대조군을 늘리는 경우. \(k=3\)까지는 검정력이 꽤 회복된다. \(k=10\)은 대체로 낭비다.
2. 치료군을 늘려 안전성 자료를 더 모을 때. 신약 시험에서 2:1 배정이 흔하다.
3. 비용이 다를 때.

비용을 넣은 최적 배분. 집단 1의 단위비용 \(c_1\), 집단 2의 \(c_2\)일 때 총비용을 최소화하는 배분은

비용비의 제곱근만큼만 기울인다. 4배 싸다고 4배 뽑지 않는다.

분산이 다르면. \(\sigma_1\ne\sigma_2\)일 때는 \(k^*=\sigma_2/\sigma_1\)(비용이 같을 때)이고, 둘을 합치면 \(k^*=(\sigma_2/\sigma_1)\sqrt{c_1/c_2}\)다.

공식. \(H_0\) 아래의 분산(합동)과 \(H_1\) 아래의 분산을 각각 쓰는 판이 표준이다.

import numpy as np
from scipy import stats

za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
print(f"{'p1':>6s} {'p2':>6s} {'차이':>6s} {'RR':>5s} "
      f"{'정밀식':>8s} {'단순식':>8s}")
for p1, p2 in [(0.30, 0.20), (0.10, 0.05), (0.60, 0.50), (0.05, 0.02)]:
    pb = (p1 + p2) / 2
    v1 = p1 * (1 - p1) + p2 * (1 - p2)
    n = (za * np.sqrt(2 * pb * (1 - pb)) + zb * np.sqrt(v1))**2 / (p1 - p2)**2
    n0 = (za + zb)**2 * v1 / (p1 - p2)**2
    print(f"{p1:6.2f} {p2:6.2f} {p1 - p2:6.2f} {p1 / p2:5.1f} "
          f"{np.ceil(n):8.0f} {np.ceil(n0):8.0f}")

    p1     p2     차이    RR      정밀식      단순식
  0.30   0.20   0.10   1.5      294      291
  0.10   0.05   0.05   2.0      435      432
  0.60   0.50   0.10   1.2      388      385
  0.05   0.02   0.03   2.5      588      586

두 공식의 차이는 1%도 안 된다. 합동분산을 쓰느냐가 큰 문제가 아니다.

무엇이 표본을 결정하는가.

1. 절대 차이가 지배한다. \((0.30,0.20)\)과 \((0.60,0.50)\)은 둘 다 차이 0.10인데 294와 388로 비슷한 규모다. 반면 \((0.10,0.05)\)는 차이가 절반이라 435로 커진다.

2. \(p\)가 0.5에 가까울수록 분산이 커져 불리하다. \((0.60,0.50)\)이 \((0.30,0.20)\)보다 32% 더 필요하다. 같은 절대 차이인데도 그렇다.

3. 상대위험은 직접적인 결정요인이 아니다. \((0.10,0.05)\)는 RR이 2.0으로 크지만 435명이 필요하고, \((0.30,0.20)\)은 RR이 1.5인데 294명이면 된다. 절대 차이가 작으면 아무리 상대위험이 커도 표본이 많이 든다.

희귀사건의 어려움. \(p\)가 작으면 차이도 작아지기 마련이라 표본이 폭발한다. 대략 \(n\approx(z_a+z_b)^2(p_1+p_2)/(p_1-p_2)^2\)이므로, 두 비율을 모두 \(f\)배로 줄이면 \(n\)이 \(1/f\)배로 커진다.

연속성 수정. 카이제곱 검정을 쓸 계획이면 수정판 공식이 있으나, 앞서 보았듯 수정 자체가 지나치게 보수적이므로 수정 없이 계산하고 우도비나 정확검정을 쓰는 편이 낫다.

슈엔펠트 공식. 비례위험 가정 아래 로그순위검정에 필요한 총 사건 수는

(1:1 배정일 때). \(n\)이 아니라 \(d\)가 나온다는 점이 핵심이다.

import numpy as np
from scipy import stats

za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
print(f"{'HR':>6s} {'필요 사건 수':>12s}")
for hr in [0.50, 0.60, 0.70, 0.75, 0.80]:
    print(f"{hr:6.2f} {np.ceil(4 * (za + zb)**2 / np.log(hr)**2):12.0f}")

    HR      필요 사건 수
  0.50           66
  0.60          121
  0.70          247
  0.75          380
  0.80          631

HR이 1에 가까울수록 급격히 커진다. 0.5면 66건, 0.8이면 631건으로 10배다. \((\log\text{HR})^2\)이 분모이기 때문이다.

\(n\)은 어떻게 나오는가. 추적기간 동안 사건이 일어날 확률 \(\psi\)를 추정해

로 환산한다. \(\psi\)는 기저 생존함수, 추적기간, 중도절단 패턴에 달려 있다.

실무적 함의 — 셋.

1. 더 많이 모집하기보다 더 오래 추적하는 것이 나을 수 있다. 사건 수만 채우면 되므로, 추적기간을 늘리면 \(\psi\)가 커져 필요한 \(n\)이 준다. 비용 구조에 따라 어느 쪽이 싼지 달라진다.

2. 예후가 좋은 집단은 연구하기 어렵다. 5년 생존율이 95%인 질환에서는 \(\psi\)가 작아 \(n\)이 폭발한다. 고위험군을 표적으로 삼는 설계가 나오는 이유다.

3. 중간분석의 시점을 사건 수로 정한다. "2년 시점"이 아니라 "사건 100건 시점"으로 계획한다. 사건이 예상보다 적게 일어나면 일정이 밀리는 것을 받아들여야 한다. 이를 "사건 주도(event-driven) 설계"라 한다.

가정의 취약성.

• 비례위험이 깨지면 로그순위검정의 검정력이 떨어진다. 치료효과가 시간에 따라 변하는 면역치료 등에서 문제가 된다.
• 중도절단이 무작위가 아니면 추정이 편향된다.
• HR을 어떻게 정할 것인가. 여기서도 임상적 근거가 필요하며, 낙관적인 HR을 넣어 표본을 줄이는 것이 흔한 실패 원인이다.

설계효과를 곱한다. 앞서 본 대로 \(\text{DEFF}=1+(m-1)\rho\)이므로

import numpy as np
from scipy import stats

za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
sigma, delta, m = 10.0, 5.0, 20
n_ind = np.ceil(2 * (za + zb)**2 * sigma**2 / delta**2)
print(f"개인무작위라면 집단당 {n_ind:.0f}명")
print(f"{'ICC':>6s} {'DEFF':>7s} {'집단당 인원':>11s} {'군집 수':>8s}")
for rho in [0.01, 0.02, 0.05, 0.10]:
    deff = 1 + (m - 1) * rho
    n = n_ind * deff
    print(f"{rho:6.2f} {deff:7.2f} {np.ceil(n):11.0f} {np.ceil(n / m):8.0f}")

개인무작위라면 집단당 63명
   ICC    DEFF      집단당 인원     군집 수
  0.01    1.19          75        4
  0.02    1.38          87        5
  0.05    1.95         123        7
  0.10    2.90         183       10

군집 수 대 군집 크기. \(\rho=0.05\)로 고정하고 \(m\)을 바꿔 보면

rho = 0.05
print(f"{'m':>5s} {'DEFF':>7s} {'총 인원':>8s} {'군집 수':>8s}")
for m in [5, 10, 20, 50, 100]:
    deff = 1 + (m - 1) * rho
    n = n_ind * deff
    print(f"{m:5d} {deff:7.2f} {np.ceil(n):8.0f} {np.ceil(n / m):8.0f}")

    m    DEFF     총 인원     군집 수
    5    1.20       76       16
   10    1.45       92       10
   20    1.95      123        7
   50    3.45      218        5
  100    5.95      375        4

군집을 작게, 많이 뽑는 것이 인원 면에서 효율적이다. \(m=5\)면 76명, \(m=100\)이면 375명으로 5배 차이다.

그러나 군집 수는 반대로 움직인다. \(m=5\)면 16개 군집, \(m=100\)이면 4개다. 학교나 병원을 섭외하는 비용이 크면 군집을 늘리기 어렵다.

최적 \(m\). 군집당 고정비용 \(c_c\), 개인당 비용 \(c_i\)일 때

군집 수의 하한. 통계적으로도 실무적으로도 군집이 너무 적으면 안 된다.

• 군집 수준 자유도가 대략 (군집 수 \(-2\))이므로, 4개 군집이면 자유도 2로 \(t\) 임계값이 4.3에 이른다. 위 계산은 정규근사를 썼으므로 군집 수가 적으면 크게 과소추정한다.
• 군집 수가 적으면 무작위배정만으로 균형이 잡히지 않는다. 층화나 짝지은 배정이 권장된다.
• 실무 지침으로 집단당 최소 10~15개 군집을 권한다.

\(\rho\)를 어디서 얻는가. 이것이 가장 큰 불확실성이다. 교육 연구의 학교 수준 \(\rho\)는 보통 0.01~0.20, 보건의료의 진료소 수준은 0.001~0.05로 보고된다. 민감도 분석이 필수다.

세 가지 설계의 가설.

import numpy as np
from scipy import stats

sigma = 10.0
za = stats.norm.ppf(0.975)          # 단측 0.025 또는 양측 0.05
zb80 = stats.norm.ppf(0.80)
zb90 = stats.norm.ppf(0.90)         # TOST는 각 검정이 더 강해야 함

print(f"{'마진 Δ':>8s} {'비열등성':>10s} {'동등성':>10s}")
for D in [2.0, 3.0, 5.0]:
    n_ni = 2 * (za + zb80)**2 * sigma**2 / D**2
    n_eq = 2 * (za + zb90)**2 * sigma**2 / D**2
    print(f"{D:8.1f} {np.ceil(n_ni):10.0f} {np.ceil(n_eq):10.0f}")

print()
print("참고: 같은 크기의 효과를 탐지하는 우월성 시험")
for d in [2.0, 3.0, 5.0]:
    print(f"  δ = {d}: 집단당 "
          f"{np.ceil(2 * (za + zb80)**2 * sigma**2 / d**2):.0f}")

  마진 Δ      비열등성       동등성
     2.0        393        526
     3.0        175        234
     5.0         63         85

참고: 같은 크기의 효과를 탐지하는 우월성 시험
  δ = 2.0: 집단당 393
  δ = 3.0: 집단당 175
  δ = 5.0: 집단당 63

비열등성과 우월성의 표본이 같다. 마진 \(\Delta\)와 효과 \(\delta\)가 같은 수라면 그렇다. 단측이라 \(\alpha\)가 절반이지만, 관례상 단측 0.025를 쓰므로 양측 0.05와 같은 임계값이다.

동등성이 더 많이 든다. 두 단측검정을 모두 통과해야 하므로 각각의 검정력을 높여야 한다(\(1-\beta/2\) 근사). 34% 더 필요하다.

마진이 결정적이다. \(\Delta\)를 5에서 2로 줄이면 표본이 63명에서 393명으로 6.2배가 된다. \(\Delta\)를 정하는 것이 이 설계의 핵심 의사결정이다.

\(\Delta\)를 어떻게 정하는가.

1. 기존 치료의 효과에서 역산. 위약 대비 기존 치료의 효과가 \(\theta\)라면, 그 효과의 일정 비율(흔히 50%)을 보존하는 것을 마진으로 삼는다.
2. 임상적 판단. "이 정도 차이는 환자에게 무의미하다"는 수준.
3. 규제 지침. 분야별로 정해진 관례가 있다.

절대 하지 말 것 — 자료를 본 뒤 마진 정하기. 어떤 결과든 "동등하다"고 만들 수 있다.

비열등성 시험의 위험 — 편향의 방향이 반대다.

• 우월성 시험에서 자료가 엉성하면(측정오차, 순응도 저하, 탈락) 효과가 희석되어 기각하기 어려워진다. 즉 보수적이다.
• 비열등성 시험에서는 반대다. 엉성한 자료는 두 군을 비슷해 보이게 만들어 비열등성을 주장하기 쉬워진다.

따라서 비열등성 시험은 더 엄격한 품질관리가 필요하고, ITT와 PP 분석을 모두 보고해야 한다.

import numpy as np
from scipy import stats

sigma, delta = 10.0, 5.0
zb = stats.norm.ppf(0.80)
n1 = None
print(f"{'비교 수':>7s} {'보정 α':>10s} {'집단당 n':>9s} {'배율':>6s}")
for m in [1, 2, 3, 5, 10, 20]:
    alpha = 0.05 / m
    za = stats.norm.ppf(1 - alpha / 2)
    n = np.ceil(2 * (za + zb)**2 * sigma**2 / delta**2)
    if n1 is None:
        n1 = n
    print(f"{m:7d} {alpha:10.4f} {n:9.0f} {n / n1:6.2f}")

   비교 수       보정 α     집단당 n     배율
      1     0.0500        63   1.00
      2     0.0250        77   1.22
      3     0.0167        84   1.33
      5     0.0100        94   1.49
     10     0.0050       107   1.70
     20     0.0025       120   1.90

생각보다 완만하다. 비교가 20개로 늘어도 표본은 1.9배다.

왜 그런가. \(n\propto(z_{1-\alpha/2m}+z_{1-\beta})^2\)인데, \(z\)는 \(\alpha\)의 로그에 대략 \(\sqrt{2\log(1/\alpha)}\)로 자란다. \(m\)이 지수적으로 늘어야 \(n\)이 선형으로 는다.

그러나 안심할 일이 아니다.

1. 검정력의 정의를 확인해야 한다. 위 계산은 각 비교의 개별 검정력 80%다. "모든 비교에서 동시에 유의" 같은 다른 정의라면 훨씬 많이 든다.

2. 주 결과를 하나로 정하는 것이 낫다. 여러 결과를 모두 주 결과로 삼으면 해석이 혼란스럽다. 주 결과 하나 + 나머지는 부수적으로 두면 보정이 필요 없다.

대안들.

계층적 검정이 실무에서 유용하다. "먼저 주 결과를 0.05로 검정하고, 유의하면 다음 결과를 0.05로 검정한다"는 규칙은 FWER을 보정 없이 유지한다. 다만 순서를 미리 정해야 하고, 앞에서 막히면 뒤를 볼 수 없다.

복합 결과변수의 함정. "사망 또는 재입원 또는 재발"처럼 묶으면 사건이 늘어 표본이 준다. 하지만 구성요소의 중요도가 다르면 해석이 왜곡된다. 재입원이 늘고 사망이 준 경우 복합 결과는 변화가 없어 보인다.

권고. 표본크기 단계에서 다중비교를 고려하는 것보다, 설계 단계에서 주 결과를 하나로 좁히는 것이 훨씬 효과적이다. 보정은 문제의 증상이지 원인이 아니다.