Welch의 일원배치 분산분석¶

import pingouin as pg
import pandas as pd

# 예시 자료
data = {
    "Group": ["A", "A", "A", "B", "B", "B", "C", "C", "C", "C"],
    "Values": [12, 14, 13, 22, 23, 19, 31, 33, 29, 35],
}
df = pd.DataFrame(data)

# Welch 분산분석 수행. scipy에는 없고 pingouin에 있다.
anova_results = pg.welch_anova(dv="Values", between="Group", data=df)

print(anova_results)

  Source  ddof1     ddof2          F     p_unc       np2
0  Group      2  4.142737  84.415503  0.000439  0.953739

# Games-Howell 사후검정. Tukey HSD와 달리 쌍마다 자유도를 따로 계산한다.
post_hoc = pg.pairwise_gameshowell(dv="Values", between="Group", data=df)
print(post_hoc.round(4).to_string(index=False))

A B  mean_A  mean_B     diff     se        T     df   pval  hedges
A B 13.0000 21.3333  -8.3333 1.3333  -6.2500 2.8764 0.0188 -4.0825
A C 13.0000 32.0000 -19.0000 1.4142 -13.4350 4.0755 0.0004 -7.6277
B C 21.3333 32.0000 -10.6667 1.7638  -6.0474 4.9154 0.0044 -3.7514

import numpy as np
from itertools import combinations

from scipy import stats
from statsmodels.stats.multicomp import pairwise_tukeyhsd

g = df.groupby('Group')['Values'].agg(['mean', 'var', 'count'])
N, k = len(df), 3
MSW = ((g['count'] - 1) * g['var']).sum() / (N - k)
print(g.round(4))
print(f"합동 MSW = {MSW:.4f}  (자유도 {N - k})")

print(f"\n{'pair':<8}{'diff':>10}{'SE(GH)':>9}{'SE(TK)':>9}{'T(GH)':>9}{'T(TK)':>9}"
      f"{'df(GH)':>9}{'p(GH)':>9}{'p(TK)':>9}")
for a, b in combinations('ABC', 2):
    ma, va, na = g.loc[a]
    mb, vb, nb = g.loc[b]
    d = ma - mb
    se_gh = np.sqrt(va / na + vb / nb)
    se_tk = np.sqrt(MSW * (1 / na + 1 / nb))
    nu = (va / na + vb / nb) ** 2 / ((va / na) ** 2 / (na - 1)
                                     + (vb / nb) ** 2 / (nb - 1))
    p_gh = stats.studentized_range.sf(np.sqrt(2) * abs(d / se_gh), k, nu)
    p_tk = stats.studentized_range.sf(np.sqrt(2) * abs(d / se_tk), k, N - k)
    print(f"{a + '-' + b:<8}{d:>10.4f}{se_gh:>9.4f}{se_tk:>9.4f}{d / se_gh:>9.4f}"
          f"{d / se_tk:>9.4f}{nu:>9.4f}{p_gh:>9.4f}{p_tk:>9.4f}")

print("\nstatsmodels 의 Tukey-Kramer")
print(pairwise_tukeyhsd(df['Values'], df['Group'], alpha=0.05))

# hedges 열 재현
print(f"\n{'pair':<8}{'s_p':>9}{'d/s_p':>9}{'J':>9}{'hedges g':>10}")
for a, b in combinations('ABC', 2):
    ma, va, na = g.loc[a]
    mb, vb, nb = g.loc[b]
    dfp = na + nb - 2
    sp = np.sqrt(((na - 1) * va + (nb - 1) * vb) / dfp)
    J = 1 - 3 / (4 * dfp - 1)
    print(f"{a + '-' + b:<8}{sp:>9.4f}{(ma - mb) / sp:>9.4f}{J:>9.4f}"
          f"{J * (ma - mb) / sp:>10.4f}")

          mean     var  count
Group
A      13.0000  1.0000      3
B      21.3333  4.3333      3
C      32.0000  6.6667      4
합동 MSW = 4.3810  (자유도 7)

