분산분석 등분산성 사전검정¶

import numpy as np
from scipy import stats

def welch_anova(groups):
    """Welch 의 일원배치 분산분석. (F, df1, df2, p) 를 돌려준다.

    등분산을 가정하지 않는다. 집단마다 분산의 역수로 가중해 평균을 내고,
    분모 자유도를 실수로 조정한다.
    """
    k = len(groups)
    n = np.array([len(g) for g in groups])
    m = np.array([g.mean() for g in groups])
    v = np.array([g.var(ddof=1) for g in groups])
    w = n / v
    m_w = np.sum(w * m) / np.sum(w)
    A = np.sum(w * (m - m_w) ** 2) / (k - 1)
    lam = np.sum((1 - w / np.sum(w)) ** 2 / (n - 1)) / (k ** 2 - 1)
    F = A / (1 + 2 * (k - 2) * lam)
    df2 = 1 / (3 * lam)
    return F, k - 1, df2, stats.f.sf(F, k - 1, df2)

# 네 집단 모두 표준편차가 5 로 같다. 곧 등분산이 참인 자료다.
rng = np.random.default_rng(42)
g1 = rng.normal(50, 5, size=15)
g2 = rng.normal(55, 5, size=15)
g3 = rng.normal(52, 5, size=15)
g4 = rng.normal(48, 5, size=15)
groups = [g1, g2, g3, g4]

# 1단계: 등분산 사전검정. 그런데 이 절차 자체가 문제를 안고 있다.
# 자료를 보고 다음 검정을 고르면, 최종 결과의 제1종 오류율이 명목수준을
# 넘어선다. 사전검정 없이 처음부터 Welch 를 쓰라는 권고가 나오는 까닭이다.
bf_stat, bf_p = stats.levene(*groups, center='median')
print(f"Brown-Forsythe: W = {bf_stat:.4f}, p = {bf_p:.4f}")

# 2단계: 두 분산분석을 모두 돌려 견준다. 등분산이 참인 자료에서는
# Welch 가 잃는 것이 거의 없다. 그래서 그냥 Welch 를 쓰면 된다.
f_std, p_std = stats.f_oneway(*groups)
f_w, df1, df2, p_w = welch_anova(groups)
print(f"Standard ANOVA: F = {f_std:.4f}, p = {p_std:.6f}")
print(f"Welch ANOVA:    F = {f_w:.4f}, df = ({df1}, {df2:.2f}), "
      f"p = {p_w:.6f}")

Brown-Forsythe: W = 0.4728, p = 0.7025
Standard ANOVA: F = 8.1365, p = 0.000138
Welch ANOVA:    F = 9.7514, df = (3, 30.81), p = 0.000112

import numpy as np
from scipy import stats


def anova_p(Y):
    """균형자료 (R, k, n) 에 대한 표준 일원분산분석 p-값을 한꺼번에."""
    R, k, n = Y.shape
    m = Y.mean(axis=2)
    gm = Y.reshape(R, -1).mean(axis=1)
    msb = n * ((m - gm[:, None]) ** 2).sum(axis=1) / (k - 1)
    mse = Y.var(axis=2, ddof=1).mean(axis=1)   # 균형자료에서는 합동분산 = 평균
    return stats.f.sf(msb / mse, k - 1, k * (n - 1)), msb, mse


def bf_p(Y):
    """Brown-Forsythe = 중앙값 중심 Levene. |y - median| 에 분산분석을 돌린다.
    scipy.stats.levene 의 center 기본값이 'median' 이므로 그것과 같은 검정이다."""
    return anova_p(np.abs(Y - np.median(Y, axis=2, keepdims=True)))[0]


def welch_p(Y):
    R, k, n = Y.shape
    m, v = Y.mean(axis=2), Y.var(axis=2, ddof=1)
    w = n / v
    sw = w.sum(axis=1)
    m_w = (w * m).sum(axis=1) / sw
    A = (w * (m - m_w[:, None]) ** 2).sum(axis=1) / (k - 1)
    lam = ((1 - w / sw[:, None]) ** 2 / (n - 1)).sum(axis=1) / (k ** 2 - 1)
    return stats.f.sf(A / (1 + 2 * (k - 2) * lam), k - 1, 1 / (3 * lam))


rng = np.random.default_rng(1)
R, n, alpha = 200_000, 15, 0.05

