Dunnett 검정 (대조군과의 비교)¶

import numpy as np
from scipy import stats
from statsmodels.stats.libqsturng import qsturng

def dunnett_crit(k, dfe, alpha=0.05, rho=0.5, M=200_000, rng=None):
    """등상관 ρ 인 다변량 t 의 최대 |t| 분위수를 모의실험으로 구한다."""
    p = k - 1
    z = rng.standard_normal((M, p))
    z0 = rng.standard_normal((M, 1))
    x = np.sqrt(rho) * z0 + np.sqrt(1 - rho) * z      # 등상관 구조
    chi = rng.chisquare(dfe, (M, 1))
    t = x / np.sqrt(chi / dfe)
    return np.quantile(np.abs(t).max(1), 1 - alpha)

print("임계값 비교 (α = 0.05, 양측)")
print(f"{'k':>3s} {'n':>4s} {'ν':>5s} {'더넷':>8s} {'본페로니':>9s} "
      f"{'튜키 q/√2':>10s} {'보정 없음':>9s}")
for k, n in [(4, 10), (5, 12), (6, 15), (8, 20)]:
    dfe = k * (n - 1)
    M = k - 1
    dc = dunnett_crit(k, dfe, rng=np.random.default_rng(7))
    bc = stats.t.ppf(1 - 0.05 / (2 * M), dfe)
    tc = qsturng(0.95, k, dfe) / np.sqrt(2)
    nc = stats.t.ppf(0.975, dfe)
    print(f"{k:3d} {n:4d} {dfe:5d} {dc:8.4f} {bc:9.4f} {tc:10.4f} {nc:9.4f}")

임계값 비교 (α = 0.05, 양측)
  k    n     ν       더넷      본페로니    튜키 q/√2     보정 없음
  4   10    36   2.4475    2.5110     2.6933    2.0281
  5   12    55   2.5174    2.5825     2.8204    2.0040
  6   15    84   2.5619    2.6356     2.9165    1.9886
  8   20   152   2.6470    2.7270     3.0735    1.9757

import numpy as np
from scipy import stats
from statsmodels.stats.libqsturng import qsturng

def dunnett_crit(k, dfe, alpha=0.05, rho=0.5, M=200_000, rng=None):
    p = k - 1
    z = rng.standard_normal((M, p))
    z0 = rng.standard_normal((M, 1))
    x = np.sqrt(rho) * z0 + np.sqrt(1 - rho) * z
    chi = rng.chisquare(dfe, (M, 1))
    return np.quantile(np.abs(x / np.sqrt(chi / dfe)).max(1), 1 - alpha)

rng = np.random.default_rng(12004)
B = 5_000
for k, n in [(4, 10), (5, 12), (6, 15)]:
    dfe = k * (n - 1)
    M = k - 1
    dc = dunnett_crit(k, dfe, rng=np.random.default_rng(7))
    bc = stats.t.ppf(1 - 0.05 / (2 * M), dfe)
    tc = qsturng(0.95, k, dfe) / np.sqrt(2)
    a = b = c = pa = pb = pc = 0
    for _ in range(B):
        gs = [rng.normal(0, 1, n) for _ in range(k)]     # 완전 귀무
        m = np.array([g.mean() for g in gs])
        se = np.sqrt(2 * np.array([g.var(ddof=1) for g in gs]).mean() / n)
        ts = np.abs(m[1:] - m[0]) / se
        a += (ts > dc).any()
        b += (ts > bc).any()
        c += (ts > tc).any()

        gs = [rng.normal(0 if i < k - 1 else 1.2, 1, n) for i in range(k)]
        m = np.array([g.mean() for g in gs])
        se = np.sqrt(2 * np.array([g.var(ddof=1) for g in gs]).mean() / n)
        ts = np.abs(m[1:] - m[0]) / se
        pa += ts[-1] > dc
        pb += ts[-1] > bc
        pc += ts[-1] > tc
    print(f"k={k}, n={n}:  임계값 더넷 {dc:.3f} / 본페로니 {bc:.3f} "
          f"/ 튜키 {tc:.3f}")
    print(f"          FWER  더넷 {a / B:.4f} / 본페로니 {b / B:.4f} "
          f"/ 튜키 {c / B:.4f}")
    print(f"          검정력 더넷 {pa / B:.4f} / 본페로니 {pb / B:.4f} "
          f"/ 튜키 {pc / B:.4f}")