pair          diff   SE(GH)   SE(TK)    T(GH)    T(TK)   df(GH)    p(GH)    p(TK)
A-B        -8.3333   1.3333   1.7090  -6.2500  -4.8762   2.8764   0.0188   0.0044
A-C       -19.0000   1.4142   1.5986 -13.4350 -11.8853   4.0755   0.0004   0.0000
B-C       -10.6667   1.7638   1.5986  -6.0474  -6.6725   4.9154   0.0044   0.0007

statsmodels 의 Tukey-Kramer
Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj  lower   upper  reject
---------------------------------------------------
     A      B   8.3333 0.0044 3.3003 13.3664   True
     A      C     19.0    0.0 14.292  23.708   True
     B      C  10.6667 0.0007 5.9587 15.3747   True
---------------------------------------------------

pair          s_p    d/s_p        J  hedges g
A-B        1.6330  -5.1031   0.8000   -4.0825
A-C        2.0976  -9.0579   0.8421   -7.6277
B-C        2.3944  -4.4548   0.8421   -3.7514

import numpy as np
from scipy import stats

n = np.array([36.0, 24.0, 48.0])
m = np.array([1.8, 1.2, 1.0])
v = np.array([12.5, 3.1, 5.8])
k = len(n)

# --- 웰치 ---
w = n / v
W = w.sum()
mt = (w * m).sum() / W
A = (w * (m - mt)**2).sum() / (k - 1)
lam = ((1 - w / W)**2 / (n - 1)).sum()
B = 1 + 2 * (k - 2) / (k**2 - 1) * lam
F_w = A / B
df2 = (3 / (k**2 - 1) * lam)**-1
print("가중치 w:", np.round(w, 4), "  W =", round(W, 4))
print(f"가중 전체평균 = {mt:.4f}")
print(f"분자 A = {A:.6f},  보정항 B = {B:.6f},  λ = {lam:.6f}")
print(f"F_W = {F_w:.4f},  df1 = {k - 1},  df2 = {df2:.4f},  "
      f"p = {stats.f.sf(F_w, k - 1, df2):.4f}")

# --- 표준 F ---
N = n.sum()
gm = (n * m).sum() / N
SSB = (n * (m - gm)**2).sum()
SSW = ((n - 1) * v).sum()
F_c = (SSB / (k - 1)) / (SSW / (N - k))
print(f"\n표준 F = {F_c:.4f},  df = ({k - 1}, {int(N - k)}),  "
      f"p = {stats.f.sf(F_c, k - 1, N - k):.4f}")

가중치 w: [2.88   7.7419 8.2759]   W = 18.8978
가중 전체평균 = 1.2039
분자 A = 0.683777,  보정항 B = 1.010600,  λ = 0.042400
F_W = 0.6766,  df1 = 2,  df2 = 62.8933,  p = 0.5120

표준 F = 0.9102,  df = (2, 105),  p = 0.4056

import numpy as np
import pandas as pd
from scipy import stats
import pingouin as pg

def welch_anova_summary(n, m, v):
    """요약통계(n, 평균, 분산)만으로 웰치 F 와 자유도를 계산한다."""
    n, m, v = map(np.asarray, (n, m, v))
    k = len(n)
    w = n / v
    W = w.sum()
    mt = (w * m).sum() / W
    lam = ((1 - w / W)**2 / (n - 1)).sum()
    F = ((w * (m - mt)**2).sum() / (k - 1)
         / (1 + 2 * (k - 2) / (k**2 - 1) * lam))
    df2 = (3 / (k**2 - 1) * lam)**-1
    return F, k - 1, df2, stats.f.sf(F, k - 1, df2)

# 검증: 평균과 분산을 정확히 지정한 자료를 만들어 대조한다
rng = np.random.default_rng(5)
rows = []
for i, (ni, mi, vi) in enumerate(zip([30, 20, 40], [0, 0.5, 1.0],
                                     [4.0, 1.0, 2.0])):
    x = rng.normal(0, 1, ni)
    x = (x - x.mean()) / x.std(ddof=1) * np.sqrt(vi) + mi   # 표준화 후 재척도
    rows.append(pd.DataFrame({"y": x, "g": f"G{i}"}))