for sds in [(5, 5, 5, 5), (5, 5, 5, 10), (2, 4, 6, 10)]:
    sd = np.array(sds, float)
    k = len(sd)
    # 평균이 모두 0 이므로 H0 이 참이다. 기각하면 모두 1종오류다.
    Y = rng.normal(0, 1, size=(R, k, n)) * sd[None, :, None]

    pre = bf_p(Y)
    p_std, msb, mse = anova_p(Y)
    p_w = welch_p(Y)

    A = pre <= alpha                      # 사전검정이 기각한 표본
    pi = A.mean()
    a_std, a_w = (p_std < alpha).mean(), (p_w < alpha).mean()
    a_two = np.where(A, p_w < alpha, p_std < alpha).mean()
    pred = (1 - pi) * a_std + pi * a_w     # (1) 의 볼록결합

    v = Y.var(axis=2, ddof=1)
    ratio = v.max(axis=1) / v.min(axis=1)

    print(f"sd = {sds}")
    print(f"   사전검정 기각률 pi      = {pi:.4f}")
    print(f"   항상 표준 분산분석      = {a_std:.4f}")
    print(f"   항상 Welch              = {a_w:.4f}")
    print(f"   두 단계 (실제)          = {a_two:.4f}")
    print(f"   두 단계 (독립가정 예측) = {pred:.4f}      차 = {a_two - pred:+.4f}")
    print(f"   P(표준이 기각 | 사전검정 통과) = {(p_std < alpha)[~A].mean():.4f}"
          f"   (무조건 {a_std:.4f})")
    print(f"   E[MSB]  전체 {msb.mean():7.2f}  ->  통과 {msb[~A].mean():7.2f}"
          f"   ({msb[~A].mean() / msb.mean():.3f} 배)")
    print(f"   E[MSE]  전체 {mse.mean():7.2f}  ->  통과 {mse[~A].mean():7.2f}"
          f"   ({mse[~A].mean() / mse.mean():.3f} 배)")
    print(f"   E[s^2_max/s^2_min] 전체 {ratio.mean():6.2f}  ->  통과 {ratio[~A].mean():6.2f}")

sd = (5, 5, 5, 5)
   사전검정 기각률 pi      = 0.0276
   항상 표준 분산분석      = 0.0488
   항상 Welch              = 0.0493
   두 단계 (실제)          = 0.0498
   두 단계 (독립가정 예측) = 0.0488      차 = +0.0010
   P(표준이 기각 | 사전검정 통과) = 0.0487   (무조건 0.0488)
   E[MSB]  전체   24.90  ->  통과   24.89   (1.000 배)
   E[MSE]  전체   25.00  ->  통과   25.00   (1.000 배)
   E[s^2_max/s^2_min] 전체   2.39  ->  통과   2.32
sd = (5, 5, 5, 10)
   사전검정 기각률 pi      = 0.6108
   항상 표준 분산분석      = 0.0663
   항상 Welch              = 0.0492
   두 단계 (실제)          = 0.0617
   두 단계 (독립가정 예측) = 0.0559      차 = +0.0059
   P(표준이 기각 | 사전검정 통과) = 0.0861   (무조건 0.0663)
   E[MSB]  전체   43.55  ->  통과   43.58   (1.001 배)
   E[MSE]  전체   43.72  ->  통과   38.84   (0.888 배)
   E[s^2_max/s^2_min] 전체   6.38  ->  통과   4.01
sd = (2, 4, 6, 10)
   사전검정 기각률 pi      = 0.9798
   항상 표준 분산분석      = 0.0777
   항상 Welch              = 0.0512
   두 단계 (실제)          = 0.0524
   두 단계 (독립가정 예측) = 0.0518      차 = +0.0007
   P(표준이 기각 | 사전검정 통과) = 0.1477   (무조건 0.0777)
   E[MSB]  전체   39.04  ->  통과   39.29   (1.006 배)
   E[MSE]  전체   39.01  ->  통과   28.47   (0.730 배)
   E[s^2_max/s^2_min] 전체  29.29  ->  통과  11.93

import numpy as np
from scipy import stats

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

def sim(sds, R=5000, n=15, seed=1):
    rng = np.random.default_rng(seed)
    a = b = c = 0
    for _ in range(R):
        gs = [rng.normal(0, s, n) for s in sds]
        p_bf = stats.levene(*gs, center='median')[1]
        p_std = stats.f_oneway(*gs)[1]
        p_w = welch_p(gs)
        a += p_std < 0.05
        b += p_w < 0.05
        c += (p_std < 0.05) if p_bf > 0.05 else (p_w < 0.05)
    return a / R, b / R, c / R

for sds in [(5, 5, 5, 5), (5, 5, 5, 10), (2, 4, 6, 10)]:
    std, welch, twostage = sim(sds)
    print(f"sd = {str(sds):>16}: standard {std:.4f}, "
          f"Welch {welch:.4f}, two-stage {twostage:.4f}")

sd =     (5, 5, 5, 5): standard 0.0452, Welch 0.0452, two-stage 0.0458
sd =    (5, 5, 5, 10): standard 0.0634, Welch 0.0440, two-stage 0.0568
sd =    (2, 4, 6, 10): standard 0.0738, Welch 0.0464, two-stage 0.0468

import numpy as np
from scipy import stats

rng = np.random.default_rng(7)
R, n = 5000, 15

for mus in [(0, 0, 0, 3), (0, 0, 0, 5)]:
    a = b = 0
    for _ in range(R):
        gs = [rng.normal(m, 5, n) for m in mus]
        a += stats.f_oneway(*gs)[1] < 0.05
        b += welch_p(gs) < 0.05          # welch_p from Exercise 1
    print(f"mu = {str(mus):>14}: standard {a/R:.4f}, Welch {b/R:.4f}")

mu =   (0, 0, 0, 3): standard 0.3368, Welch 0.3190
mu =   (0, 0, 0, 5): standard 0.7804, Welch 0.7536