등분산성 확인¶

import numpy as np
import pandas as pd
from statsmodels.formula.api import ols

# 이 페이지의 모든 진단은 아래 모형 하나를 놓고 수행한다.
# 집단마다 표준편차를 1.0, 1.3, 1.6으로 다르게 주어 진단이 무엇을 잡아내는지
# (그리고 무엇을 못 잡아내는지) 볼 수 있게 했다.
rng = np.random.default_rng(42)
n = 20
data = pd.DataFrame({
    "group": np.repeat(["A", "B", "C"], n),
    "response": np.concatenate([
        rng.normal(10.0, 1.0, n),
        rng.normal(10.8, 1.3, n),
        rng.normal(12.0, 1.6, n),
    ]),
})
model = ols("response ~ C(group)", data=data).fit()

print(data.groupby("group").response.agg(["count", "mean", "std"]).round(3))
print(f"\nF = {model.fvalue:.4f}, p = {model.f_pvalue:.4f}")

e = model.resid
print(f"\n잔차 전체의 합 = {e.sum():+.3e}")
print("집단별 잔차 합")
for g, v in e.groupby(data["group"]):
    print(f"  {g}: {v.sum():+.3e}")

# 설계행렬의 세 열(절편, B 더미, C 더미)에 잔차를 직접 내적해 본다.
X = model.model.exog
print("\n설계행렬 열이름:", model.model.exog_names)
print("X^T e        =", np.array2string(X.T @ e.values, precision=3))

# 이 쪽에서 쓸 양: 집단별 잔차 분산과 그 합동값
s2 = e.groupby(data["group"]).var(ddof=1)
print("\n집단별 잔차 분산 s_g^2")
print(s2.round(6).to_string())
print(f"\n합동분산 (s_A^2+s_B^2+s_C^2)/3 = {s2.mean():.10f}")
print(f"model.mse_resid               = {model.mse_resid:.10f}")
print(f"분산의 최대/최소 = {s2.max() / s2.min():.4f},  "
      f"표준편차의 최대/최소 = {np.sqrt(s2.max() / s2.min()):.4f}")

       count    mean    std
group                      
A         20   9.967  0.870
B         20  10.942  1.034
C         20  12.191  1.145

F = 23.7708, p = 0.0000

잔차 전체의 합 = -2.469e-13
집단별 잔차 합
  A: -7.105e-14
  B: -5.329e-15
  C: -1.705e-13

설계행렬 열이름: ['Intercept', 'C(group)[T.B]', 'C(group)[T.C]']
X^T e        = [-2.469e-13 -5.329e-15 -1.705e-13]

집단별 잔차 분산 s_g^2
group
A    0.757209
B    1.068200
C    1.311068

합동분산 (s_A^2+s_B^2+s_C^2)/3 = 1.0454926260
model.mse_resid               = 1.0454926260
분산의 최대/최소 = 1.7314,  표준편차의 최대/최소 = 1.3158

import numpy as np
from scipy.stats import f_oneway, levene

group1 = data[data['group'] == 'A']['response']
group2 = data[data['group'] == 'B']['response']
group3 = data[data['group'] == 'C']['response']

# scipy의 기본값은 center='median'(Brown-Forsythe 변형)이다.
# 원래의 Levene(1960)은 평균을 쓰지만 이상점에 약해 기본값이 중앙값으로 바뀌었다.
stat, p_value = levene(group1, group2, group3)
print(f"Levene's Test: F = {stat:.4f}, p-value = {p_value:.4f}")

# Levene 검정은 z = |y - 집단중위수| 에 돌린 일원배치 F-검정과 같은 것이다.
zs = [np.abs(g - np.median(g)) for g in (group1, group2, group3)]
F_hand, p_hand = f_oneway(*zs)
print(f"|y - 중위수| 의 F-검정: F = {F_hand:.4f}, p-value = {p_hand:.4f}")
print(f"두 F 의 차이: {abs(F_hand - stat):.3e}")

print(f"\n{'z 평균':>8s} {'0.7979 s':>9s} {'s':>8s}  집단")
for name, g, z in zip("ABC", (group1, group2, group3), zs):
    s = g.std(ddof=1)
    print(f"{z.mean():8.4f} {np.sqrt(2 / np.pi) * s:9.4f} {s:8.4f}  {name}")

Levene's Test: F = 0.4091, p-value = 0.6662
|y - 중위수| 의 F-검정: F = 0.4091, p-value = 0.6662
두 F 의 차이: 5.551e-17

    z 평균  0.7979 s        s  집단
  0.6952    0.6943   0.8702  A
  0.8180    0.8246   1.0335  B
  0.8773    0.9136   1.1450  C

import numpy as np
from scipy import stats