df = pd.concat(rows, ignore_index=True)

g = df.groupby("g").y.agg(["size", "mean", "var"])
F, d1, d2, p = welch_anova_summary(g["size"], g["mean"], g["var"])
print(f"자작 구현:   F = {F:.6f},  df2 = {d2:.6f},  p = {p:.6f}")

r = pg.welch_anova(dv="y", between="g", data=df)
print(f"pingouin:  F = {r.F[0]:.6f},  df2 = {r.ddof2[0]:.6f},  "
      f"p = {r.p_unc[0]:.6f}")

자작 구현:   F = 2.988446,  df2 = 52.608843,  p = 0.058981
pingouin:  F = 2.988446,  df2 = 52.608843,  p = 0.058981

x ← (x - x̄) / s × √v + m

import numpy as np
from scipy import stats

def welch_p(groups):
    n = np.array([len(g) for g in groups], float)
    m = np.array([g.mean() for g in groups])
    v = np.array([g.var(ddof=1) for g in groups])
    k = len(n)
    w = n / v
    W = w.sum()
    mt = (w * m).sum() / W
    lam = ((1 - w / W)**2 / (n - 1)).sum()
    F = ((w * (m - mt)**2).sum() / (k - 1)
         / (1 + 2 * (k - 2) / (k**2 - 1) * lam))
    return stats.f.sf(F, k - 1, (3 / (k**2 - 1) * lam)**-1)

rng = np.random.default_rng(20240)
M = 2_000
print(f"{'설계':>32s} {'표준 F':>8s} {'웰치':>8s}")
cases = [
    ("등분산 n=(20,20,20) σ=(1,1,1)", [20, 20, 20], [1, 1, 1]),
    ("이분산 n=(20,20,20) σ=(1,1,3)", [20, 20, 20], [1, 1, 3]),
    ("정페어링 n=(10,20,40) σ=(1,2,3)", [10, 20, 40], [1, 2, 3]),
    ("역페어링 n=(10,20,40) σ=(3,2,1)", [10, 20, 40], [3, 2, 1]),
    ("역페어링 n=(5,10,40)  σ=(4,2,1)", [5, 10, 40], [4, 2, 1]),
]
for lab, ns, sds in cases:
    a = b = 0
    for _ in range(M):
        gs = [rng.normal(0, s, n) for n, s in zip(ns, sds)]
        a += stats.f_oneway(*gs).pvalue < 0.05
        b += welch_p(gs) < 0.05
    print(f"{lab:>32s} {a / M:8.4f} {b / M:8.4f}")

                              설계     표준 F       웰치
      등분산 n=(20,20,20) σ=(1,1,1)   0.0525   0.0540
      이분산 n=(20,20,20) σ=(1,1,3)   0.0780   0.0565
     정페어링 n=(10,20,40) σ=(1,2,3)   0.0150   0.0480
     역페어링 n=(10,20,40) σ=(3,2,1)   0.1930   0.0475
     역페어링 n=(5,10,40)  σ=(4,2,1)   0.3745   0.0635

import numpy as np
from scipy import stats

def welch_p(groups):
    n = np.array([len(g) for g in groups], float)
    m = np.array([g.mean() for g in groups])
    v = np.array([g.var(ddof=1) for g in groups])
    k = len(n)
    w = n / v
    W = w.sum()
    mt = (w * m).sum() / W
    lam = ((1 - w / W)**2 / (n - 1)).sum()
    F = ((w * (m - mt)**2).sum() / (k - 1)
         / (1 + 2 * (k - 2) / (k**2 - 1) * lam))
    return stats.f.sf(F, k - 1, (3 / (k**2 - 1) * lam)**-1)

