Fisher 방법을 이용한 분산분석 수동 계산¶

import numpy as np
from scipy import stats

def manual_anova(groups):
    all_data = np.concatenate(list(groups.values()))
    grand_mean = all_data.mean()
    N = len(all_data)
    k = len(groups)

    # SST: 집단평균이 전체평균에서 얼마나 떨어져 있는가. n_i로 가중한다.
    # 큰 집단의 평균이 어긋나는 것이 더 무겁게 세어져야 하기 때문이다.
    SST = sum(len(g) * (g.mean() - grand_mean) ** 2
              for g in groups.values())
    # SSE: 각 관측값이 **자기 집단의** 평균에서 얼마나 떨어져 있는가.
    SSE = sum(np.sum((g - g.mean()) ** 2)
              for g in groups.values())

    MST = SST / (k - 1)
    MSE = SSE / (N - k)
    F = MST / MSE                    # 신호 대 잡음
    p_value = 1 - stats.f.cdf(F, k - 1, N - k)
    return SST, SSE, MST, MSE, F, p_value


# R의 PlantGrowth 자료 (대조군과 두 처리, 각 10개)
groups = {
    "ctrl": np.array([4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14]),
    "trt1": np.array([4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69]),
    "trt2": np.array([6.31, 5.12, 5.54, 5.50, 5.37, 5.29, 4.92, 6.15, 5.80, 5.26]),
}

SST, SSE, MST, MSE, F, p = manual_anova(groups)
print(f"SST = {SST:.4f}, SSE = {SSE:.4f}")
print(f"MST = {MST:.4f}, MSE = {MSE:.4f}")
print(f"F   = {F:.4f}, p = {p:.4f}")

SST = 3.7663, SSE = 10.4921
MST = 1.8832, MSE = 0.3886
F   = 4.8461, p = 0.0159

y = np.concatenate(list(groups.values()))
N, k = len(y), len(groups)
gbar = y.mean()
n_i = np.array([len(g) for g in groups.values()])
ybar_i = np.array([g.mean() for g in groups.values()])

# 정의대로 잰 것과 계산용 공식으로 잰 것을 나란히 둔다.
tot_def = ((y - gbar) ** 2).sum()
tot_mac = (y ** 2).sum() - N * gbar ** 2
sst_def = (n_i * (ybar_i - gbar) ** 2).sum()
sst_mac = (n_i * ybar_i ** 2).sum() - N * gbar ** 2
print(f"SS_total  정의 {tot_def!r}")
print(f"          계산용 {tot_mac!r}   차이 {abs(tot_def - tot_mac):.2e}")
print(f"SST       정의 {sst_def!r}")
print(f"          계산용 {sst_mac!r}   차이 {abs(sst_def - sst_mac):.2e}")

# R^2 과 F 의 일대일 관계
R2 = sst_def / tot_def
print(f"\nR^2 = SST / SS_total = {R2:.7f}")
print(f"F from R^2 = {(R2 / (k - 1)) / ((1 - R2) / (N - k)):.10f}")
print(f"F from SS  = {(sst_def / (k - 1)) / ((tot_def - sst_def) / (N - k)):.10f}")

# 자료를 통째로 1000 만큼 옮기면 제곱합은 그대로여야 한다.
z = y + 1000.0
tot_def2 = ((z - z.mean()) ** 2).sum()
tot_mac2 = (z ** 2).sum() - N * z.mean() ** 2
print(f"\ny + 1000 으로 옮긴 뒤")
print(f"  정의    {tot_def2:.9f}   (상대오차 {abs(tot_def2 / tot_def - 1):.1e})")
print(f"  계산용  {tot_mac2:.9f}   (상대오차 {abs(tot_mac2 / tot_def - 1):.1e})")

SS_total  정의 14.258429999999999
          계산용 14.258430000000203   차이 2.04e-13
SST       정의 3.766340000000002
          계산용 3.7663399999997864   차이 2.15e-13

R^2 = SST / SS_total = 0.2641483
F from R^2 = 4.8460878624
F from SS  = 4.8460878624

y + 1000 으로 옮긴 뒤
  정의    14.258430000   (상대오차 2.4e-14)
  계산용  14.258430004   (상대오차 2.9e-10)