k=4, n=10:  임계값 더넷 2.447 / 본페로니 2.511 / 튜키 2.693
          FWER  더넷 0.0536 / 본페로니 0.0458 / 튜키 0.0300
          검정력 더넷 0.5938 / 본페로니 0.5688 / 튜키 0.5006
k=5, n=12:  임계값 더넷 2.517 / 본페로니 2.583 / 튜키 2.820
          FWER  더넷 0.0518 / 본페로니 0.0446 / 튜키 0.0238
          검정력 더넷 0.6540 / 본페로니 0.6284 / 튜키 0.5410
k=6, n=15:  임계값 더넷 2.562 / 본페로니 2.636 / 튜키 2.916
          FWER  더넷 0.0470 / 본페로니 0.0380 / 튜키 0.0164
          검정력 더넷 0.7640 / 본페로니 0.7408 / 튜키 0.6482

import numpy as np
from scipy import stats

def dunnett_crit(k, dfe, alpha=0.05, rho=0.5, two_sided=True,
                 M=200_000, rng=None):
    p = k - 1
    z = rng.standard_normal((M, p))
    z0 = rng.standard_normal((M, 1))
    x = np.sqrt(rho) * z0 + np.sqrt(1 - rho) * z
    t = x / np.sqrt(rng.chisquare(dfe, (M, 1)) / dfe)
    s = np.abs(t) if two_sided else t
    return np.quantile(s.max(1), 1 - alpha)

print("단측 대 양측 더넷 (α = 0.05)")
print(f"{'k':>3s} {'ν':>5s} {'양측':>8s} {'단측':>8s} {'절감':>7s} "
      f"{'단측 본페로니':>12s}")
for k, n in [(3, 12), (4, 10), (5, 12), (6, 15)]:
    dfe = k * (n - 1)
    M = k - 1
    d2 = dunnett_crit(k, dfe, two_sided=True, rng=np.random.default_rng(11))
    d1 = dunnett_crit(k, dfe, two_sided=False, rng=np.random.default_rng(11))
    b1 = stats.t.ppf(1 - 0.05 / M, dfe)
    print(f"{k:3d} {dfe:5d} {d2:8.4f} {d1:8.4f} "
          f"{100 * (d2 - d1) / d2:6.2f}% {b1:12.4f}")

print("\n표본 크기로 환산 (같은 검정력에 필요한 n 의 비)")
for k, n in [(4, 10), (6, 15)]:
    dfe = k * (n - 1)
    d2 = dunnett_crit(k, dfe, two_sided=True, rng=np.random.default_rng(11))
    d1 = dunnett_crit(k, dfe, two_sided=False, rng=np.random.default_rng(11))
    print(f"  k={k}: (단측/양측)² = {(d1 / d2)**2:.4f}  "
          f"→ 표본을 {100 * (1 - (d1 / d2)**2):.1f}% 줄일 수 있다")

단측 대 양측 더넷 (α = 0.05)
  k     ν       양측       단측      절감      단측 본페로니
  3    33   2.3070   1.9862  13.90%       2.0345
  4    36   2.4459   2.1328  12.80%       2.2129
  5    55   2.5155   2.2115  12.08%       2.3044
  6    84   2.5615   2.2671  11.50%       2.3716

표본 크기로 환산 (같은 검정력에 필요한 n 의 비)
  k=4: (단측/양측)² = 0.7604  → 표본을 24.0% 줄일 수 있다
  k=6: (단측/양측)² = 0.7833  → 표본을 21.7% 줄일 수 있다