rng = np.random.default_rng(555)
M = 2_000
print(f"{'설계 (모두 등분산 σ=1)':>32s} {'표준 F':>8s} {'웰치':>8s} {'손실':>8s}")
for lab, ns, mus in [
    ("n=(10,10,10) μ=(0,0.5,1)", [10, 10, 10], [0, 0.5, 1.0]),
    ("n=(20,20,20) μ=(0,0.5,1)", [20, 20, 20], [0, 0.5, 1.0]),
    ("n=(5,5,5)    μ=(0,1,2)", [5, 5, 5], [0, 1, 2]),
    ("n=(10,10,10,10,10) μ=0..1", [10] * 5, [0, 0.25, 0.5, 0.75, 1.0]),
]:
    a = b = 0
    for _ in range(M):
        gs = [rng.normal(u, 1, n) for n, u in zip(ns, mus)]
        a += stats.f_oneway(*gs).pvalue < 0.05
        b += welch_p(gs) < 0.05
    print(f"{lab:>32s} {a / M:8.4f} {b / M:8.4f} {(a - b) / max(a, 1):8.4f}")

                 설계 (모두 등분산 σ=1)     표준 F       웰치       손실
        n=(10,10,10) μ=(0,0.5,1)   0.4690   0.4455   0.0501
        n=(20,20,20) μ=(0,0.5,1)   0.7830   0.7780   0.0064
          n=(5,5,5)    μ=(0,1,2)   0.7060   0.6210   0.1204
       n=(10,10,10,10,10) μ=0..1   0.4430   0.4090   0.0767

import numpy as np
from scipy import stats

def welch_F(groups):
    n = np.array([len(g) for g in groups], float)
    m = np.array([g.mean() for g in groups])
    v = np.array([g.var(ddof=1) for g in groups])
    k = len(n)
    w = n / v
    W = w.sum()
    mt = (w * m).sum() / W
    lam = ((1 - w / W)**2 / (n - 1)).sum()
    F = ((w * (m - mt)**2).sum() / (k - 1)
         / (1 + 2 * (k - 2) / (k**2 - 1) * lam))
    df2 = (3 / (k**2 - 1) * lam)**-1
    return F, df2, stats.f.sf(F, k - 1, df2)

rng = np.random.default_rng(777)
print(f"{'':4s} {'F_W':>10s} {'t²':>10s} {'df2':>9s} "
      f"{'t 의 df':>9s} {'p(F)':>9s} {'p(t)':>9s}")
for i in range(4):
    a = rng.normal(0, 1, rng.integers(8, 25))
    b = rng.normal(0.6, 2.5, rng.integers(8, 25))
    F, df2, p = welch_F([a, b])
    t = stats.ttest_ind(a, b, equal_var=False)
    print(f"{i:4d} {F:10.6f} {t.statistic**2:10.6f} {df2:9.4f} "
          f"{t.df:9.4f} {p:9.6f} {t.pvalue:9.6f}")

            F_W         t²       df2    t 의 df      p(F)      p(t)
   0   8.778756   8.778756   31.1663   31.1663  0.005792  0.005792
   1   1.767110   1.767110   17.8449   17.8449  0.200486  0.200486
   2   0.145205   0.145205   14.9190   14.9190  0.708532  0.708532
   3   2.267246   2.267246   19.6365   19.6365  0.148052  0.148052

import numpy as np

print(f"{'k':>3s} {'df2/(N-k)':>11s}   해석")
for k in [2, 3, 4, 5, 6, 10, 20]:
    r = (k + 1) / (3 * (k - 1))
    note = "표준보다 큼" if r > 1 else ("같음" if abs(r - 1) < 1e-9 else "표준보다 작음")
    print(f"{k:3d} {r:11.4f}   {note}")