# 손으로 구한 값과 맞는지 확인한다. 한 줄이면 되는 계산을 굳이 풀어 쓴 까닭은
# 제곱합이 어떻게 갈라지는지를 보이기 위해서다.
F_scipy, p_scipy = stats.f_oneway(*groups.values())
print(f"scipy: F = {F_scipy:.4f}, p = {p_scipy:.4f}")

scipy: F = 4.8461, p = 0.0159

m = 27
print(f"1 바로 아래의 배정밀도 수 = 1 - {1 - np.nextafter(1.0, 0.0):.6e}")
print(f"{'F':>9}{'1 - cdf':>15}{'sf':>15}{'closed':>15}{'relerr':>10}")
for f in [4.8461, 20, 100, 200, 500, 1000]:
    a = 1 - stats.f.cdf(f, 2, m)
    b = stats.f(2, m).sf(f)
    c = (1 + 2 * f / m) ** (-m / 2)   # d1=2 의 닫힌 꼴. 이것을 참값으로 삼는다.
    print(f"{f:>9.4f}{a:>15.6e}{b:>15.6e}{c:>15.6e}{abs(a - c) / c:>10.1e}")

1 바로 아래의 배정밀도 수 = 1 - 1.110223e-16
        F        1 - cdf             sf         closed    relerr
   4.8461   1.590982e-02   1.590982e-02   1.590982e-02   2.4e-15
  20.0000   4.692464e-06   4.692464e-06   4.692464e-06   1.7e-12
 100.0000   3.288481e-13   3.288984e-13   3.288984e-13   1.5e-04
 200.0000   1.110223e-16   6.495828e-17   6.495828e-17   7.1e-01
 500.0000   1.110223e-16   4.647398e-22   4.647398e-22   2.4e+05
1000.0000   1.110223e-16   4.796067e-26   4.796067e-26   2.3e+09

from itertools import combinations

def fisher_lsd(groups, MSE, alpha=0.05):
    """Fisher 의 최소유의차로 쌍별 비교를 한다.

    쌍마다 t-검정을 하되 표준오차를 그 두 집단이 아니라 전체 MSE 로 만든다.
    모든 집단의 정보를 쓰므로 자유도가 커지는 것이 이점이다.
    다만 다중비교를 보정하지 않으므로, 분산분석이 유의할 때만 쓴다.
    """
    names = list(groups.keys())
    N_total = sum(len(g) for g in groups.values())
    k = len(groups)
    df_within = N_total - k
    results = []
    for (n1, g1), (n2, g2) in combinations(groups.items(), 2):
        t_crit = stats.t.ppf(1 - alpha / 2, df_within)
        lsd_val = t_crit * np.sqrt(MSE * (1/len(groups[n1]) + 1/len(groups[n2])))
        diff = abs(groups[n1].mean() - groups[n2].mean())
        results.append({"pair": f"{n1} vs {n2}",
                        "diff": diff, "LSD": lsd_val,
                        "significant": diff > lsd_val})
    return results


for r in fisher_lsd(groups, MSE):
    print(f"{r['pair']:<14} diff = {r['diff']:.4f}  LSD = {r['LSD']:.4f}  {r['significant']}")

ctrl vs trt1   diff = 0.3710  LSD = 0.5720  False
ctrl vs trt2   diff = 0.4940  LSD = 0.5720  False
trt1 vs trt2   diff = 0.8650  LSD = 0.5720  True

t_crit = stats.t.ppf(0.975, 27)
q_crit = stats.studentized_range.ppf(0.95, 3, 27)
LSD = t_crit * np.sqrt(2 * MSE / 10)
HSD = q_crit * np.sqrt(MSE / 10)
print(f"t(0.975, 27) = {t_crit:.4f}   LSD = {LSD:.4f}")
print(f"q(0.05, 3, 27) = {q_crit:.4f}   HSD = {HSD:.4f}")
print(f"HSD / LSD = {HSD / LSD:.4f}   q/(sqrt(2) t) = {q_crit / (np.sqrt(2) * t_crit):.4f}")


