분산분석의 가정¶

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)

rng = np.random.default_rng(1001)
B = 20_000

print("균형 설계 n = (10, 10, 10)")
for lab, sds in [("등분산 σ²=(1,1,1)", [1, 1, 1]),
                 ("본문 σ²=(1,1,9)", [1, 1, 3])]:
    a = b = 0
    for _ in range(B):
        gs = [rng.normal(0, s, 10) for s in sds]
        a += stats.f_oneway(*gs).pvalue < 0.05
        b += welch_p(gs) < 0.05
    print(f"  {lab:20s} 표준 F {a / B:.4f}   웰치 {b / B:.4f}")

print("\n불균형 n = (20, 20, 5), σ²=(1,1,9)")
for lab, ns, sds in [("역페어링", [20, 20, 5], [1, 1, 3]),
                     ("정페어링", [5, 20, 20], [1, 1, 3])]:
    a = b = 0
    for _ in range(B):
        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:20s} 표준 F {a / B:.4f}   웰치 {b / B:.4f}")

균형 설계 n = (10, 10, 10)
  등분산 σ²=(1,1,1)       표준 F 0.0498   웰치 0.0472
  본문 σ²=(1,1,9)        표준 F 0.0801   웰치 0.0534

불균형 n = (20, 20, 5), σ²=(1,1,9)
  역페어링                 표준 F 0.2858   웰치 0.0587
  정페어링                 표준 F 0.0328   웰치 0.0568

import numpy as np
from scipy import stats

rng = np.random.default_rng(606)
B = 8_000

def ar1(n, rho, rng):
    """1차 자기상관을 가진 계열. 주변분산은 1 로 맞춘다."""
    e = rng.normal(0, 1, n)
    x = np.empty(n)
    x[0] = e[0] / np.sqrt(1 - rho**2)
    for t in range(1, n):
        x[t] = rho * x[t - 1] + e[t]
    return x

cases = {
    "기준 (모든 가정 성립)": lambda: [rng.normal(0, 1, 15) for _ in range(3)],
    "등분산 위반 (σ=1,1,3)": lambda: [rng.normal(0, s, 15) for s in [1, 1, 3]],
    "등분산 위반 (σ=1,3,5)": lambda: [rng.normal(0, s, 15) for s in [1, 3, 5]],
    "정규성 위반 (지수)": lambda: [rng.exponential(1, 15) - 1 for _ in range(3)],
    "정규성 위반 (로그정규)":
        lambda: [np.exp(rng.normal(0, 1, 15)) - np.exp(.5) for _ in range(3)],
    "독립성 위반 (AR1 ρ=0.3)": lambda: [ar1(15, 0.3, rng) for _ in range(3)],
    "독립성 위반 (AR1 ρ=0.6)": lambda: [ar1(15, 0.6, rng) for _ in range(3)],
}
print("k=3, 집단당 n=15, 균형, 명목 0.05")
print(f"{'상황':>26s} {'표준 F 오류율':>12s}")
for lab, f in cases.items():
    a = sum(stats.f_oneway(*f()).pvalue < 0.05 for _ in range(B))
    print(f"{lab:>26s} {a / B:12.4f}")

k=3, 집단당 n=15, 균형, 명목 0.05
                        상황     표준 F 오류율
             기준 (모든 가정 성립)       0.0548
          등분산 위반 (σ=1,1,3)       0.0720
          등분산 위반 (σ=1,3,5)       0.0670
               정규성 위반 (지수)       0.0408
             정규성 위반 (로그정규)       0.0365
        독립성 위반 (AR1 ρ=0.3)       0.1963
        독립성 위반 (AR1 ρ=0.6)       0.4754

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)

rng = np.random.default_rng(717)
B = 10_000
print(f"{'상황':>24s} {'항상 F':>8s} {'항상 웰치':>9s} "
      f"{'2단계':>8s} {'레빈 기각률':>10s}")
for lab, ns, sds in [
    ("등분산 n=(15,15,15)", [15, 15, 15], [1, 1, 1]),
    ("이분산 n=(15,15,15) σ=(1,1,3)", [15, 15, 15], [1, 1, 3]),
    ("역페어링 n=(20,20,6) σ=(1,1,4)", [20, 20, 6], [1, 1, 4]),
]:
    a = b = c = lev = 0
    for _ in range(B):
        gs = [rng.normal(0, s, n) for n, s in zip(ns, sds)]
        pf = stats.f_oneway(*gs).pvalue
        pw = welch_p(gs)
        pl = stats.levene(*gs, center="median").pvalue
        a += pf < 0.05
        b += pw < 0.05
        lev += pl < 0.05
        c += (pw < 0.05) if pl < 0.05 else (pf < 0.05)   # 2단계 절차
    print(f"{lab:>24s} {a / B:8.4f} {b / B:9.4f} {c / B:8.4f} {lev / B:10.4f}")

                      상황     항상 F     항상 웰치      2단계     레빈 기각률
        등분산 n=(15,15,15)   0.0493    0.0487   0.0495     0.0263
이분산 n=(15,15,15) σ=(1,1,3)   0.0734    0.0471   0.0497     0.9373
역페어링 n=(20,20,6) σ=(1,1,4)   0.3142    0.0551   0.0842     0.8926

import numpy as np
from scipy import stats