print(f"\n{'n':>5s} {'k':>3s} {'df2':>9s} {'N-k':>6s}")
for ns, vs in [([10, 10, 10], [1, 1, 1]),
               ([10, 10, 10], [1, 1, 100]),
               ([5, 5, 100], [1, 1, 1]),
               ([100, 100, 100], [1, 1, 1]),
               ([5, 5, 5], [1, 1, 1])]:
    n = np.array(ns, float)
    v = np.array(vs, float)
    k = len(n)
    w = n / v
    W = w.sum()
    lam = ((1 - w / W)**2 / (n - 1)).sum()
    df2 = (3 / (k**2 - 1) * lam)**-1
    print(f"{str(ns):>16s} {str(vs):>12s} → df2 = {df2:9.4f}   "
          f"N-k = {int(n.sum() - k):4d}")

  k   df2/(N-k)   해석
  2      1.0000   같음
  3      0.6667   표준보다 작음
  4      0.5556   표준보다 작음
  5      0.5000   표준보다 작음
  6      0.4667   표준보다 작음
 10      0.4074   표준보다 작음
 20      0.3684   표준보다 작음

    n   k       df2    N-k
    [10, 10, 10]    [1, 1, 1] → df2 =   18.0000   N-k =   27
    [10, 10, 10]  [1, 1, 100] → df2 =   16.0528   N-k =   27
     [5, 5, 100]    [1, 1, 1] → df2 =    5.8523   N-k =  107
 [100, 100, 100]    [1, 1, 1] → df2 =  198.0000   N-k =  297
       [5, 5, 5]    [1, 1, 1] → df2 =    8.0000   N-k =   12

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

def tukey_any(groups, qcrit):
    """Tukey-Kramer: 합동 MSE 와 공통 자유도를 쓴다."""
    n = np.array([len(g) for g in groups], float)
    m = np.array([g.mean() for g in groups])
    v = np.array([g.var(ddof=1) for g in groups])
    MSE = ((n - 1) * v).sum() / (n.sum() - len(n))
    for i in range(len(n)):
        for j in range(i + 1, len(n)):
            if abs(m[i] - m[j]) / np.sqrt(MSE / 2 * (1 / n[i] + 1 / n[j])) > qcrit:
                return True
    return False

def gh_any(groups, alpha=0.05):
    """게임스-하웰: 쌍마다 웰치형 표준오차와 자유도를 쓴다."""
    k = len(groups)
    n = np.array([len(g) for g in groups], float)
    m = np.array([g.mean() for g in groups])
    v = np.array([g.var(ddof=1) for g in groups])
    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

rng = np.random.default_rng(8642)
M = 1_500
print("완전 귀무 하의 FWER (세 집단 평균 모두 0, 명목 0.05)")
print(f"{'설계':>32s} {'Tukey':>8s} {'G-H':>8s}")
for lab, ns, sds in [
    ("등분산 n=(20,20,20) σ=(1,1,1)", [20, 20, 20], [1, 1, 1]),
    ("이분산 n=(20,20,20) σ=(1,1,4)", [20, 20, 20], [1, 1, 4]),
    ("역페어링 n=(10,20,40) σ=(4,2,1)", [10, 20, 40], [4, 2, 1]),
    ("정페어링 n=(10,20,40) σ=(1,2,4)", [10, 20, 40], [1, 2, 4]),
]:
    qc = qsturng(0.95, 3, sum(ns) - 3)
    a = b = 0
    for _ in range(M):
        gs = [rng.normal(0, s, n) for n, s in zip(ns, sds)]
        a += tukey_any(gs, qc)
        b += gh_any(gs)
    print(f"{lab:>32s} {a / M:8.4f} {b / M:8.4f}")

완전 귀무 하의 FWER (세 집단 평균 모두 0, 명목 0.05)
                              설계    Tukey      G-H
      등분산 n=(20,20,20) σ=(1,1,1)   0.0460   0.0473
      이분산 n=(20,20,20) σ=(1,1,4)   0.0760   0.0487
     역페어링 n=(10,20,40) σ=(4,2,1)   0.2640   0.0373
     정페어링 n=(10,20,40) σ=(1,2,4)   0.0067   0.0393

import numpy as np
from scipy import stats