def protected_lsd_fwer(mu, B=40000, n=10, alpha=0.05, seed=0):
    """평균이 mu 인 k 개 집단에서, 전역 F 가 기각한 뒤 LSD 를 돌렸을 때
    '참으로 같은 쌍' 가운데 하나라도 유의하다고 선언될 확률을 센다."""
    mu = np.asarray(mu, float)
    k = len(mu)
    N = n * k
    rng = np.random.default_rng(seed)
    Y = rng.normal(mu[None, :, None], 1.0, size=(B, k, n))
    gmean = Y.mean(axis=2)
    grand = Y.mean(axis=(1, 2))
    SSB = n * ((gmean - grand[:, None]) ** 2).sum(axis=1)
    SSE_ = ((Y - gmean[:, :, None]) ** 2).sum(axis=(1, 2))
    Fsim = (SSB / (k - 1)) / (SSE_ / (N - k))
    gate = Fsim > stats.f(k - 1, N - k).ppf(1 - alpha)
    thr = stats.t.ppf(1 - alpha / 2, N - k) * np.sqrt((SSE_ / (N - k)) * 2 / n)
    null_pairs = [(i, j) for i, j in combinations(range(k), 2) if mu[i] == mu[j]]
    any_false = np.zeros(B, bool)
    for i, j in null_pairs:
        any_false |= np.abs(gmean[:, i] - gmean[:, j]) > thr
    return np.mean(gate & any_false), gate.mean(), len(null_pairs)


print(f"\n{'design':>22}{'nullpairs':>11}{'global.rej':>12}{'FWER':>9}")
for mu, label in [((0, 0, 3), "k=3, mu=(0,0,3)"),
                  ((0, 0, 0, 0, 5), "k=5, mu=(0,0,0,0,5)"),
                  ((0, 0, 0, 0, 0, 0, 8), "k=7, mu=(0,...,0,8)")]:
    rate, gate, npair = protected_lsd_fwer(mu)
    print(f"{label:>22}{npair:>11}{gate:>12.4f}{rate:>9.4f}")
print(f"\nB=40000 에서 FWER 추정의 MC 표준오차 ~ {np.sqrt(0.05 * 0.95 / 40000):.4f}")

t(0.975, 27) = 2.0518   LSD = 0.5720
q(0.05, 3, 27) = 3.5064   HSD = 0.6912
HSD / LSD = 1.2084   q/(sqrt(2) t) = 1.2084

                design  nullpairs  global.rej     FWER
       k=3, mu=(0,0,3)          1      1.0000   0.0493
   k=5, mu=(0,0,0,0,5)          6      1.0000   0.1974
   k=7, mu=(0,...,0,8)         15      1.0000   0.3545

B=40000 에서 FWER 추정의 MC 표준오차 ~ 0.0011

import numpy as np
from scipy import stats
from itertools import combinations
from statsmodels.stats.multitest import multipletests
from statsmodels.stats.multicomp import pairwise_tukeyhsd

groups = {
    "ctrl": np.array([4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14]),
    "trt1": np.array([4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69]),
    "trt2": np.array([6.31, 5.12, 5.54, 5.50, 5.37, 5.29, 4.92, 6.15, 5.80, 5.26]),
}
names = list(groups)
values = np.concatenate([groups[n] for n in names])
N, k = len(values), len(names)
SSE = sum(((g - g.mean())**2).sum() for g in groups.values())
MSE = SSE / (N - k)
t_crit = stats.t.ppf(0.975, N - k)
print(f"MSE = {MSE:.4f},  df = {N - k},  t_crit = {t_crit:.4f}\n")

print("피셔 LSD")
pvals, labels = [], []
for a, b in combinations(range(k), 2):
    ga, gb = groups[names[a]], groups[names[b]]
    diff = ga.mean() - gb.mean()
    se = np.sqrt(MSE * (1 / len(ga) + 1 / len(gb)))
    t = diff / se
    p = 2 * stats.t.sf(abs(t), N - k)
    pvals.append(p)
    labels.append(f"{names[a]}-{names[b]}")
    print(f"  {names[a]}-{names[b]}:  차이 {diff:+.4f},  LSD {t_crit * se:.4f},"
          f"  t = {t:+.4f},  p = {p:.4f}"
          f"   {'유의' if abs(diff) > t_crit * se else ''}")