import numpy as np
from scipy import stats

def dunnett_crit(k, dfe, ns, alpha=0.05, M=100_000, rng=None):
    """불균형 설계의 더넷 임계값. 상관이 쌍마다 다르다."""
    p = k - 1
    n0 = ns[0]
    rho = np.array([[np.sqrt(ns[i + 1] * ns[j + 1]
                             / ((n0 + ns[i + 1]) * (n0 + ns[j + 1])))
                     if i != j else 1.0 for j in range(p)] for i in range(p)])
    L = np.linalg.cholesky(rho + 1e-12 * np.eye(p))
    x = rng.standard_normal((M, p)) @ L.T
    t = x / np.sqrt(rng.chisquare(dfe, (M, 1)) / dfe)
    return np.quantile(np.abs(t).max(1), 1 - alpha)

k, N = 5, 100
print(f"k={k}개 집단(대조 1 + 처치 4), 총 N={N}. 최적 비 √(k-1) = {np.sqrt(k-1):.3f}")
print(f"{'배분 (n0, n_i)':>20s} {'비':>6s} {'임계값':>8s} {'SE 배수':>8s} "
      f"{'최소 탐지 차이':>13s}")
best = None
for n0 in [12, 16, 20, 24, 28, 32, 36, 40]:
    ni = (N - n0) // (k - 1)
    if ni < 4:
        continue
    ns = [n0] + [ni] * (k - 1)
    dfe = sum(ns) - k
    dc = dunnett_crit(k, dfe, ns, rng=np.random.default_rng(21))
    se = np.sqrt(1 / n0 + 1 / ni)          # σ=1 기준
    mdd = dc * se
    star = " ←" if best is None or mdd < best[1] else ""
    if best is None or mdd < best[1]:
        best = (n0, mdd)
    print(f"{str((n0, ni)):>20s} {n0 / ni:6.2f} {dc:8.4f} {se:8.4f} "
          f"{mdd:13.4f}{star}")
print(f"\n최적 배분: n0 = {best[0]}, 최소 탐지 차이 = {best[1]:.4f}σ")

k=5개 집단(대조 1 + 처치 4), 총 N=100. 최적 비 √(k-1) = 2.000
        배분 (n0, n_i)      비      임계값    SE 배수      최소 탐지 차이
            (12, 22)   0.55   2.4329   0.3589        0.8731 ←
            (16, 21)   0.76   2.4609   0.3318        0.8166 ←
            (20, 20)   1.00   2.4786   0.3162        0.7838 ←
            (24, 19)   1.26   2.4937   0.3071        0.7658 ←
            (28, 18)   1.56   2.5044   0.3021        0.7566 ←
            (32, 17)   1.88   2.5123   0.3001        0.7540 ←
            (36, 16)   2.25   2.5196   0.3005        0.7570
            (40, 15)   2.67   2.5243   0.3028        0.7643

최적 배분: n0 = 32, 최소 탐지 차이 = 0.7540σ

import numpy as np
from scipy import stats

def dunnett_crit(k, dfe, alpha=0.05, rho=0.5, M=200_000, rng=None):
    p = k - 1
    z = rng.standard_normal((M, p))
    z0 = rng.standard_normal((M, 1))
    x = np.sqrt(rho) * z0 + np.sqrt(1 - rho) * z
    t = x / np.sqrt(rng.chisquare(dfe, (M, 1)) / dfe)
    return np.quantile(np.abs(t).max(1), 1 - alpha)

k, n, MSW = 4, 10, 5.2
dfe = k * (n - 1)
means = np.array([12.0, 14.8, 11.5, 16.1])
names = ["대조", "처치1", "처치2", "처치3"]
se = np.sqrt(2 * MSW / n)
dc = dunnett_crit(k, dfe, rng=np.random.default_rng(31))
bc = stats.t.ppf(1 - 0.05 / (2 * (k - 1)), dfe)

