사후비교 보기¶

import numpy as np
import pandas as pd
from statsmodels.stats.multicomp import pairwise_tukeyhsd

# 세 집단의 참 평균을 10, 12, 15 로 두고 표준편차는 모두 3.5 로 맞췄다.
# A와 B의 차이는 표준편차보다 작고 A와 C의 차이는 그보다 크다.
# Tukey 가 어느 쌍을 갈라내고 어느 쌍을 갈라내지 못하는지 보게 된다.
rng = np.random.default_rng(42)
n = 15
df = pd.DataFrame({
    "response": np.concatenate([
        rng.normal(10.0, 3.5, n),
        rng.normal(12.0, 3.5, n),
        rng.normal(15.0, 3.5, n),
    ]),
    "group": ["A"] * n + ["B"] * n + ["C"] * n,
})

# Tukey HSD 는 쌍 세 개를 한꺼번에 견주면서 전체 오류율을 0.05 로 묶는다.
# 쌍마다 t-검정을 따로 하면 이 통제가 무너진다.
print(pairwise_tukeyhsd(endog=df["response"], groups=df["group"], alpha=0.05))

Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj   lower  upper  reject
---------------------------------------------------
     A      B   2.1353 0.1115 -0.3876 4.6582  False
     A      C   5.4745    0.0  2.9516 7.9975   True
     B      C   3.3392 0.0069  0.8163 5.8621   True
---------------------------------------------------

from scipy import stats

k, sigma = 3, 3.5
nu = k * n - k                      # 42
delta = 2.0                         # A 와 B 의 참 평균차
se_true = sigma * np.sqrt(2 / n)    # 차이의 '참' 표준오차
ncp = delta / se_true
print(f"참 표준오차 = {se_true:.4f},  비중심모수 delta/SE = {ncp:.4f}")

thr = {
    "보정 없음": stats.t.ppf(0.975, nu),
    "Tukey": stats.studentized_range.ppf(0.95, k, nu) / np.sqrt(2),
    "Bonferroni": stats.t.ppf(1 - 0.05 / 6, nu),
    "Scheffe": np.sqrt((k - 1) * stats.f(k - 1, nu).ppf(0.95)),
}

# 모의실험으로 확인한다 (통계량을 직접 만들어 문턱과 견준다).
B = 200_000
sim = np.random.default_rng(7)
mu = np.array([10.0, 12.0, 15.0])
Y = sim.normal(mu[None, :, None], sigma, size=(B, k, n))
gm = Y.mean(axis=2)
MSW = ((Y - gm[:, :, None]) ** 2).sum(axis=(1, 2)) / nu
T_AB = (gm[:, 1] - gm[:, 0]) / np.sqrt(MSW * 2 / n)

print(f"\n{'방법':<12}{'문턱(SE 단위)':>14}{'검정력(이론)':>14}{'검정력(모의)':>14}")
for name, c in thr.items():
    th = stats.nct(nu, ncp).sf(c) + stats.nct(nu, ncp).cdf(-c)
    print(f"{name:<12}{c:>14.4f}{th:>14.4f}{np.mean(np.abs(T_AB) > c):>14.4f}")

# Tukey 로 A-B 를 80% 확률로 잡으려면 집단당 몇 개가 필요한가.
print(f"\n{'n':>5}{'nu':>6}{'Tukey 문턱':>12}{'검정력':>10}")
for m in (15, 30, 50, 60, 63, 70):
    nu_m = k * m - k
    c = stats.studentized_range.ppf(0.95, k, nu_m) / np.sqrt(2)
    nc = delta / (sigma * np.sqrt(2 / m))
    pw = stats.nct(nu_m, nc).sf(c) + stats.nct(nu_m, nc).cdf(-c)
    print(f"{m:>5}{nu_m:>6}{c:>12.4f}{pw:>10.4f}")

print(f"\n이 표본의 MSW = 8.088 로 참 분산 sigma^2 = {sigma ** 2:.3f} 보다 작다")

참 표준오차 = 1.2780,  비중심모수 delta/SE = 1.5649

방법               문턱(SE 단위)       검정력(이론)       검정력(모의)
보정 없음               2.0181        0.3336        0.3340
Tukey               2.4295        0.2056        0.2059
Bonferroni          2.4937        0.1889        0.1895
Scheffe             2.5377        0.1780        0.1780

    n    nu    Tukey 문턱       검정력
   15    42      2.4295    0.2056
   30    87      2.3845    0.4357
   50   147      2.3677    0.6875
   60   177      2.3636    0.7774
   63   186      2.3626    0.7999
   70   207      2.3607    0.8452

이 표본의 MSW = 8.088 로 참 분산 sigma^2 = 12.250 보다 작다

import warnings
warnings.filterwarnings("ignore")

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