for method, name in [("bonferroni", "본페로니"), ("holm", "홀름  ")]:
    adj = multipletests(pvals, method=method)[1]
    print(f"  {name}: " + "   ".join(f"{labels[i]} {adj[i]:.4f}"
                                     for i in range(3)))

lab = np.repeat(names, [len(groups[n]) for n in names])
print("\n투키 HSD")
print(pairwise_tukeyhsd(values, lab, alpha=0.05))

MSE = 0.3886,  df = 27,  t_crit = 2.0518

피셔 LSD
  ctrl-trt1:  차이 +0.3710,  LSD 0.5720,  t = +1.3308,  p = 0.1944
  ctrl-trt2:  차이 -0.4940,  LSD 0.5720,  t = -1.7720,  p = 0.0877
  trt1-trt2:  차이 -0.8650,  LSD 0.5720,  t = -3.1028,  p = 0.0045   유의
  본페로니: ctrl-trt1 0.5832   ctrl-trt2 0.2630   trt1-trt2 0.0134
  홀름  : ctrl-trt1 0.1944   ctrl-trt2 0.1754   trt1-trt2 0.0134

투키 HSD
Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj   lower  upper  reject
---------------------------------------------------
  ctrl   trt1   -0.371 0.3909 -1.0622 0.3202  False
  ctrl   trt2    0.494  0.198 -0.1972 1.1852  False
  trt1   trt2    0.865  0.012  0.1738 1.5562   True
---------------------------------------------------

import numpy as np
from scipy import stats
from itertools import combinations
from statsmodels.stats.libqsturng import qsturng

rng = np.random.default_rng(1212)
M, n = 5_000, 10

print("① 완전 귀무: 모든 평균이 같다")
print(f"{'k':>3s} {'쌍':>4s} {'보호 LSD':>10s} {'무보호 LSD':>11s} "
      f"{'본페로니':>10s} {'투키':>8s}")
for k in [3, 4, 6, 10]:
    n_pair = k * (k - 1) // 2
    a = b = c = d = 0
    for _ in range(M):
        g = [rng.standard_normal(n) for _ in range(k)]
        N = k * n
        MSE = sum(((x - x.mean())**2).sum() for x in g) / (N - k)
        se = np.sqrt(MSE * 2 / n)
        diffs = [abs(g[i].mean() - g[j].mean())
                 for i, j in combinations(range(k), 2)]
        lsd_hit = any(dd > stats.t.ppf(0.975, N - k) * se for dd in diffs)
        b += lsd_hit
        a += (stats.f_oneway(*g).pvalue < 0.05) and lsd_hit
        c += any(dd > stats.t.ppf(1 - 0.025 / n_pair, N - k) * se
                 for dd in diffs)
        d += any(dd > qsturng(0.95, k, N - k) * np.sqrt(MSE / n)
                 for dd in diffs)
    print(f"{k:3d} {n_pair:4d} {a / M:10.4f} {b / M:11.4f} "
          f"{c / M:10.4f} {d / M:8.4f}")

① 완전 귀무: 모든 평균이 같다
  k    쌍     보호 LSD     무보호 LSD       본페로니       투키
  3    3     0.0520      0.1224     0.0442   0.0516
  4    6     0.0470      0.1954     0.0384   0.0476
  6   15     0.0468      0.3550     0.0362   0.0524
 10   45     0.0498      0.6098     0.0362   0.0536

rng = np.random.default_rng(3434)
print("② 부분 귀무: 한 집단만 멀리 떨어지고 나머지는 모두 같다")
print(f"{'k':>3s} {'F 기각률':>10s} {'보호 LSD':>10s} {'본페로니':>10s} "
      f"{'투키':>8s}")