# 등분산이 **참**인 자료를 만들어 기각률을 잰다. 이론값은 명목수준 0.05 다.
rng = np.random.default_rng(2026)
k, n, B = 3, 30, 2000
mc = np.sqrt(0.05 * 0.95 / B)      # 몬테카를로 오차

cases = [
    ("정규 (대칭)", lambda: rng.normal(0, 1, n)),
    ("지수 (치우침)", lambda: rng.exponential(1.0, n)),
    ("로그정규 (심한 치우침)", lambda: rng.lognormal(0, 1, n)),
]

print(f"집단 {k}개, 각 n = {n}, 반복 {B}회, 명목수준 0.05, "
      f"몬테카를로 오차 {mc:.4f}")
print(f"\n{'median':>8s} {'mean':>8s} {'왜도':>6s}  자료 분포")
for name, draw in cases:
    cm = ca = 0
    skew = []
    for _ in range(B):
        gs = [draw() for _ in range(k)]
        cm += stats.levene(*gs, center="median").pvalue < 0.05
        ca += stats.levene(*gs, center="mean").pvalue < 0.05
        skew.append(stats.skew(gs[0]))
    print(f"{cm / B:8.4f} {ca / B:8.4f} {np.mean(skew):6.2f}  {name}")

집단 3개, 각 n = 30, 반복 2000회, 명목수준 0.05, 몬테카를로 오차 0.0049

  median     mean     왜도  자료 분포
  0.0370   0.0460  -0.01  정규 (대칭)
  0.0465   0.1765   1.47  지수 (치우침)
  0.0340   0.2505   2.19  로그정규 (심한 치우침)

import numpy as np
from scipy.stats import bartlett, chi2

# Bartlett 검정은 Levene 보다 검정력이 높지만 정규성을 전제한다.
# 자료가 정규에서 조금만 벗어나도 등분산을 지나치게 자주 기각한다.
stat, p_value = bartlett(group1, group2, group3)
print(f"Bartlett's Test: chi2 = {stat:.4f}, p-value = {p_value:.4f}")

# 닫힌 꼴을 손으로 계산해 맞춰 본다.
gs = [group1.values, group2.values, group3.values]
ns = np.array([len(g) for g in gs])
s2 = np.array([g.var(ddof=1) for g in gs])
k, N = len(gs), ns.sum()
sp2 = ((ns - 1) * s2).sum() / (N - k)
num = (N - k) * np.log(sp2) - ((ns - 1) * np.log(s2)).sum()
corr = 1 + (1 / (3 * (k - 1))) * ((1 / (ns - 1)).sum() - 1 / (N - k))
T = num / corr
print(f"\ns_g^2        = {np.array2string(s2, precision=6)}")
print(f"s_p^2        = {sp2:.10f}   (model.mse_resid = {model.mse_resid:.10f})")
print(f"분자 M       = {num:.10f}")
print(f"보정계수 C   = {corr:.10f}")
print(f"T = M / C    = {T:.10f}")
print(f"p = P(chi2_2 > T) = {chi2.sf(T, k - 1):.10f}")
print(f"scipy 와의 차이: {abs(T - stat):.3e}")

Bartlett's Test: chi2 = 1.3880, p-value = 0.4996

s_g^2        = [0.757209 1.0682   1.311068]
s_p^2        = 1.0454926260   (model.mse_resid = 1.0454926260)
분자 M       = 1.4204884460
보정계수 C   = 1.0233918129
T = M / C    = 1.3880201386
p = P(chi2_2 > T) = 0.4995687417
scipy 와의 차이: 0.000e+00

import numpy as np
from scipy import stats

# 분산이 **모두 같은** 자료를 만들어 기각률을 잰다. 이론값은 0.05 다.
rng = np.random.default_rng(2026)
k, n, B = 3, 30, 2000
mc = np.sqrt(0.05 * 0.95 / B)

# t_3 는 분산이 3/(3-2) = 3 이다. 세 집단에 같은 분포를 주었으니
# 등분산 귀무가설은 참이다. 다른 것은 꼬리 두께뿐이다.
cases = [
    ("정규 — Bartlett 의 가정이 성립", lambda: rng.normal(0, 1, n)),
    ("t_3 — 대칭이지만 꼬리가 두껍다", lambda: rng.standard_t(3, n)),
]

print(f"집단 {k}개, 각 n = {n}, 반복 {B}회, 명목수준 0.05, "
      f"몬테카를로 오차 {mc:.4f}")
