Bonferroni와 Scheffé 방법¶

import numpy as np
from scipy import stats

print("본페로니 임계값과 셰페 임계값 (α=0.05)")
Ms = [3, 5, 10, 20, 50, 100]
print(f"{'k':>3s} {'ν':>5s} {'셰페 배수':>9s} "
      + " ".join(f"{'M=' + str(m):>8s}" for m in Ms))
for k, nu in [(3, 27), (4, 36), (5, 45), (6, 54), (10, 90)]:
    sc = np.sqrt((k - 1) * stats.f.ppf(0.95, k - 1, nu))
    row = [stats.t.ppf(1 - 0.05 / (2 * m), nu) for m in Ms]
    print(f"{k:3d} {nu:5d} {sc:9.4f} " + " ".join(f"{r:8.4f}" for r in row))

print("\n각 k 에서 본페로니가 셰페를 넘어서는 최소 M")
for k, nu in [(3, 27), (4, 36), (5, 45), (6, 54), (10, 90), (20, 190)]:
    sc = np.sqrt((k - 1) * stats.f.ppf(0.95, k - 1, nu))
    M = 1
    while stats.t.ppf(1 - 0.05 / (2 * M), nu) < sc:
        M += 1
    print(f"  k={k:3d}: 셰페 배수 {sc:.4f}  →  M ≥ {M} 이면 셰페가 유리 "
          f"(쌍별 비교 수는 {k * (k - 1) // 2})")

본페로니 임계값과 셰페 임계값 (α=0.05)
  k     ν     셰페 배수      M=3      M=5     M=10     M=20     M=50    M=100
  3    27    2.5900   2.5525   2.7707   3.0565   3.3334   3.6896   3.9538
  4    36    2.9324   2.5110   2.7195   2.9905   3.2507   3.5821   3.8255
  5    45    3.2117   2.4868   2.6896   2.9521   3.2028   3.5203   3.7519
  6    54    3.4540   2.4708   2.6700   2.9270   3.1716   3.4800   3.7042
 10    90    4.2273   2.4395   2.6316   2.8779   3.1108   3.4019   3.6118

각 k 에서 본페로니가 셰페를 넘어서는 최소 M
  k=  3: 셰페 배수 2.5900  →  M ≥ 4 이면 셰페가 유리 (쌍별 비교 수는 3)
  k=  4: 셰페 배수 2.9324  →  M ≥ 9 이면 셰페가 유리 (쌍별 비교 수는 6)
  k=  5: 셰페 배수 3.2117  →  M ≥ 21 이면 셰페가 유리 (쌍별 비교 수는 10)
  k=  6: 셰페 배수 3.4540  →  M ≥ 47 이면 셰페가 유리 (쌍별 비교 수는 15)
  k= 10: 셰페 배수 4.2273  →  M ≥ 883 이면 셰페가 유리 (쌍별 비교 수는 45)
  k= 20: 셰페 배수 5.5848  →  M ≥ 624455 이면 셰페가 유리 (쌍별 비교 수는 190)

print("시닥 대 본페로니 (α = 0.05)")
print(f"{'M':>5s} {'본페로니 α/M':>12s} {'시닥 1-(1-α)^(1/M)':>18s} {'차이':>8s}")
for M in [2, 3, 5, 10, 20, 50, 100]:
    b = 0.05 / M
    s = 1 - (1 - 0.05)**(1 / M)
    print(f"{M:5d} {b:12.6f} {s:18.6f} {100 * (s - b) / b:7.2f}%")

시닥 대 본페로니 (α = 0.05)
    M     본페로니 α/M   시닥 1-(1-α)^(1/M)       차이
    2     0.025000           0.025321    1.28%
    3     0.016667           0.016952    1.71%
    5     0.010000           0.010206    2.06%
   10     0.005000           0.005116    2.32%
   20     0.002500           0.002561    2.46%
   50     0.001000           0.001025    2.53%
  100     0.000500           0.000513    2.56%

import numpy as np
from scipy import stats

k, n, MSW = 5, 10, 4.0
nu = k * (n - 1)
means = np.array([12.0, 9.5, 11.2, 10.0, 8.6])
contrasts = {
    "C1: μ1-μ2": np.array([1, -1, 0, 0, 0], float),
    "C2: μ3-(μ4+μ5)/2": np.array([0, 0, 1, -0.5, -0.5]),
    "C3: μ4-μ5": np.array([0, 0, 0, 1, -1], float),
}

print(f"MSW = {MSW}, ν = {nu}, 집단당 n = {n}")
print(f"{'대비':>20s} {'추정값':>8s} {'SE':>7s} {'t':>7s} {'원 p':>8s}")
rows = []
for lab, c in contrasts.items():
    est = c @ means
    se = np.sqrt(MSW * (c**2 / n).sum())
    t = est / se
    p = 2 * stats.t.sf(abs(t), nu)
    rows.append((lab, est, se, t, p))
    print(f"{lab:>20s} {est:8.4f} {se:7.4f} {t:7.4f} {p:8.5f}")

