Welch 분산분석의 제1종 오류와 검정력 모의실험¶

import numpy as np
import pandas as pd
import pingouin as pg

rng = np.random.default_rng(0)

def simulate_once(null=True):
    ns = [10, 18, 7]
    sigmas = [1.0, 3.0, 6.0]
    means = [10.0, 10.0, 10.0] if null else [10.0, 10.0, 12.0]

    rows = []
    for i, (n, mu, sd) in enumerate(zip(ns, means, sigmas), start=1):
        x = rng.normal(mu, sd, size=n)
        rows += [{"Group": f"G{i}", "Values": v} for v in x]
    df = pd.DataFrame(rows)

    aov = pg.welch_anova(dv="Values", between="Group", data=df)
    # pingouin 0.6부터 열 이름이 "p-unc"에서 "p_unc"로 바뀌었다.
    # 두 이름을 모두 받아들여 버전에 무관하게 동작하도록 한다.
    col = "p_unc" if "p_unc" in aov.columns else "p-unc"
    return float(aov[col].iloc[0])

# 한 번 돌려 형태를 확인한다.
print(f"single run p-value (null) = {simulate_once(null=True):.4f}")

single run p-value (null) = 0.4545

from scipy import stats


def welch_f_by_hand(groups):
    """공식 그대로 F_W 와 두 자유도를 만든다."""
    k = len(groups)
    n = np.array([len(g) for g in groups], float)
    m = np.array([g.mean() for g in groups])
    s2 = np.array([g.var(ddof=1) for g in groups])
    w = n / s2
    W = w.sum()
    y_tilde = (w * m).sum() / W
    Lam = (((1 - w / W) ** 2) / (n - 1)).sum()
    num = (w * (m - y_tilde) ** 2).sum() / (k - 1)
    den = 1 + 2 * (k - 2) / (k ** 2 - 1) * Lam
    return num / den, k - 1, (k ** 2 - 1) / (3 * Lam)


# 세 집단에서 pingouin 과 맞는지 먼저 확인한다 (본문의 rng 은 건드리지 않는다).
chk = np.random.default_rng(99)
g3 = [chk.normal(10, s, n) for n, s in zip((10, 18, 7), (1.0, 3.0, 6.0))]
d3 = pd.DataFrame({"Values": np.concatenate(g3),
                   "Group": sum([[f"G{i+1}"] * len(g) for i, g in enumerate(g3)], [])})
a3 = pg.welch_anova(dv="Values", between="Group", data=d3)
col = "p_unc" if "p_unc" in a3.columns else "p-unc"
F3, d1_3, d2_3 = welch_f_by_hand(g3)
print(f"손계산    F = {F3:.6f},  ddof1 = {d1_3},  ddof2 = {d2_3:.6f}")
print(f"pingouin  F = {a3['F'].iloc[0]:.6f},  ddof1 = {a3['ddof1'].iloc[0]},  "
      f"ddof2 = {a3['ddof2'].iloc[0]:.6f}")

# k=2 이면 보정항이 0 이 되어 F_W = t_W^2 여야 한다.
chk2 = np.random.default_rng(5)
x = chk2.normal(10, 1.0, 12)
y = chk2.normal(11, 4.0, 9)
d2 = pd.DataFrame({"Values": np.r_[x, y], "Group": ["A"] * 12 + ["B"] * 9})
a2 = pg.welch_anova(dv="Values", between="Group", data=d2)
tt = stats.ttest_ind(x, y, equal_var=False)
print(f"\nk=2:  welch_anova  F = {a2['F'].iloc[0]:.6f}, "
      f"ddof2 = {a2['ddof2'].iloc[0]:.6f}, p = {a2[col].iloc[0]:.6f}")
print(f"      ttest_ind    t^2 = {tt.statistic ** 2:.6f}, df = {tt.df:.6f}, "
      f"p = {tt.pvalue:.6f}")
print(f"      보정항 2(k-2)/(k^2-1) = {2 * (2 - 2) / (2 ** 2 - 1)}")

손계산    F = 1.455237,  ddof1 = 2,  ddof2 = 13.231803
pingouin  F = 1.455237,  ddof1 = 2,  ddof2 = 13.231803

k=2:  welch_anova  F = 0.298550, ddof2 = 8.727233, p = 0.598482
      ttest_ind    t^2 = 0.298550, df = 8.727233, p = 0.598482
      보정항 2(k-2)/(k^2-1) = 0.0

def run(n_sims=500, alpha=0.05):
    """귀무가 참인 경우와 거짓인 경우를 함께 돌려 제1종 오류율과 검정력을 잰다.

    분산이 다른 설계에서 보통의 분산분석은 제1종 오류율이 명목수준을 넘는다.
    Welch 는 그것을 0.05 근처로 지켜 준다.
    """
    pvals_null = [simulate_once(null=True) for _ in range(n_sims)]
    pvals_alt  = [simulate_once(null=False) for _ in range(n_sims)]
    type1 = np.mean(np.array(pvals_null) < alpha)
    power = np.mean(np.array(pvals_alt)  < alpha)
    return type1, power

