Games-Howell 검정 (분산이 다를 때)¶

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
from statsmodels.stats.libqsturng import qsturng, psturng
import pingouin as pg

def games_howell(gs, names=None, alpha=0.05):
    """게임스-하웰: 쌍마다 웰치형 표준오차와 자유도를 쓴다."""
    k = len(gs)
    names = names or [f"G{i}" for i in range(k)]
    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])
    rows = []
    for i in range(k):
        for j in range(i + 1, k):
            s = v[i] / n[i] + v[j] / n[j]
            se = np.sqrt(s / 2)
            df = s**2 / (v[i]**2 / (n[i]**2 * (n[i] - 1))
                         + v[j]**2 / (n[j]**2 * (n[j] - 1)))
            diff = m[j] - m[i]
            p = float(np.atleast_1d(psturng(abs(diff) / se, k, df))[0])
            h = qsturng(1 - alpha, k, df) * se
            rows.append([names[i], names[j], m[i], m[j], diff,
                         se * np.sqrt(2), df, p, diff - h, diff + h])
    return pd.DataFrame(rows, columns=["A", "B", "mean_A", "mean_B", "diff",
                                       "se", "df", "pval", "lo", "hi"])

rng = np.random.default_rng(14001)
gs = [rng.normal(mu, sd, n) for n, mu, sd in
      [(12, 10.0, 1.0), (15, 12.0, 2.5), (10, 11.0, 4.0), (14, 13.5, 1.5)]]
names = ["A", "B", "C", "D"]
print("자작 구현")
print(games_howell(gs, names).round(4).to_string(index=False))

df = pd.DataFrame({"y": np.concatenate(gs),
                   "g": np.concatenate([[names[i]] * len(g)
                                        for i, g in enumerate(gs)])})
print("\npingouin")
print(pg.pairwise_gameshowell(dv="y", between="g", data=df)
      .round(4).to_string(index=False))

자작 구현
A B  mean_A  mean_B    diff     se      df   pval      lo     hi
A B 10.1831 11.8290  1.6459 0.6694 21.6293 0.0955 -0.2154 3.5071
A C 10.1831 10.9416  0.7585 0.8585 12.1079 0.7968 -1.7873 3.3043
A D 10.1831 13.6749  3.4918 0.5182 23.7497 0.0010  2.0611 4.9225
B C 11.8290 10.9416 -0.8874 0.9831 17.9694 0.7843 -3.6663 1.8915
B D 11.8290 13.6749  1.8459 0.7056 24.4388 0.0671 -0.0982 3.7900
C D 10.9416 13.6749  2.7333 0.8871 13.5407 0.0372  0.1436 5.3230

pingouin
A B  mean_A  mean_B    diff     se       T      df   pval  hedges
A B 10.1831 11.8290 -1.6459 0.6694 -2.4589 21.6293 0.0956 -0.8627
A C 10.1831 10.9416 -0.7585 0.8585 -0.8835 12.1079 0.8134 -0.3875
A D 10.1831 13.6749 -3.4918 0.5182 -6.7383 23.7497 0.0000 -2.5144
B C 11.8290 10.9416  0.8874 0.9831  0.9027 17.9694 0.8036  0.3643
B D 11.8290 13.6749 -1.8459 0.7056 -2.6160 24.4388 0.0670 -0.9318
C D 10.9416 13.6749 -2.7333 0.8871 -3.0812 13.5407 0.0372 -1.3389

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):
    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)
    return any(abs(m[i] - m[j]) / np.sqrt(MSE / 2 * (1 / n[i] + 1 / n[j])) > qc
               for i in range(k) for j in range(i + 1, k))

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

import warnings
warnings.filterwarnings("ignore")

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

def gh_rej_last(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])
    i, j = 0, k - 1
    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)))
    return abs(m[i] - m[j]) / np.sqrt(s / 2) > qsturng(1 - alpha, k, df)