M = len(rows)
ps = np.array([r[4] for r in rows])
a_bon = 0.05 / M
a_sid = 1 - 0.95**(1 / M)
sc = np.sqrt((k - 1) * stats.f.ppf(0.95, k - 1, nu))

print(f"\n임계값·수준")
print(f"  본페로니: α' = {a_bon:.5f}, 임계 t = {stats.t.ppf(1 - a_bon / 2, nu):.4f}")
print(f"  시닥    : α' = {a_sid:.5f}, 임계 t = {stats.t.ppf(1 - a_sid / 2, nu):.4f}")
print(f"  셰페    : 임계 배수 = {sc:.4f}")

order = np.argsort(ps)
holm = np.zeros(M, bool)
for r, idx in enumerate(order):
    if ps[idx] < 0.05 / (M - r):
        holm[idx] = True
    else:
        break

print(f"\n{'대비':>20s} {'본페로니':>9s} {'시닥':>7s} {'홀름':>7s} {'셰페':>7s}")
for i, (lab, est, se, t, p) in enumerate(rows):
    print(f"{lab:>20s} {'기각' if p < a_bon else '  -':>9s} "
          f"{'기각' if p < a_sid else '  -':>7s} "
          f"{'기각' if holm[i] else '  -':>7s} "
          f"{'기각' if abs(t) > sc else '  -':>7s}")

MSW = 4.0, ν = 45, 집단당 n = 10
                  대비      추정값      SE       t      원 p
           C1: μ1-μ2   2.5000  0.8944  2.7951  0.00760
    C2: μ3-(μ4+μ5)/2   1.9000  0.7746  2.4529  0.01811
           C3: μ4-μ5   1.4000  0.8944  1.5652  0.12453

임계값·수준
  본페로니: α' = 0.01667, 임계 t = 2.4868
  시닥    : α' = 0.01695, 임계 t = 2.4799
  셰페    : 임계 배수 = 3.2117

                  대비      본페로니      시닥      홀름      셰페
           C1: μ1-μ2        기각      기각      기각       -
    C2: μ3-(μ4+μ5)/2         -       -      기각       -
           C3: μ4-μ5         -       -       -       -

정렬:  p(1)=0.0076,  p(2)=0.0181,  p(3)=0.1246

홀름 1단계: 0.0076 < 0.05/3 = 0.0167  ✓ 기각
홀름 2단계: 0.0181 < 0.05/2 = 0.0250  ✓ 기각
홀름 3단계: 0.1246 < 0.05/1 = 0.0500  ✗ 중단

import warnings
warnings.filterwarnings("ignore")

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

rng = np.random.default_rng(12005)
B = 4_000
k, n = 5, 12
dfe = k * (n - 1)
tc = qsturng(0.95, k, dfe) / np.sqrt(2)
sc = np.sqrt((k - 1) * stats.f.ppf(0.95, k - 1, dfe))
bc = stats.t.ppf(1 - 0.05 / (2 * 10), dfe)

a = b = c = d = 0
for _ in range(B):
    gs = [rng.normal(0, 1, n) for _ in range(k)]
    m = np.array([g.mean() for g in gs])
    MSE = np.array([g.var(ddof=1) for g in gs]).mean()
    tmax = max(abs(m[i] - m[j]) / np.sqrt(2 * MSE / n)
               for i in range(k) for j in range(i + 1, k))
    a += tmax > tc
    b += tmax > bc
    c += tmax > sc
    d += np.sqrt(n * ((m - m.mean())**2).sum() / MSE) > sc

print("k=5, n=12, 완전 귀무. '가장 큰 차이를 내는 것'을 자료에서 찾아 검정")
print(f"{'쌍별 최대 — 튜키 기준':>28s} {a / B:8.4f}")
print(f"{'쌍별 최대 — 본페로니(10쌍)':>28s} {b / B:8.4f}")
print(f"{'쌍별 최대 — 셰페 기준':>28s} {c / B:8.4f}")
print(f"{'최적 대비 — 셰페 기준':>28s} {d / B:8.4f}")
print(f"\n임계값: 튜키 {tc:.4f}, 본페로니 {bc:.4f}, 셰페 {sc:.4f}")

rng = np.random.default_rng(12005)
e = sum(stats.f_oneway(*[rng.normal(0, 1, n) for _ in range(k)]).pvalue < 0.05
        for _ in range(B))
print(f"F 검정의 기각률(참고): {e / B:.4f}")