print(f"\n{'Bartlett':>9s} {'Levene':>8s} {'초과첨도':>8s}  자료 분포")
for name, draw in cases:
    cb = cl = 0
    kurt = []
    for _ in range(B):
        gs = [draw() for _ in range(k)]
        cb += stats.bartlett(*gs).pvalue < 0.05
        cl += stats.levene(*gs).pvalue < 0.05
        kurt.append(stats.kurtosis(gs[0]))
    print(f"{cb / B:9.4f} {cl / B:8.4f} {np.mean(kurt):8.2f}  {name}")

집단 3개, 각 n = 30, 반복 2000회, 명목수준 0.05, 몬테카를로 오차 0.0049

 Bartlett   Levene     초과첨도  자료 분포
   0.0435   0.0370    -0.22  정규 — Bartlett 의 가정이 성립
   0.4720   0.0440     2.60  t_3 — 대칭이지만 꼬리가 두껍다

import matplotlib.pyplot as plt
import pandas as pd

# 그림에서 읽으려는 것을 먼저 수로 적어 둔다.
e, g = model.resid, data["group"]
tab = pd.DataFrame({
    "적합값": model.fittedvalues.groupby(g).first().round(4),
    "잔차 s": e.groupby(g).std(ddof=1).round(4),
    "최소": e.groupby(g).min().round(4),
    "최대": e.groupby(g).max().round(4),
})
tab["범위"] = (tab["최대"] - tab["최소"]).round(4)
print(tab.to_string())

sd = e.groupby(g).std(ddof=1)
print(f"\n서로 다른 적합값이 {model.fittedvalues.nunique()}개뿐이다 (세로 띠 세 줄)")
print(f"잔차 s 의 최대/최소 = {sd.max() / sd.min():.4f}   (모표준편차의 참 비는 1.6)")
print(f"범위의 최대/최소    = {tab['범위'].max() / tab['범위'].min():.4f}")

# 일원배치 분산분석에서 적합값은 집단평균뿐이므로 세로줄이 집단 수만큼만 생긴다.
# 회귀분석의 잔차 그림처럼 연속적으로 퍼지지 않는다.
plt.scatter(model.fittedvalues, model.resid, alpha=0.6)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel("Fitted Values")
plt.ylabel("Residuals")
plt.title("Residuals vs. Fitted Values")
plt.show()

           적합값    잔차 s      최소      최대      범위
group                                         
A       9.9671  0.8702 -1.9181  1.1602  3.0783
B      10.9423  1.0335 -1.2345  2.6418  3.8763
C      12.1909  1.1450 -2.5224  2.2010  4.7234

서로 다른 적합값이 3개뿐이다 (세로 띠 세 줄)
잔차 s 의 최대/최소 = 1.3158   (모표준편차의 참 비는 1.6)
범위의 최대/최소    = 1.5344

import numpy as np

rng = np.random.default_rng(808)

def report(name, groups, tfs):
    print(f"[{name}]")
    print(f"{'변환':>12s} "
          + " ".join(f"{'g' + str(i + 1):>9s}" for i in range(len(groups)))
          + f" {'최대/최소':>9s}")
    for tn, tf in tfs:
        sds = [tf(g).std(ddof=1) for g in groups]
        print(f"{tn:>12s} " + " ".join(f"{s:9.4f}" for s in sds)
              + f" {max(sds) / min(sds):9.2f}")
    print()

mus = [5.0, 20.0, 80.0]

gs = [rng.poisson(m, 4_000).astype(float) for m in mus]
report("포아송  σ² = μ", gs,
       [("원 척도", lambda x: x), ("√x", np.sqrt),
        ("log(x+1)", lambda x: np.log(x + 1))])

gs = [rng.gamma(4, m / 4, 4_000) for m in mus]
report("감마  σ ∝ μ (CV 일정)", gs,
       [("원 척도", lambda x: x), ("√x", np.sqrt), ("log x", np.log)])

ps = [0.05, 0.3, 0.5]
gs = [rng.binomial(50, p, 4_000) / 50 for p in ps]
report("이항 비율  σ² = p(1-p)/n", gs,
       [("원 척도", lambda x: x),
        ("arcsin √p", lambda x: np.arcsin(np.sqrt(x))),
        ("logit", lambda x: np.log((x + 0.01) / (1 - x + 0.01)))])

[포아송  σ² = μ]
          변환        g1        g2        g3     최대/최소
        원 척도    2.2708    4.5119    8.9901      3.96
          √x    0.5476    0.5080    0.5038      1.09
    log(x+1)    0.4257    0.2213    0.1120      3.80

[감마  σ ∝ μ (CV 일정)]
          변환        g1        g2        g3     최대/최소
        원 척도    2.4679   10.4018   40.1318     16.26
          √x    0.5435    1.1302    2.1982      4.04
       log x    0.5253    0.5413    0.5304      1.03