def fwer_procs(gs, alpha=0.05):
    k = len(gs)
    n = np.array([len(g) for g in gs], float)
    m = np.array([g.mean() for g in gs])
    v = np.array([g.var(ddof=1) for g in gs])
    MSE = ((n - 1) * v).sum() / (n.sum() - k)
    dfe = int(n.sum() - k)
    ts = np.array([abs(m[i] - m[j]) / np.sqrt(MSE * (1 / n[i] + 1 / n[j]))
                   for i in range(k) for j in range(i + 1, k)])
    ps = 2 * stats.t.sf(ts, dfe)
    M = len(ts)

    srt = np.sort(ps)
    holm = any(srt[idx] < alpha / (M - idx)
               for idx in range(next((i for i in range(M)
                                      if srt[i] >= alpha / (M - i)), M)))
    return {
        "보정 없음": (ps < alpha).any(),
        "본페로니": (ps < alpha / M).any(),
        "홀름": holm,
        "튜키": (ts * np.sqrt(2) > qsturng(1 - alpha, k, dfe)).any(),
        "셰페": (ts**2 > (k - 1) * stats.f.ppf(1 - alpha, k - 1, dfe)).any(),
    }

rng = np.random.default_rng(12001)
B = 4_000
keys = ["보정 없음", "본페로니", "홀름", "튜키", "셰페"]
print("완전 귀무, 정규·등분산, 명목 0.05")
print(f"{'설계':>14s} " + " ".join(f"{s:>9s}" for s in keys))
for k, n in [(3, 15), (4, 15), (6, 12), (10, 10)]:
    acc = {s: 0 for s in keys}
    for _ in range(B):
        r = fwer_procs([rng.normal(0, 1, n) for _ in range(k)])
        for s in keys:
            acc[s] += r[s]
    print(f"k={k}, n={n:<6d} " + " ".join(f"{acc[s] / B:9.4f}" for s in keys))

완전 귀무, 정규·등분산, 명목 0.05
            설계     보정 없음      본페로니        홀름        튜키        셰페
k=3, n=15        0.1230    0.0387    0.0387    0.0470    0.0348
k=4, n=15        0.2050    0.0450    0.0450    0.0537    0.0307
k=6, n=12        0.3560    0.0333    0.0333    0.0475    0.0107
k=10, n=10        0.5920    0.0328    0.0328    0.0460    0.0018

import warnings
warnings.filterwarnings("ignore")

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

def proc_rejects(gs, alpha=0.05):
    """각 절차가 어느 쌍을 기각하는지."""
    k = len(gs)
    n = np.array([len(g) for g in gs], float)
    m = np.array([g.mean() for g in gs])
    v = np.array([g.var(ddof=1) for g in gs])
    MSE = ((n - 1) * v).sum() / (n.sum() - k)
    dfe = int(n.sum() - k)
    pairs = [(i, j) for i in range(k) for j in range(i + 1, k)]
    ts = np.array([abs(m[i] - m[j]) / np.sqrt(MSE * (1 / n[i] + 1 / n[j]))
                   for i, j in pairs])
    ps = 2 * stats.t.sf(ts, dfe)
    M = len(pairs)
    rej = np.zeros(M, bool)
    for r, idx in enumerate(np.argsort(ps)):
        if ps[idx] < alpha / (M - r):
            rej[idx] = True
        else:
            break
    return pairs, {
        "본페로니": ps < alpha / M,
        "홀름": rej,
        "튜키": ts * np.sqrt(2) > qsturng(1 - alpha, k, dfe),
        "셰페": ts**2 > (k - 1) * stats.f.ppf(1 - alpha, k - 1, dfe),
    }

procs = ["본페로니", "홀름", "튜키", "셰페"]
B = 4_000

rng = np.random.default_rng(12002)
print("검정력: k=4, n=15, μ=(0,0,0,δ) — 참인 차이는 3쌍")
print(f"{'δ':>5s} " + " ".join(f"{s + ' 전체':>11s}" for s in procs))
for d in [0.8, 1.2, 1.6]:
    acc = {s: 0 for s in procs}
    for _ in range(B):
        pairs, out = proc_rejects([rng.normal(mu, 1, 15)
                                   for mu in [0, 0, 0, d]])
        idx = [i for i, p in enumerate(pairs) if p[1] == 3]
        for s in procs:
            acc[s] += out[s][idx].all()
    print(f"{d:5.1f} " + " ".join(f"{acc[s] / B:11.4f}" for s in procs))

rng = np.random.default_rng(12003)
print("\n부분 귀무에서의 FWER (참 차이가 없는 3쌍 중 하나라도 기각)")
print(f"{'δ':>5s} " + " ".join(f"{s:>11s}" for s in procs))
for d in [0, 2, 4, 8]:
    acc = {s: 0 for s in procs}
    for _ in range(B):
        pairs, out = proc_rejects([rng.normal(mu, 1, 15)
                                   for mu in [0, 0, 0, d]])
        idx = [i for i, p in enumerate(pairs) if p[1] != 3]
        for s in procs:
            acc[s] += out[s][idx].any()
    print(f"{d:5.1f} " + " ".join(f"{acc[s] / B:11.4f}" for s in procs))