rng = np.random.default_rng(818)
B = 8_000
dists = {
    "정규 (치우침 0)": lambda n: rng.normal(0, 1, n),
    "지수 (치우침 2)": lambda n: rng.exponential(1, n) - 1,
    "카이제곱 1 (치우침 2.8)": lambda n: rng.chisquare(1, n) - 1,
    "로그정규 (치우침 6.2)": lambda n: np.exp(rng.normal(0, 1, n)) - np.exp(.5),
    "파레토 a=3 (치우침 ∞)": lambda n: (rng.pareto(3, n) + 1) - 1.5,
}
ns = [5, 10, 20, 50, 100]
print("k=3, 등분산·균형, 명목 0.05")
print(f"{'분포':>22s} " + " ".join(f"{'n=' + str(n):>8s}" for n in ns))
for lab, f in dists.items():
    row = [sum(stats.f_oneway(f(n), f(n), f(n)).pvalue < 0.05
               for _ in range(B)) / B for n in ns]
    print(f"{lab:>22s} " + " ".join(f"{r:8.4f}" for r in row))

k=3, 등분산·균형, 명목 0.05
                    분포      n=5     n=10     n=20     n=50    n=100
            정규 (치우침 0)   0.0522   0.0476   0.0469   0.0471   0.0535
            지수 (치우침 2)   0.0393   0.0429   0.0508   0.0480   0.0530
      카이제곱 1 (치우침 2.8)   0.0346   0.0334   0.0447   0.0454   0.0418
        로그정규 (치우침 6.2)   0.0314   0.0341   0.0385   0.0420   0.0446
       파레토 a=3 (치우침 ∞)   0.0301   0.0336   0.0357   0.0410   0.0401

rng = np.random.default_rng(8181)
B = 8_000
lg = lambda n: np.exp(rng.normal(0, 1, n)) - np.exp(.5)
print("로그정규 기저, 척도를 곱해 이분산을 만듦 (명목 0.05)")
print(f"{'설계':>32s} {'표준 F':>8s}")
for lab, ns, sc in [
    ("균형·등척도 n=(20,20,20)", [20, 20, 20], [1, 1, 1]),
    ("불균형 n=(40,20,8) 등척도", [40, 20, 8], [1, 1, 1]),
    ("역페어링 n=(40,20,8) 척도=(1,2,4)", [40, 20, 8], [1, 2, 4]),
    ("정페어링 n=(8,20,40) 척도=(1,2,4)", [8, 20, 40], [1, 2, 4]),
]:
    a = sum(stats.f_oneway(*[lg(n) * s for n, s in zip(ns, sc)]).pvalue < 0.05
            for _ in range(B))
    print(f"{lab:>32s} {a / B:8.4f}")

로그정규 기저, 척도를 곱해 이분산을 만듦 (명목 0.05)
                              설계     표준 F
             균형·등척도 n=(20,20,20)   0.0356
             불균형 n=(40,20,8) 등척도   0.0489
     역페어링 n=(40,20,8) 척도=(1,2,4)   0.3003
     정페어링 n=(8,20,40) 척도=(1,2,4)   0.0333

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(919)
B = 6_000
print("(가) 자료가 정말 정규일 때 샤피로 검정의 기각률 (명목 0.05)")
print(f"{'집단당 n':>10s} {'참 오차에 적용':>14s} {'잔차에 적용':>12s}")
for n in [10, 20, 50]:
    a = b = 0
    for _ in range(B):
        gs = [rng.normal(0, 1, n) for _ in range(3)]
        eps = np.concatenate(gs)                          # 참 오차
        res = np.concatenate([g - g.mean() for g in gs])  # 잔차
        a += stats.shapiro(eps).pvalue < 0.05
        b += stats.shapiro(res).pvalue < 0.05
    print(f"{n:10d} {a / B:14.4f} {b / B:12.4f}")

print("\n(나) 약간만 치우친 분포(감마 shape=9, 치우침 0.67)에서")
print(f"{'집단당 n':>10s} {'샤피로 기각률':>13s} {'F 의 실제 오류율':>15s}")
for n in [10, 20, 50, 100, 200]:
    a = b = 0
    for _ in range(B // 2):
        gs = [rng.gamma(9, 1, n) - 9 for _ in range(3)]
        res = np.concatenate([g - g.mean() for g in gs])
        a += stats.shapiro(res).pvalue < 0.05
        b += stats.f_oneway(*gs).pvalue < 0.05
    print(f"{n:10d} {a / (B // 2):13.4f} {b / (B // 2):15.4f}")

(가) 자료가 정말 정규일 때 샤피로 검정의 기각률 (명목 0.05)
     집단당 n       참 오차에 적용       잔차에 적용
        10         0.0565       0.0495
        20         0.0503       0.0492
        50         0.0493       0.0493

(나) 약간만 치우친 분포(감마 shape=9, 치우침 0.67)에서
     집단당 n       샤피로 기각률      F 의 실제 오류율
        10        0.1867          0.0503
        20        0.3837          0.0470
        50        0.8233          0.0583
       100        0.9903          0.0533
       200        1.0000          0.0533

0. 설계를 검토한다 ── 독립성은 여기서 끝난다
     │  (반복측정·군집이면 분산분석을 쓰지 않는다)
     ↓
1. 칸별 n, 평균, 표준편차를 본다
     │  불균형인가? 분산비는? 역페어링인가?
     ↓
2. 처음부터 웰치를 쓴다 ── 레빈으로 분기하지 않는다
     ↓
3. 잔차를 그림으로 본다 ── Q-Q, 잔차 대 적합값
     │  치우침이 심하면 변환이나 부트스트랩
     ↓
4. 공변량이 있으면 선형성 확인
     ↓
5. 보고: 방법의 선택 이유와 가정 점검 결과를 함께