[이항 비율  σ² = p(1-p)/n]
          변환        g1        g2        g3     최대/최소
        원 척도    0.0309    0.0652    0.0713      2.30
   arcsin √p    0.0852    0.0722    0.0720      1.18
       logit    0.6669    0.3128    0.2851      2.34

import numpy as np

rng = np.random.default_rng(808)
print("자료를 σ ∝ μ^θ 로 만들고, log(IQR) 대 log(중앙값)의 기울기를 잰다")
print(f"{'참 θ':>6s} {'기울기 b':>9s} {'추정 지수 1-b':>12s} {'권장 변환':>12s}")
for theta in [0.0, 0.5, 1.0, 1.5]:
    mus = np.array([5.0, 10.0, 20.0, 40.0, 80.0])
    med, iqr = [], []
    for m in mus:
        x = rng.normal(m, 0.3 * m**theta, 3_000)
        med.append(np.median(x))
        iqr.append(np.subtract(*np.percentile(x, [75, 25])))
    b = np.polyfit(np.log(med), np.log(iqr), 1)[0]
    p = 1 - b
    rec = ("로그" if abs(p) < 0.25 else
           "제곱근" if abs(p - 0.5) < 0.25 else
           "없음" if abs(p - 1) < 0.25 else
           "역수" if abs(p + 1) < 0.3 else f"x^{p:.2f}")
    print(f"{theta:6.1f} {b:9.4f} {p:12.4f} {rec:>12s}")

자료를 σ ∝ μ^θ 로 만들고, log(IQR) 대 log(중앙값)의 기울기를 잰다
   참 θ     기울기 b    추정 지수 1-b        권장 변환
   0.0    0.0037       0.9963           없음
   0.5    0.5215       0.4785          제곱근
   1.0    1.0087      -0.0087           로그
   1.5    1.5161      -0.5161      x^-0.52

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)

def perm_F_p(gs, rng, R=400):
    obs = stats.f_oneway(*gs).statistic
    z = np.concatenate(gs)
    ns = [len(g) for g in gs]
    cnt = 0
    for _ in range(R):
        p = rng.permutation(z)
        idx, parts = 0, []
        for n in ns:
            parts.append(p[idx:idx + n])
            idx += n
        cnt += stats.f_oneway(*parts).statistic >= obs
    return (cnt + 1) / (R + 1)

rng = np.random.default_rng(909)
NS = [10, 20, 40]
B = 2_500

def draw(mus, shape=4.0):
    return [rng.gamma(shape, m / shape, n) for m, n in zip(mus, NS)]

print("감마 자료(CV 일정) — 평균이 커지면 분산도 커진다, n=(10,20,40)")
print(f"{'상황':>22s} {'표준 F':>8s} {'웰치':>8s} {'로그변환 F':>10s} {'순열 F':>8s}")
for lab, mus in [("귀무: μ=(20,20,20)", [20.0, 20.0, 20.0]),
                 ("대립: μ=(20,26,32)", [20.0, 26.0, 32.0])]:
    a = b = c = d = 0
    for _ in range(B):
        gs = draw(mus)
        a += stats.f_oneway(*gs).pvalue < 0.05
        b += welch_p(gs) < 0.05
        c += stats.f_oneway(*[np.log(g) for g in gs]).pvalue < 0.05
        d += perm_F_p(gs, rng) < 0.05
    print(f"{lab:>22s} {a / B:8.4f} {b / B:8.4f} {c / B:10.4f} {d / B:8.4f}")

감마 자료(CV 일정) — 평균이 커지면 분산도 커진다, n=(10,20,40)
                    상황     표준 F       웰치     로그변환 F     순열 F
      귀무: μ=(20,20,20)   0.0524   0.0604     0.0532   0.0532
      대립: μ=(20,26,32)   0.5868   0.7156     0.6068   0.5860

이분산이 확인되었다
    │
    ├─ 산포가 평균과 함께 커지는가? (스프레드-레벨 기울기 b)
    │     │
    │     ├─ 예, b ≈ 1  ──→ 로그 변환 (해석도 개선될 수 있음)
    │     ├─ 예, b ≈ 0.5 ──→ 제곱근 변환 (계수 자료)
    │     └─ 예, 그 외   ──→ x^(1-b) 또는 박스-콕스
    │
    └─ 아니오 (평균과 무관한 이분산)
          └─→ 웰치 분산분석  (+ 게임스-하웰 사후검정)

집단별 표준편차: 1.17, 3.67, 0.32  (최대/최소 = 11.5)
잔차 대 적합값 그림에서 뚜렷한 깔때기 모양
스프레드-레벨 기울기 b = 0.98 → 로그 변환 권장

→ 로그 척도에서 분산분석을 수행 (기하평균 비교로 해석)
   또는 원 척도에서 웰치 분산분석 (산술평균 비교)