검정력: k=4, n=15, μ=(0,0,0,δ) — 참인 차이는 3쌍
    δ     본페로니 전체       홀름 전체       튜키 전체       셰페 전체
  0.8      0.1067      0.1373      0.1250      0.0815
  1.2      0.4763      0.5417      0.5150      0.4200
  1.6      0.8682      0.8982      0.8838      0.8297

부분 귀무에서의 FWER (참 차이가 없는 3쌍 중 하나라도 기각)
    δ        본페로니          홀름          튜키          셰페
  0.0      0.0255      0.0262      0.0323      0.0177
  2.0      0.0198      0.0398      0.0253      0.0112
  4.0      0.0210      0.0413      0.0255      0.0155
  8.0      0.0222      0.0475      0.0272      0.0127

p(1) ≤ p(2) ≤ ... ≤ p(M) 으로 정렬
  p(1) < α/M     인가?  아니면 모두 기각하지 않고 종료
  p(2) < α/(M-1) 인가?  아니면 여기서 종료
  p(3) < α/(M-2) 인가?  ...

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
    # 자료가 알려 주는 '최적' 대비: c ∝ (m - m̄)
    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

import warnings
warnings.filterwarnings("ignore")

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

def gh_any(gs, alpha=0.05):
    """게임스-하웰: 쌍마다 웰치형 표준오차와 자유도."""
    k = len(gs)
    n = np.array([len(g) for g in gs], float)
    m = np.array([g.mean() for g in gs])
    v = np.array([g.var(ddof=1) for g in gs])
    for i in range(k):
        for j in range(i + 1, k):
            s = v[i] / n[i] + v[j] / n[j]
            df = s**2 / (v[i]**2 / (n[i]**2 * (n[i] - 1))
                         + v[j]**2 / (n[j]**2 * (n[j] - 1)))
            if abs(m[i] - m[j]) / np.sqrt(s / 2) > qsturng(1 - alpha, k, df):
                return True
    return False

def tk_any(gs, alpha=0.05):
    """Tukey-Kramer: 합동 MSE 와 공통 자유도."""
    k = len(gs)
    n = np.array([len(g) for g in gs], float)
    m = np.array([g.mean() for g in gs])
    v = np.array([g.var(ddof=1) for g in gs])
    MSE = ((n - 1) * v).sum() / (n.sum() - k)
    qc = qsturng(1 - alpha, k, n.sum() - k)
    for i in range(k):
        for j in range(i + 1, k):
            if abs(m[i] - m[j]) / np.sqrt(MSE / 2 * (1 / n[i] + 1 / n[j])) > qc:
                return True
    return False

rng = np.random.default_rng(12006)
B = 3_000
print("완전 귀무에서의 FWER (k=4), 명목 0.05")
print(f"{'설계':>30s} {'튜키':>8s} {'G-H':>8s}")
for lab, ns, sds in [("등분산·균형 n=15, σ=1", [15] * 4, [1, 1, 1, 1]),
                     ("이분산·균형 σ=(1,1,1,4)", [15] * 4, [1, 1, 1, 4]),
                     ("역페어링 n=(6,15,15,30)", [6, 15, 15, 30], [4, 2, 1, 1]),
                     ("정페어링 n=(30,15,15,6)", [30, 15, 15, 6], [4, 2, 1, 1]),
                     ("극단 역페어링 n=(5,10,20,40)", [5, 10, 20, 40],
                      [5, 3, 2, 1])]:
    a = b = 0
    for _ in range(B):
        gs = [rng.normal(0, s, n) for n, s in zip(ns, sds)]
        a += tk_any(gs)
        b += gh_any(gs)
    print(f"{lab:>30s} {a / B:8.4f} {b / B:8.4f}")

완전 귀무에서의 FWER (k=4), 명목 0.05
                            설계       튜키      G-H
              등분산·균형 n=15, σ=1   0.0510   0.0520
            이분산·균형 σ=(1,1,1,4)   0.1037   0.0510
           역페어링 n=(6,15,15,30)   0.2960   0.0480
           정페어링 n=(30,15,15,6)   0.0120   0.0523
        극단 역페어링 n=(5,10,20,40)   0.3620   0.0493

무엇을 비교하려는가?
    │
    ├─ 모든 쌍 (k(k-1)/2 개)
    │     ├─ 등분산·균형 ──→ 튜키 HSD
    │     └─ 이분산      ──→ 게임스-하웰
    │
    ├─ 대조군과의 비교만 (k-1 개)
    │     └─→ 더넷 (튜키보다 검정력 18% 높다)
    │
    ├─ 미리 정한 소수의 대비
    │     └─→ 홀름 (본페로니를 지배)
    │
    └─ 자료를 보고 고른 대비
          └─→ 셰페 (다른 절차는 무효)

전체 검정: Welch F(3, 21.4) = 5.82, p = 0.0045
사후검정: 게임스-하웰 (등분산을 가정하지 않으므로)

  비교          차이     95% 구간        보정 p
  A vs B       2.31   [0.42, 4.20]     0.012
  A vs C      -0.14   [-2.03, 1.75]    0.998
  ...