k=5, n=12, 완전 귀무. '가장 큰 차이를 내는 것'을 자료에서 찾아 검정
               쌍별 최대 — 튜키 기준   0.0590
           쌍별 최대 — 본페로니(10쌍)   0.0435
               쌍별 최대 — 셰페 기준   0.0195
               최적 대비 — 셰페 기준   0.0542

임계값: 튜키 2.8204, 본페로니 2.9247, 셰페 3.1873
F 검정의 기각률(참고): 0.0542

1. 다섯 집단의 평균을 본다
2. "A, B 가 높고 C, D, E 가 낮아 보인다"
3. 대비 (A+B)/2 - (C+D+E)/3 를 검정한다
4. p = 0.02 → "유의하다"고 보고

→ 3단계의 대비를 2단계에서 자료를 보고 골랐다
→ 본페로니(M=1)나 t 검정은 무효
→ 셰페 기준으로 다시 검정해야 한다

import numpy as np
from scipy import stats

k, n, MSW = 5, 10, 4.0
nu = k * (n - 1)
means = np.array([12.0, 9.5, 11.2, 10.0, 8.6])
contrasts = {
    "C1: μ1-μ2": np.array([1, -1, 0, 0, 0], float),
    "C2: μ3-(μ4+μ5)/2": np.array([0, 0, 1, -0.5, -0.5]),
    "C3: μ4-μ5": np.array([0, 0, 0, 1, -1], float),
    "C4: (μ1+μ3)/2-(μ2+μ4+μ5)/3":
        np.array([0.5, -1 / 3, 0.5, -1 / 3, -1 / 3]),
}
S = np.sqrt((k - 1) * stats.f.ppf(0.95, k - 1, nu))
tb = stats.t.ppf(1 - 0.05 / (2 * 3), nu)      # 사전에 정한 3개 기준

print(f"셰페 배수 S = {S:.4f},  본페로니(M=3) t = {tb:.4f}")
print(f"{'대비':>30s} {'추정값':>8s} {'SE':>7s} "
      f"{'셰페 구간':>22s} {'본페로니 구간':>22s}")
for lab, c in contrasts.items():
    est = c @ means
    se = np.sqrt(MSW * (c**2 / n).sum())
    print(f"{lab:>30s} {est:8.3f} {se:7.4f} "
          f"[{est - S * se:8.3f},{est + S * se:8.3f}] "
          f"[{est - tb * se:8.3f},{est + tb * se:8.3f}]")
print(f"\n구간 폭의 비 (셰페/본페로니) = {S / tb:.4f}")

셰페 배수 S = 3.2117,  본페로니(M=3) t = 2.4868
                            대비      추정값      SE                  셰페 구간                본페로니 구간
                     C1: μ1-μ2    2.500  0.8944 [  -0.373,   5.373] [   0.276,   4.724]
              C2: μ3-(μ4+μ5)/2    1.900  0.7746 [  -0.588,   4.388] [  -0.026,   3.826]
                     C3: μ4-μ5    1.400  0.8944 [  -1.473,   4.273] [  -0.824,   3.624]
    C4: (μ1+μ3)/2-(μ2+μ4+μ5)/3    2.233  0.5774 [   0.379,   4.088] [   0.798,   3.669]

구간 폭의 비 (셰페/본페로니) = 1.2915

보고서 구조

[확증적 분석]  사전에 정한 대비 3개, 홀름 보정
   C1: 처리 A 대 대조   p_adj = 0.008  ✓
   C2: 처리 B 대 대조   p_adj = 0.041  ✓
   C3: 처리 A 대 B      p_adj = 0.310

[탐색적 분석]  자료를 보고 발견한 패턴, 셰페 기준
   (A+B)/2 대 (C+D)/2   셰페 구간 [0.29, 4.21], 0 을 포함하지 않음
   → 가설 생성용. 확증에는 새 자료가 필요하다.

import numpy as np
from scipy import stats

k, n, MSW = 4, 10, 4.0
nu = k * (n - 1)
C = {
    "A: μ1-μ2": np.array([1, -1, 0, 0], float),
    "B: μ3-μ4": np.array([0, 0, 1, -1], float),
    "C: (μ1+μ2)-(μ3+μ4)": np.array([1, 1, -1, -1], float),
}
names = list(C)
print("대비 사이의 내적 (0 이면 직교)")
for i in range(3):
    for j in range(i + 1, 3):
        print(f"  {names[i]:>22s} · {names[j]:<22s} = "
              f"{C[names[i]] @ C[names[j]]:.1f}")

print("\n직교 대비의 제곱합이 SSB 를 분할하는가")
rng = np.random.default_rng(15001)
gs = [rng.normal(m, 2.0, n) for m in [10, 12, 9, 11]]
mm = np.array([g.mean() for g in gs])
SSB = n * ((mm - mm.mean())**2).sum()
tot = 0
for lab, c in C.items():
    ss = (c @ mm)**2 / ((c**2).sum() / n)
    tot += ss
    print(f"  {lab:>22s}: SS = {ss:9.4f}")
