일원배치 분산분석: 모형과 가정¶

import numpy as np
from scipy import stats

rng = np.random.default_rng(4242)
print("k = 2 일 때  F = t²")
for n1, n2 in [(10, 10), (12, 8), (30, 25)]:
    x = rng.normal(0, 1, n1)
    y = rng.normal(0.5, 1, n2)
    F, p_F = stats.f_oneway(x, y)
    t, p_t = stats.ttest_ind(x, y, equal_var=True)
    print(f"  n=({n1},{n2}):  F = {F:.6f},  t² = {t**2:.6f},  "
          f"차이 = {abs(F - t**2):.2e}")
    print(f"              p_F = {p_F:.6f},  p_t = {p_t:.6f}")

k = 2 일 때  F = t²
  n=(10,10):  F = 3.611317,  t² = 3.611317,  차이 = 0.00e+00
              p_F = 0.073526,  p_t = 0.073526
  n=(12,8):  F = 6.773407,  t² = 6.773407,  차이 = 2.66e-15
              p_F = 0.018002,  p_t = 0.018002
  n=(30,25):  F = 2.855746,  t² = 2.855746,  차이 = 4.44e-16
              p_F = 0.096923,  p_t = 0.096923

import numpy as np
from scipy import stats
from itertools import combinations

rng = np.random.default_rng(4242)
M, n = 5_000, 15

print(f"{'k':>3s} {'쌍 m':>5s} {'상한 mα':>9s} "
      f"{'독립가정 1-(1-α)^m':>19s} {'실제':>9s}")
for k in [3, 4, 5, 6]:
    m = k * (k - 1) // 2
    hit = 0
    for _ in range(M):
        g = [rng.standard_normal(n) for _ in range(k)]   # 평균이 모두 같다
        hit += any(stats.ttest_ind(g[i], g[j]).pvalue < 0.05
                   for i, j in combinations(range(k), 2))
    print(f"{k:3d} {m:5d} {min(1, 0.05 * m):9.4f} "
          f"{1 - 0.95**m:19.4f} {hit / M:9.4f}")

  k   쌍 m     상한 mα      독립가정 1-(1-α)^m        실제
  3     3    0.1500              0.1426    0.1176
  4     6    0.3000              0.2649    0.2088
  5    10    0.5000              0.4013    0.2940
  6    15    0.7500              0.5367    0.3776

import numpy as np
from scipy import stats

rng = np.random.default_rng(6363)

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

M = 10_000
sigma = np.sqrt([4.2, 18.7, 5.1])
ns = (15, 10, 20)
a = b = 0
for _ in range(M):
    g = [rng.normal(0, s, n) for s, n in zip(sigma, ns)]   # 평균 모두 0
    a += stats.f_oneway(*g).pvalue < 0.05
    b += welch_anova(*g)[1] < 0.05
print(f"σ² = (4.2, 18.7, 5.1),  n = {ns}")
print(f"  고전 F 의 실제 수준   {a / M:.4f}")
print(f"  Welch 의 실제 수준    {b / M:.4f}")
print(f"  분산비 최대/최소 = {max(sigma**2) / min(sigma**2):.2f}")

σ² = (4.2, 18.7, 5.1),  n = (15, 10, 20)
  고전 F 의 실제 수준   0.1138
  Welch 의 실제 수준    0.0512
  분산비 최대/최소 = 4.45

print(f"\n{'n':>16s} {'고전 F':>9s} {'Welch':>9s}")
for ns2 in [(15, 15, 15), (15, 10, 20), (10, 20, 15), (20, 5, 20)]:
    a = b = 0
    for _ in range(M):
        g = [rng.normal(0, s, n) for s, n in zip(sigma, ns2)]
        a += stats.f_oneway(*g).pvalue < 0.05
        b += welch_anova(*g)[1] < 0.05
    print(f"{str(ns2):>16s} {a / M:9.4f} {b / M:9.4f}")

               n      고전 F     Welch
    (15, 15, 15)    0.0642    0.0524
    (15, 10, 20)    0.1096    0.0505
    (10, 20, 15)    0.0313    0.0482
     (20, 5, 20)    0.1789    0.0567

import numpy as np
from scipy import stats

rng = np.random.default_rng(6363)

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

M = 10_000
gens = {"정규": lambda n: rng.standard_normal(n),
        "균등": lambda n: rng.uniform(-1, 1, n),
        "t(3)": lambda n: rng.standard_t(3, n),
        "지수": lambda n: rng.exponential(1, n),
        "로그정규": lambda n: rng.lognormal(0, 1, n),
        "코시": lambda n: rng.standard_cauchy(n)}