for k in [3, 4, 6, 10]:
    mu = np.zeros(k)
    mu[0] = 4.0
    n_pair = k * (k - 1) // 2
    f_rej = a = c = d = 0
    for _ in range(M):
        g = [rng.normal(m, 1, n) for m in mu]
        N = k * n
        MSE = sum(((x - x.mean())**2).sum() for x in g) / (N - k)
        se = np.sqrt(MSE * 2 / n)
        p_f = stats.f_oneway(*g).pvalue
        f_rej += p_f < 0.05
        # 참으로 같은 집단들(1..k-1) 사이에서 거짓 발견이 있는가
        idx = list(range(1, k))
        eq = [abs(g[i].mean() - g[j].mean())
              for i, j in combinations(idx, 2)]
        if p_f < 0.05:
            a += any(dd > stats.t.ppf(0.975, N - k) * se for dd in eq)
        c += any(dd > stats.t.ppf(1 - 0.025 / n_pair, N - k) * se
                 for dd in eq)
        d += any(dd > qsturng(0.95, k, N - k) * np.sqrt(MSE / n)
                 for dd in eq)
    print(f"{k:3d} {f_rej / M:10.4f} {a / M:10.4f} {c / M:10.4f} "
          f"{d / M:8.4f}")

② 부분 귀무: 한 집단만 멀리 떨어지고 나머지는 모두 같다
  k      F 기각률     보호 LSD       본페로니       투키
  3     1.0000     0.0540     0.0192   0.0222
  4     1.0000     0.1204     0.0246   0.0290
  6     1.0000     0.2864     0.0272   0.0352
 10     1.0000     0.5598     0.0280   0.0406

import numpy as np
from scipy import stats

groups = {
    "ctrl": np.array([4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14]),
    "trt1": np.array([4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69]),
    "trt2": np.array([6.31, 5.12, 5.54, 5.50, 5.37, 5.29, 4.92, 6.15, 5.80, 5.26]),
}
names = list(groups)
y = np.concatenate([groups[n] for n in names])
N, k = len(y), len(names)

# 처리(treatment) 코딩: 절편 + 두 더미
X = np.zeros((N, k))
X[:, 0] = 1
for i in range(1, k):
    X[i * 10:(i + 1) * 10, i] = 1

beta, *_ = np.linalg.lstsq(X, y, rcond=None)
resid = y - X @ beta
SSE = resid @ resid
SS_total = ((y - y.mean())**2).sum()
SSR = SS_total - SSE
F = (SSR / (k - 1)) / (SSE / (N - k))
MSE = SSE / (N - k)
se = np.sqrt(np.diag(MSE * np.linalg.inv(X.T @ X)))

print(f"회귀계수 {np.round(beta, 4).tolist()}")
print(f"  절편   = ctrl 평균        {groups['ctrl'].mean():.4f}")
print(f"  β1     = trt1 − ctrl     "
      f"{groups['trt1'].mean() - groups['ctrl'].mean():+.4f}")
print(f"  β2     = trt2 − ctrl     "
      f"{groups['trt2'].mean() - groups['ctrl'].mean():+.4f}")
print(f"\nSSR = {SSR:.4f}  (= SST),  SSE = {SSE:.4f}")
print(f"F = {F:.4f},  p = {stats.f.sf(F, k - 1, N - k):.4f}")
print(f"scipy F = {stats.f_oneway(*groups.values()).statistic:.4f}")
print(f"\nR² = η² = {SSR / SS_total:.4f}")
omega2 = (SSR - (k - 1) * MSE) / (SS_total + MSE)
print(f"ω²      = {omega2:.4f}   (편향 보정한 효과크기)")

print(f"\n계수별 검정 (= ctrl 대비 비교)")
for i, label in enumerate(["절편(ctrl)", "trt1−ctrl ", "trt2−ctrl "]):
    t = beta[i] / se[i]
    print(f"  {label}  β = {beta[i]:+.4f},  SE = {se[i]:.4f},  "
          f"t = {t:+.4f},  p = {2 * stats.t.sf(abs(t), N - k):.4f}")

회귀계수 [5.032, -0.371, 0.494]
  절편   = ctrl 평균        5.0320
  β1     = trt1 − ctrl     -0.3710
  β2     = trt2 − ctrl     +0.4940

SSR = 3.7663  (= SST),  SSE = 10.4921
F = 4.8461,  p = 0.0159
scipy F = 4.8461

R² = η² = 0.2641
ω²      = 0.2041   (편향 보정한 효과크기)