print(f"  {'합':>22s}: {tot:9.4f}   SSB = {SSB:9.4f}")

print("\n비직교 대비를 섞으면")
D = {"A: μ1-μ2": np.array([1, -1, 0, 0], float),
     "D: μ1-μ3": np.array([1, 0, -1, 0], float),
     "E: μ1-μ4": np.array([1, 0, 0, -1], float)}
tot = 0
for lab, c in D.items():
    ss = (c @ mm)**2 / ((c**2).sum() / n)
    tot += ss
    print(f"  {lab:>22s}: SS = {ss:9.4f}")
print(f"  {'합':>22s}: {tot:9.4f}   SSB = {SSB:9.4f}  (분할되지 않는다)")

대비 사이의 내적 (0 이면 직교)
                A: μ1-μ2 · B: μ3-μ4               = 0.0
                A: μ1-μ2 · C: (μ1+μ2)-(μ3+μ4)     = 0.0
                B: μ3-μ4 · C: (μ1+μ2)-(μ3+μ4)     = 0.0

직교 대비의 제곱합이 SSB 를 분할하는가
                A: μ1-μ2: SS =    6.9669
                B: μ3-μ4: SS =   29.6221
      C: (μ1+μ2)-(μ3+μ4): SS =    6.1408
                       합:   42.7299   SSB =   42.7299

비직교 대비를 섞으면
                A: μ1-μ2: SS =    6.9669
                D: μ1-μ3: SS =    9.9466
                E: μ1-μ4: SS =    5.2386
                       합:   22.1521   SSB =   42.7299  (분할되지 않는다)

import numpy as np
from scipy import stats

rng = np.random.default_rng(15002)
B = 20_000
k, n, M = 4, 12, 6
nu = k * (n - 1)
tb = stats.t.ppf(1 - 0.05 / (2 * M), nu)
t1 = stats.t.ppf(0.975, nu)
pairs = [(i, j) for i in range(k) for j in range(i + 1, k)]

cov_all = cov_each = 0
for _ in range(B):
    gs = [rng.normal(0, 1, n) for _ in range(k)]
    m = np.array([g.mean() for g in gs])
    MSE = np.array([g.var(ddof=1) for g in gs]).mean()
    se = np.sqrt(2 * MSE / n)
    ok_all = ok_each = True
    for i, j in pairs:
        d = m[i] - m[j]
        if abs(d) > tb * se:
            ok_all = False
        if abs(d) > t1 * se:
            ok_each = False
    cov_all += ok_all
    cov_each += ok_each

print(f"k={k}, n={n}, M={M} 쌍, 참 차이는 모두 0")
print(f"  본페로니 t = {tb:.4f},  개별 t = {t1:.4f}")
print(f"  본페로니 동시구간이 여섯 개를 모두 덮을 확률 = {cov_all / B:.4f}")
print(f"  개별 95% 구간이 여섯 개를 모두 덮을 확률   = {cov_each / B:.4f}")
print(f"  (참고: 여섯 검정이 독립이면 0.95^6 = {0.95**6:.4f})")

k=4, n=12, M=6 쌍, 참 차이는 모두 0
  본페로니 t = 2.7628,  개별 t = 2.0154
  본페로니 동시구간이 여섯 개를 모두 덮을 확률 = 0.9592
  개별 95% 구간이 여섯 개를 모두 덮을 확률   = 0.8070
  (참고: 여섯 검정이 독립이면 0.95^6 = 0.7351)

대비를 언제 정했는가?
    │
    ├─ 자료 수집 전 (계획된 대비)
    │     ├─ 검정만 필요 ──→ 홀름 (본페로니를 지배)
    │     └─ 구간도 필요 ──→ 본페로니
    │
    ├─ 자료를 본 뒤 (사후 대비)
    │     └─→ 셰페 (다른 절차는 무효)
    │
    └─ 모든 쌍별 비교
          ├─ 등분산 ──→ 튜키 (가장 정확)
          └─ 이분산 ──→ 게임스-하웰

[확증적] 사전에 정한 대비 3개, 홀름 보정 (FWER = 0.05)
  C1  추정 +2.50  95% 본페로니 구간 [ 0.28,  4.72]  p_adj = 0.023
  C2  추정 +1.90  95% 본페로니 구간 [-0.03,  3.83]  p_adj = 0.036
  C3  추정 +1.40  95% 본페로니 구간 [-0.82,  3.62]  p_adj = 0.125

[탐색적] 자료를 보고 발견한 대비, 셰페 기준
  C4  추정 +2.23  95% 셰페 구간 [ 0.38,  4.09]
  → 가설 생성용이며 확증에는 새 자료가 필요하다.