type1, power = run(n_sims=300, alpha=0.05)
print(f"Estimated Type I error: {type1:.3f}")
print(f"Estimated Power:        {power:.3f}")

Estimated Type I error: 0.047
Estimated Power:        0.100

for B in (300, 10000):
    se = np.sqrt(0.05 * 0.95 / B)
    print(f"B = {B:>5}:  제1종 오류 추정의 MC 표준오차 = {se:.4f},  "
          f"95% 폭 = +-{1.96 * se:.4f}")

ns, sigmas = (10, 18, 7), (1.0, 3.0, 6.0)
B = 10000
sim = np.random.default_rng(2026)


def simulate(delta):
    """mu3 를 delta 만큼 올린 설계에서 두 검정의 통계량을 모은다."""
    mus = (10.0, 10.0, 10.0 + delta)
    p_welch = np.empty(B)
    F_classic = np.empty(B)
    for b in range(B):
        gs = [sim.normal(mu, sd, n) for n, mu, sd in zip(ns, mus, sigmas)]
        Fw, d1, d2 = welch_f_by_hand(gs)
        p_welch[b] = stats.f(d1, d2).sf(Fw)
        F_classic[b] = stats.f_oneway(*gs).statistic
    return p_welch, F_classic


pw0, Fc0 = simulate(0.0)
pw1, Fc1 = simulate(2.0)
crit_nominal = stats.f(2, sum(ns) - 3).ppf(0.95)
crit_empirical = np.quantile(Fc0, 0.95)   # 수준을 0.05 로 맞춘 임계값

half = lambda r: 1.96 * np.sqrt(r * (1 - r) / B)
size_w, size_c = np.mean(pw0 < 0.05), np.mean(Fc0 > crit_nominal)
pow_w, pow_c = np.mean(pw1 < 0.05), np.mean(Fc1 > crit_nominal)
pow_c_adj = np.mean(Fc1 > crit_empirical)
print(f"\n제1종 오류율   Welch {size_w:.4f} (+-{half(size_w):.4f})   "
      f"고전 F {size_c:.4f} (+-{half(size_c):.4f})")
print(f"검정력         Welch {pow_w:.4f} (+-{half(pow_w):.4f})   "
      f"고전 F {pow_c:.4f} (+-{half(pow_c):.4f})")
print(f"\n고전 F 의 임계값:  이론 {crit_nominal:.4f}  vs  수준 0.05 를 맞춘 "
      f"경험적 임계값 {crit_empirical:.4f}")
print(f"수준을 맞춘 뒤의 고전 F 검정력 = {pow_c_adj:.4f} (+-{half(pow_c_adj):.4f})")

B =   300:  제1종 오류 추정의 MC 표준오차 = 0.0126,  95% 폭 = +-0.0247
B = 10000:  제1종 오류 추정의 MC 표준오차 = 0.0022,  95% 폭 = +-0.0043

제1종 오류율   Welch 0.0532 (+-0.0044)   고전 F 0.1418 (+-0.0068)
검정력         Welch 0.1059 (+-0.0060)   고전 F 0.2659 (+-0.0087)

고전 F 의 임계값:  이론 3.2945  vs  수준 0.05 를 맞춘 경험적 임계값 5.9757
수준을 맞춘 뒤의 고전 F 검정력 = 0.1222 (+-0.0064)

import numpy as np
from scipy import stats

def welch_p(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])
    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)

NS, SD = [10, 18, 7], [1.0, 3.0, 6.0]

def draw(rng, mus):
    return [rng.normal(u, s, n) for n, s, u in zip(NS, SD, mus)]

def wilson(x, n, z=1.96):
    """이항 비율의 윌슨 구간. p̂ 이 0 이나 1 에 가까워도 무너지지 않는다."""
    p = x / n
    d = 1 + z * z / n
    c = p + z * z / (2 * n)
    h = z * np.sqrt(p * (1 - p) / n + z * z / (4 * n * n))
    return (c - h) / d, (c + h) / d

rng = np.random.default_rng(11)
print(f"{'B':>7s} {'추정 α':>8s} {'95% 윌슨 구간':>22s} {'폭':>7s}")
for B in [300, 3_000, 30_000]:
    x = sum(welch_p(draw(rng, [10, 10, 10])) < 0.05 for _ in range(B))
    lo, hi = wilson(x, B)
    print(f"{B:7d} {x / B:8.4f}   [{lo:.4f}, {hi:.4f}] {hi - lo:8.4f}")

for tol in [0.01, 0.005, 0.002]:
    B = int(np.ceil(1.96**2 * 0.05 * 0.95 / tol**2))
    print(f"  반폭 {tol} 을 얻으려면 B ≈ {B:,}")

      B     추정 α              95% 윌슨 구간       폭
    300   0.0900   [0.0626, 0.1278]   0.0652
   3000   0.0493   [0.0421, 0.0577]   0.0155
  30000   0.0525   [0.0500, 0.0551]   0.0050
  반폭 0.01 을 얻으려면 B ≈ 1,825
  반폭 0.005 을 얻으려면 B ≈ 7,300
  반폭 0.002 을 얻으려면 B ≈ 45,619