def welch_p(groups):
    n = np.array([len(g) for g in groups], float)
    m = np.array([g.mean() for g in groups])
    v = np.array([g.var(ddof=1) for g in groups])
    k = len(n)
    w = n / v
    W = w.sum()
    mt = (w * m).sum() / W
    lam = ((1 - w / W)**2 / (n - 1)).sum()
    F = ((w * (m - mt)**2).sum() / (k - 1)
         / (1 + 2 * (k - 2) / (k**2 - 1) * lam))
    return stats.f.sf(F, k - 1, (3 / (k**2 - 1) * lam)**-1)

def trim_stats(x, tr=0.2):
    """절사평균과 윈저화 분산(Yuen 의 표준오차용)."""
    x = np.sort(x)
    n = len(x)
    g = int(np.floor(tr * n))
    win = np.clip(x, x[g], x[n - g - 1])
    d = win.var(ddof=1) * (n - 1) / ((n - 2 * g - 1) * (n - 2 * g))
    return x[g:n - g].mean(), d, n - 2 * g

def yuen_welch(groups, tr=0.2):
    """절사평균에 웰치 구조를 씌운 검정 (Yuen)."""
    st = [trim_stats(g, tr) for g in groups]
    m = np.array([s[0] for s in st])
    d = np.array([s[1] for s in st])
    h = np.array([s[2] for s in st], float)
    k = len(groups)
    w = 1 / d
    W = w.sum()
    mt = (w * m).sum() / W
    lam = ((1 - w / W)**2 / (h - 1)).sum()
    F = ((w * (m - mt)**2).sum() / (k - 1)
         / (1 + 2 * (k - 2) / (k**2 - 1) * lam))
    return stats.f.sf(F, k - 1, (3 / (k**2 - 1) * lam)**-1)

rng = np.random.default_rng(3141)
M = 1_500
print("제1종 오류율 — 세 집단 평균이 모두 0, n=(10,20,40), σ 배율=(3,2,1)")
print(f"{'분포':>24s} {'표준 F':>8s} {'웰치':>8s} {'절사웰치':>9s}")
dists = {
    "정규": lambda n, s: rng.normal(0, s, n),
    "로그정규": lambda n, s: (np.exp(rng.normal(0, 1, n)) - np.exp(0.5)) * s,
    "지수": lambda n, s: (rng.exponential(1, n) - 1) * s,
    "t(3)": lambda n, s: rng.standard_t(3, n) * s,
}
for lab, f in dists.items():
    a = b = c = 0
    for _ in range(M):
        gs = [f(n, s) for n, s in zip([10, 20, 40], [3, 2, 1])]
        a += stats.f_oneway(*gs).pvalue < 0.05
        b += welch_p(gs) < 0.05
        c += yuen_welch(gs) < 0.05
    print(f"{lab:>24s} {a / M:8.4f} {b / M:8.4f} {c / M:9.4f}")

제1종 오류율 — 세 집단 평균이 모두 0, n=(10,20,40), σ 배율=(3,2,1)
                      분포     표준 F       웰치      절사웰치
                      정규   0.2067   0.0427    0.0600
                    로그정규   0.1973   0.1753    0.3120
                      지수   0.2033   0.1140    0.1833
                    t(3)   0.1893   0.0320    0.0533

기저 로그정규의  평균  =  0.000
기저 로그정규의 20% 절사평균 = -0.537    ← 0 이 아니다!

따라서 배율 s 를 곱하면
  s=3: 평균 0,  절사평균 -1.610
  s=2: 평균 0,  절사평균 -1.073
  s=1: 평균 0,  절사평균 -0.537

1. 집단별 n, 평균, 표준편차를 먼저 본다
     ↓
2. 분산비와 n 의 짝짓기 방향을 확인한다
     (역페어링이면 표준 F 를 절대 쓰지 않는다)
     ↓
3. 웰치 F 를 수행한다  →  pg.welch_anova
     ↓
4. 잔차의 치우침을 확인한다
     (심하면 변환·부트스트랩)
     ↓
5. 유의하면 게임스-하웰로 사후검정
     ↓
6. 보고: F(df1, df2) 를 소수점까지, 효과크기, 집단별 요약