def tk_rej_last(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)
    qc = qsturng(1 - alpha, k, n.sum() - k)
    i, j = 0, k - 1
    return abs(m[i] - m[j]) / np.sqrt(MSE / 2 * (1 / n[i] + 1 / n[j])) > qc

rng = np.random.default_rng(12006)
B = 3_000
print("등분산(σ=1) 자료에서 G0 대 G3 비교의 검정력")
print(f"{'설계':>30s} {'튜키':>8s} {'G-H':>8s} {'손실':>7s}")
for lab, n, d in [("n=10, μ=(0,0,0,1.5)", 10, 1.5),
                  ("n=20, μ=(0,0,0,1.0)", 20, 1.0),
                  ("n=30, μ=(0,0,0,0.8)", 30, 0.8)]:
    a = b = 0
    for _ in range(B):
        gs = [rng.normal(m, 1, n) for m in [0, 0, 0, d]]
        a += tk_rej_last(gs)
        b += gh_rej_last(gs)
    print(f"{lab:>30s} {a / B:8.4f} {b / B:8.4f} "
          f"{100 * (b / B - a / B) / (a / B):6.1f}%")

등분산(σ=1) 자료에서 G0 대 G3 비교의 검정력
                            설계       튜키      G-H      손실
           n=10, μ=(0,0,0,1.5)   0.7337   0.6743   -8.1%
           n=20, μ=(0,0,0,1.0)   0.7010   0.6763   -3.5%
           n=30, μ=(0,0,0,0.8)   0.7070   0.6880   -2.7%

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats
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 t3_any(gs, alpha=0.05):
    """던넷 T3 (시닥 보정 웰치 t 로 근사)."""
    k = len(gs)
    M = k * (k - 1) // 2
    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])
    a2 = 1 - (1 - alpha)**(1 / M)
    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) > stats.t.ppf(1 - a2 / 2, df):
                return True
    return False

rng = np.random.default_rng(14002)
B = 4_000
print("완전 귀무에서의 FWER, 명목 0.05")
print(f"{'설계':>28s} {'게임스-하웰':>11s} {'던넷 T3':>9s}")
for lab, ns, sds in [("k=3, n=5, σ=(1,1,1)", [5] * 3, [1, 1, 1]),
                     ("k=3, n=5, σ=(1,2,4)", [5] * 3, [1, 2, 4]),
                     ("k=3, n=10, σ=(1,2,4)", [10] * 3, [1, 2, 4]),
                     ("k=4, n=6, σ=(1,1,2,4)", [6] * 4, [1, 1, 2, 4]),
                     ("k=4, n=20, σ=(1,1,2,4)", [20] * 4, [1, 1, 2, 4])]:
    a = b = 0
    for _ in range(B):
        gs = [rng.normal(0, s, n) for n, s in zip(ns, sds)]
        a += gh_any(gs)
        b += t3_any(gs)
    print(f"{lab:>28s} {a / B:11.4f} {b / B:9.4f}")

print("\n검정력 (이분산, k=4, n=15, σ=(1,1,1,3), μ4 만 이동)")
print(f"{'δ':>5s} {'게임스-하웰':>11s} {'던넷 T3':>9s}")
for d in [2.0, 3.0, 4.0]:
    a = b = 0
    for _ in range(B):
        gs = [rng.normal(0 if i < 3 else d, s, 15)
              for i, s in enumerate([1, 1, 1, 3])]
        a += gh_any(gs)
        b += t3_any(gs)
    print(f"{d:5.1f} {a / B:11.4f} {b / B:9.4f}")

완전 귀무에서의 FWER, 명목 0.05
                          설계      게임스-하웰     던넷 T3
         k=3, n=5, σ=(1,1,1)      0.0493    0.0355
         k=3, n=5, σ=(1,2,4)      0.0565    0.0435
        k=3, n=10, σ=(1,2,4)      0.0498    0.0430
       k=4, n=6, σ=(1,1,2,4)      0.0517    0.0385
      k=4, n=20, σ=(1,1,2,4)      0.0455    0.0380