print(f"MSW = {MSW}, ν = {dfe}, SE(차이) = {se:.4f}")
print(f"더넷 임계값 = {dc:.4f},  본페로니 = {bc:.4f}")
print(f"\n{'비교':>16s} {'차이':>8s} {'t':>8s} "
      f"{'더넷 95% 구간':>22s} {'판정':>6s}")
for i in range(1, k):
    d = means[i] - means[0]
    t = d / se
    print(f"{names[i] + ' - 대조':>16s} {d:8.2f} {t:8.4f} "
          f"[{d - dc * se:9.3f},{d + dc * se:9.3f}] "
          f"{'기각' if abs(t) > dc else '  -':>6s}")

MSW = 5.2, ν = 36, SE(차이) = 1.0198
더넷 임계값 = 2.4549,  본페로니 = 2.5110

              비교       차이        t              더넷 95% 구간     판정
        처치1 - 대조     2.80   2.7456 [    0.296,    5.304]     기각
        처치2 - 대조    -0.50  -0.4903 [   -3.004,    2.004]      -
        처치3 - 대조     4.10   4.0204 [    1.596,    6.604]     기각

위약군(대조)과 세 가지 처치군을 비교했다(각 군 n = 10).

일원배치 분산분석:  MSW = 5.2,  ν = 36

대조군과의 비교에만 관심이 있으므로 더넷 검정을 사용했다
(모든 쌍을 비교하는 튜키 HSD 보다 검정력이 약 19% 높다).

  비교              차이    95% 동시구간        판정
  처치1 vs 대조     +2.80   [ 0.30,  5.30]     유의
  처치2 vs 대조     -0.50   [-3.00,  2.00]     유의하지 않음
  처치3 vs 대조     +4.10   [ 1.60,  6.60]     유의

처치 1 과 3 이 대조군보다 반응이 높았다. 처치 2 는 대조군과
구별되지 않았으나, 구간이 [-2.99, 3.00]{=5.0} 로 넓어
3 수준의 차이를 배제하지는 못한다.

위약 vs 약1, 약2, 약3        (3개)
표준치료 vs 약1, 약2, 약3     (3개)
                              ─────
                               6개

잘못된 방법: 각 시점마다 더넷을 따로 수행
   → 시점 간 상관을 무시
   → 다중성도 통제되지 않음

올바른 방법:
   1. 주요 시점 하나를 사전에 지정 (주 평가변수)
   2. 나머지는 부차 분석으로 표시
   또는
   3. 혼합효과 모형 + 시점별 대비

H1: 고용량 vs 위약   →  p = 0.003  ✓ 계속
H2: 중용량 vs 위약   →  p = 0.021  ✓ 계속
H3: 저용량 vs 위약   →  p = 0.180  ✗ 중단

→ H1, H2 를 α = 0.05 그대로 검정해도 FWER 이 0.05 로 통제된다

초기 배분:  H1: 0.03,  H2: 0.02,  H3: 0.00

H1 기각 → H1 의 0.03 을 H2, H3 에 분배
          H2: 0.02 + 0.015 = 0.035
          H3: 0.00 + 0.015 = 0.015

scipy >= 1.11:  scipy.stats.dunnett(*treatments, control=ctrl)
R:              multcomp::glht(fit, linfct = mcp(group = "Dunnett"))
직접 계산:       다변량 t 의 최대 |t| 분위수를 모의실험으로
                (연습문제 4의 dunnett_crit)

위약군과 세 처치군을 비교했다(각 n = 10).
대조군과의 비교에만 관심이 있으므로 더넷 검정(양측, FWER = 0.05)을
사용했다. 모든 쌍을 비교하는 튜키 HSD 보다 검정력이 약 19% 높다.

  처치1 vs 위약   +2.80  [ 0.30,  5.30]   유의
  처치2 vs 위약   -0.50  [-3.00,  2.00]   유의하지 않음
  처치3 vs 위약   +4.10  [ 1.60,  6.60]   유의