계수별 검정 (= ctrl 대비 비교)
  절편(ctrl)  β = +5.0320,  SE = 0.1971,  t = +25.5265,  p = 0.0000
  trt1−ctrl   β = -0.3710,  SE = 0.2788,  t = -1.3308,  p = 0.1944
  trt2−ctrl   β = +0.4940,  SE = 0.2788,  t = +1.7720,  p = 0.0877

import numpy as np
from scipy import stats

rng = np.random.default_rng(5656)
M, sigma2 = 20_000, 4.0

print(f"{'μ':>18s} {'n':>4s} {'E[MSE]':>9s} {'σ²':>6s} "
      f"{'E[MST]':>10s} {'이론':>10s}")
for mu, n in [([0, 0, 0], 10), ([0, 1, 2], 10),
              ([0, 3, 6], 10), ([0, 1, 2], 30)]:
    mu = np.array(mu, float)
    k, N = len(mu), len(mu) * n
    mst, mse = [], []
    for _ in range(M):
        g = [rng.normal(m, np.sqrt(sigma2), n) for m in mu]
        allv = np.concatenate(g)
        gm = allv.mean()
        mst.append(sum(len(x) * (x.mean() - gm)**2 for x in g) / (k - 1))
        mse.append(sum(((x - x.mean())**2).sum() for x in g) / (N - k))
    theory = sigma2 + n * ((mu - mu.mean())**2).sum() / (k - 1)
    print(f"{str(mu.tolist()):>18s} {n:4d} {np.mean(mse):9.4f} "
          f"{sigma2:6.1f} {np.mean(mst):10.4f} {theory:10.4f}")

                 μ    n    E[MSE]     σ²     E[MST]         이론
   [0.0, 0.0, 0.0]   10    4.0092    4.0     3.9900     4.0000
   [0.0, 1.0, 2.0]   10    3.9855    4.0    13.9298    14.0000
   [0.0, 3.0, 6.0]   10    4.0061    4.0    94.0162    94.0000
   [0.0, 1.0, 2.0]   30    4.0061    4.0    34.0595    34.0000

import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)

sigma, k = 1.0, 3
mu = np.array([0, 0.5, 1.0])
print(f"\n총 N = 60 고정,  μ = {mu.tolist()},  σ = {sigma}")
print(f"{'배분':>16s} {'λ':>8s} {'검정력':>8s} {'1-3 쌍의 LSD 문턱':>18s}")
for ns in [(20, 20, 20), (30, 25, 5), (10, 20, 30), (40, 15, 5), (26, 26, 8)]:
    ns = np.array(ns)
    N = ns.sum()
    m_bar = (ns * mu).sum() / N
    lam = (ns * (mu - m_bar)**2).sum() / sigma**2
    crit = stats.f.ppf(0.95, k - 1, N - k)
    power = stats.ncf.sf(crit, k - 1, N - k, lam)
    lsd = stats.t.ppf(0.975, N - k) * sigma * np.sqrt(1 / ns[0] + 1 / ns[2])
    print(f"{str(ns.tolist()):>16s} {lam:8.4f} {power:8.4f} {lsd:18.4f}")

총 N = 60 고정,  μ = [0.0, 0.5, 1.0],  σ = 1.0
              배분        λ      검정력      1-3 쌍의 LSD 문턱
    [20, 20, 20]  10.0000   0.7933             0.6332
     [30, 25, 5]   6.1458   0.5710             0.9673
    [10, 20, 30]   8.3333   0.7120             0.7312
     [40, 15, 5]   6.1458   0.5710             0.9499
     [26, 26, 8]   7.1500   0.6407             0.8096

① 전체 평균 ȳ 와 집단 평균 ȳ_i 를 구한다
          ↓
② SST = Σ n_i (ȳ_i − ȳ)²        ← 집단 간
   SSE = Σ Σ (y_ij − ȳ_i)²       ← 집단 내
          ↓ 검산
   SS_total = Σ (y_ij − ȳ)²  =  SST + SSE
          ↓
③ MST = SST/(k−1),  MSE = SSE/(N−k)
          ↓
④ F = MST / MSE,  p = P(F_{k−1,N−k} > F)   ← 오른쪽 꼬리
          ↓
⑤ 효과크기 η² = SST/SS_total,  ω² (편향 보정)
          ↓
⑥ 유의하면 사후비교 (k≥4 이면 투키)