검정력 (이분산, k=4, n=15, σ=(1,1,1,3), μ4 만 이동)
    δ      게임스-하웰     던넷 T3
  2.0      0.4970    0.4555
  3.0      0.8525    0.8195
  4.0      0.9838    0.9768

A 대 D:  차이 +3.49,  95% 동시구간 [2.06, 4.92],  p < 0.001
A 대 B:  차이 +1.65,  95% 동시구간 [-0.22, 3.51], p = 0.096

설계: k=4, n=(6,15,15,30), σ=(4,2,1,1)  — 역페어링

웰치 F 검정:      오류율 0.05 근처  ✓  (정확)
튜키 사후검정:     FWER 0.296       ✗  (여섯 배)

→ "전체 검정은 믿을 만한데 사후검정은 여섯 번에 한 번 헛것을 본다"

F 검정: p = 0.02  (유의)
게임스-하웰: 여섯 쌍 모두 p > 0.3

→ "전체는 다른데 어느 쌍도 다르지 않다"는 모순된 보고

1. 집단별 n 과 s 를 본다
     ↓
2. 분산비와 짝짓기 방향을 확인
     │
     ├─ 등분산 · 균형 ──→ 표준 F  +  튜키 HSD
     │
     └─ 그 외 ──────────→ 웰치 F  +  게임스-하웰
     ↓
3. 두 단계를 같은 가정 위에서 수행했음을 보고서에 명시

1. 디자인별 노출 수 n, 클릭 수, 클릭률 p̂, 표준편차를 표로
     ↓
2. 분산비와 짝짓기 방향 확인
     ↓
3. 웰치 분산분석  (또는 이분 자료이므로 로지스틱 회귀)
     ↓
4. 게임스-하웰 사후검정 + 신뢰구간
     ↓
5. 실질적 유의성: 클릭률 차이가 비즈니스에 의미 있는가
     ↓
6. 보고

네 가지 광고 디자인의 클릭률을 비교했다.

디자인   노출      클릭      클릭률    표준오차
  A     12,400    248      2.00%     0.13%p
  B     11,800    401      3.40%     0.17%p
  C     12,100    290      2.40%     0.14%p
  D     12,300    455      3.70%     0.17%p

클릭률이 디자인마다 다르므로 분산도 다르다(분산비 1.9 배).
따라서 웰치 분산분석과 게임스-하웰 사후검정을 사용했다.

Welch F(3, 24_xxx) = 5.12,  p = 0.008

게임스-하웰 (FWER = 0.05)
  비교      차이(%p)   95% 동시구간      p-adj
  B vs A     +1.40   [ 0.83,  1.97]    <0.001
  D vs A     +1.70   [ 1.13,  2.27]    <0.001
  D vs C     +1.30   [ 0.72,  1.88]    <0.001
  B vs C     +1.00   [ 0.43,  1.57]     0.002
  D vs B     +0.30   [-0.26,  0.86]     0.42
  C vs A     +0.40   [-0.15,  0.95]     0.21

결론: B 와 D 가 A·C 보다 명확히 높고, B 와 D 사이에는
      차이의 증거가 없다(구간 [-0.26, 0.86]%p).

실질적 함의: D 를 채택하면 A 대비 클릭률이 1.7%p(상대 85%)
             오르며, 월 100만 노출 기준 클릭 17,000 건 증가에
             해당한다.

pingouin.pairwise_gameshowell(dv="y", between="g", data=df)
  → 차이, se, T, df, pval, hedges
  → 신뢰구간은 주지 않으므로 직접 계산 (연습문제 2)

p < 0.001 은 psturng 의 하한이므로 "0.001 미만"으로만 읽는다

집단별 표준편차: 1.02, 2.48, 4.13, 1.51 (분산비 16 배)
→ 등분산을 가정하지 않고 웰치 분산분석 + 게임스-하웰 사용

  비교      차이    95% 동시구간     p-adj    자유도
  A vs D   +3.49   [2.06, 4.92]    <0.001   23.7
  ...