import numpy as np
from scipy import stats

def welch_p(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])
    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)

NS, SD = [10, 18, 7], [1.0, 3.0, 6.0]

def draw(rng, mus):
    return [rng.normal(u, s, n) for n, s, u in zip(NS, SD, mus)]

rng = np.random.default_rng(2024)
B = 20_000
c0 = w0 = c1 = w1 = 0
for _ in range(B):
    gs = draw(rng, [10, 10, 10])            # 귀무
    c0 += stats.f_oneway(*gs).pvalue < 0.05
    w0 += welch_p(gs) < 0.05
    gs = draw(rng, [10, 10, 12])            # 대립
    c1 += stats.f_oneway(*gs).pvalue < 0.05
    w1 += welch_p(gs) < 0.05

print(f"{'':10s} {'제1종 오류':>10s} {'검정력':>9s}")
print(f"{'고전 F':>10s} {c0 / B:10.4f} {c1 / B:9.4f}")
print(f"{'웰치':>10s} {w0 / B:10.4f} {w1 / B:9.4f}")

               제1종 오류       검정력
      고전 F     0.1419    0.2700
        웰치     0.0533    0.1091

import numpy as np
from scipy import stats

def welch_p(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])
    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)

SD = [1.0, 3.0, 6.0]
rng = np.random.default_rng(707)
B = 6_000
allocs = [
    ("본문 (10,18,7)", [10, 18, 7]),
    ("균등 (12,12,11)", [12, 12, 11]),
    ("σ 비례 (4,10,21)", [4, 10, 21]),
    ("σ² 비례 (1,7,27)", [2, 7, 26]),
    ("작은집단 보강 (7,10,18)", [7, 10, 18]),
]
print(f"{'배분 (총 N=35)':>22s} {'제1종 오류':>10s} {'검정력':>9s}")
for lab, ns in allocs:
    a = b = 0
    for _ in range(B):
        gs = [rng.normal(10, s, n) for n, s in zip(ns, SD)]
        a += welch_p(gs) < 0.05
        gs = [rng.normal(u, s, n)
              for n, s, u in zip(ns, SD, [10, 10, 12])]
        b += welch_p(gs) < 0.05
    print(f"{lab:>22s} {a / B:10.4f} {b / B:9.4f}")

           배분 (총 N=35)     제1종 오류       검정력
          본문 (10,18,7)     0.0520    0.1103
         균등 (12,12,11)     0.0503    0.1492
        σ 비례 (4,10,21)     0.0480    0.1945
        σ² 비례 (1,7,27)     0.0485    0.1708
     작은집단 보강 (7,10,18)     0.0502    0.1795

import numpy as np
from scipy import stats

def welch_p(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])
    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)

NS, SD = [10, 18, 7], [1.0, 3.0, 6.0]
rng = np.random.default_rng(31337)
R, B = 60, 800
d_crn, d_ind = [], []
for _ in range(R):
    a = b = 0
    for _ in range(B):
        gs = [rng.normal(10, s, n) for n, s in zip(NS, SD)]
        a += stats.f_oneway(*gs).pvalue < 0.05
        b += welch_p(gs) < 0.05                 # 같은 자료에 두 검정
    d_crn.append((a - b) / B)

    a = b = 0
    for _ in range(B):
        gs = [rng.normal(10, s, n) for n, s in zip(NS, SD)]
        a += stats.f_oneway(*gs).pvalue < 0.05
        gs = [rng.normal(10, s, n) for n, s in zip(NS, SD)]
        b += welch_p(gs) < 0.05                 # 각각 새 자료
    d_ind.append((a - b) / B)

d_crn, d_ind = np.array(d_crn), np.array(d_ind)
print(f"차이(고전 - 웰치) 추정, R={R} 회 반복, 회당 B={B}")
print(f"  공통난수:  평균 {d_crn.mean():.4f},  "
      f"표준편차 {d_crn.std(ddof=1):.5f}")
print(f"  독립난수:  평균 {d_ind.mean():.4f},  "
      f"표준편차 {d_ind.std(ddof=1):.5f}")
r = d_ind.var(ddof=1) / d_crn.var(ddof=1)
print(f"  분산비 = {r:.2f} 배")
print(f"  → 같은 정밀도를 얻는 데 필요한 B 가 {r:.1f} 배 차이")

차이(고전 - 웰치) 추정, R=60 회 반복, 회당 B=800
  공통난수:  평균 0.0882,  표준편차 0.01006
  독립난수:  평균 0.0874,  표준편차 0.01341
  분산비 = 1.78 배
  → 같은 정밀도를 얻는 데 필요한 B 가 1.8 배 차이

제1종 오류율 (명목 0.05, B = 20,000, 씨앗 2024)

  설계                  고전 F            웰치
  σ=(1,3,6), n=(10,18,7)  0.142 [0.137,0.147]  0.053 [0.050,0.056]