print(f"{'분포':>10s} {'고전 F':>9s} {'Welch':>9s} {'크러스컬·월리스':>15s}")
for label, gen in gens.items():
    a = b = c = 0
    for _ in range(M):
        g = [gen(15) for _ in range(3)]      # 세 집단이 같은 분포
        a += stats.f_oneway(*g).pvalue < 0.05
        b += welch_anova(*g)[1] < 0.05
        c += stats.kruskal(*g).pvalue < 0.05
    print(f"{label:>10s} {a / M:9.4f} {b / M:9.4f} {c / M:15.4f}")

        분포      고전 F     Welch        크러스컬·월리스
        정규    0.0522    0.0513          0.0477
        균등    0.0508    0.0535          0.0492
      t(3)    0.0485    0.0416          0.0497
        지수    0.0434    0.0493          0.0460
      로그정규    0.0353    0.0423          0.0472
        코시    0.0178    0.0134          0.0514

print(f"\n검정력 (참 차이 있음)")
print(f"{'분포':>10s} {'고전 F':>9s} {'크러스컬·월리스':>15s}")
shifts = {"정규": lambda n, d: rng.standard_normal(n) + d,
          "지수": lambda n, d: rng.exponential(1, n) + d,
          "로그정규": lambda n, d: rng.lognormal(0, 1, n) + d,
          "코시": lambda n, d: rng.standard_cauchy(n) + d}
for label, gen in shifts.items():
    a = c = 0
    for _ in range(M):
        g = [gen(15, d) for d in (0, 0.8, 1.6)]
        a += stats.f_oneway(*g).pvalue < 0.05
        c += stats.kruskal(*g).pvalue < 0.05
    print(f"{label:>10s} {a / M:9.4f} {c / M:15.4f}")

검정력 (참 차이 있음)
        분포      고전 F        크러스컬·월리스
        정규    0.9752          0.9673
        지수    0.9615          0.9935
      로그정규    0.5607          0.9281
        코시    0.0793          0.4556

import numpy as np
from scipy import stats

rng = np.random.default_rng(9898)
M, n_cluster = 5_000, 6      # 집단마다 군집 6개

print(f"{'군집 크기':>9s} {'ICC':>5s} {'DEFF':>6s} {'고전 F':>9s} "
      f"{'군집평균 F':>11s}")
for m, icc in [(1, 0.0), (5, 0.1), (5, 0.3), (10, 0.3), (10, 0.5)]:
    deff = 1 + (m - 1) * icc
    s_between = np.sqrt(icc)
    s_within = np.sqrt(1 - icc) if icc > 0 else 1.0
    a = c = 0
    for _ in range(M):
        groups, cluster_means = [], []
        for _ in range(3):
            cl = rng.normal(0, s_between, n_cluster) if icc > 0 \
                else np.zeros(n_cluster)
            vals = np.concatenate([cl[j] + rng.normal(0, s_within, m)
                                   for j in range(n_cluster)])
            groups.append(vals)
        a += stats.f_oneway(*groups).pvalue < 0.05
        # 군집 평균을 분석단위로 삼는다
        cm = [np.array([x[j * m:(j + 1) * m].mean()
                        for j in range(n_cluster)]) for x in groups]
        c += stats.f_oneway(*cm).pvalue < 0.05
    print(f"{m:9d} {icc:5.1f} {deff:6.2f} {a / M:9.4f} {c / M:11.4f}")

    군집 크기   ICC   DEFF      고전 F      군집평균 F
        1   0.0   1.00    0.0464      0.0464
        5   0.1   1.40    0.1166      0.0476
        5   0.3   2.20    0.2704      0.0498
       10   0.3   3.70    0.4552      0.0498
       10   0.5   5.50    0.6078      0.0526

① 독립성
    자료가 어떻게 수집되었는가?
    ├─ 군집·반복측정 → 혼합모형 / 군집 단위 분석
    └─ 단순 확률표본 → 계속
          ↓
② 등분산
    집단별 표준편차를 본다
    ├─ 비슷하고 n 도 균형 → 고전 F
    └─ 그 밖 → Welch 분산분석   ← 기본으로 권장
          ↓
③ 정규성
    정규분위수그림, 상자그림
    ├─ 대체로 정규 → 그대로
    ├─ 꼬리가 두꺼움 → 순위 검정 / 절사평균 (검정력 회복)
    └─ 심하게 치우침 → 변환 검토
          ↓
④ 검정 + 효과크기 + 사후비교