일원배치 분산분석 파이프라인¶

import pandas as pd
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm

url = ('https://raw.githubusercontent.com/vincentarelbundock/'
       'Rdatasets/1dcc2bf5f955cc1224a3e1307256e1fe86b68dae/csv/datasets/PlantGrowth.csv')
df = pd.read_csv(url, usecols=[1, 2])

# C()로 감싸지 않으면 group을 숫자처럼 취급해 회귀직선을 적합해 버린다.
# 문자열 열이면 statsmodels가 알아서 범주형으로 보지만, 수준이 0/1/2 같은
# 숫자로 코딩되어 있으면 조용히 틀린 모형이 된다. 습관적으로 감싸는 편이 안전하다.
model = ols('weight ~ C(group)', data=df).fit()
aov = anova_lm(model)
print(aov)

            df    sum_sq   mean_sq         F   PR(>F)
C(group)   2.0   3.76634  1.883170  4.846088  0.01591
Residual  27.0  10.49209  0.388596       NaN      NaN

means = df.groupby('group')['weight'].mean()
print(model.params)
print(f"\n절편             = ctrl 평균     = {means['ctrl']:.3f}")
print(f"C(group)[T.trt1] = trt1 - ctrl  = {means['trt1'] - means['ctrl']:.3f}")
print(f"C(group)[T.trt2] = trt2 - ctrl  = {means['trt2'] - means['ctrl']:.3f}")
print(f"\nR^2 = SSB/SST = {model.rsquared:.7f}"
      f"   (표에서 {3.76634 / (3.76634 + 10.49209):.7f})")
print(f"모형 자유도 {model.df_model:.0f},  잔차 자유도 {model.df_resid:.0f}")

# C() 를 빠뜨리면 어떻게 되는가. 수준을 0/1/2 로 코딩한 열을 만들어 넣어 본다.
df_num = df.assign(gcode=df['group'].map({'ctrl': 0, 'trt1': 1, 'trt2': 2}))
bad = ols('weight ~ gcode', data=df_num).fit()
print("\nC() 없이 숫자 코드를 그대로 넣으면")
print(anova_lm(bad))
print(f"기울기 = {bad.params['gcode']:.3f}  (집단당 '한 칸'씩 올라가는 직선)")

Intercept           5.032
C(group)[T.trt1]   -0.371
C(group)[T.trt2]    0.494
dtype: float64

절편             = ctrl 평균     = 5.032
C(group)[T.trt1] = trt1 - ctrl  = -0.371
C(group)[T.trt2] = trt2 - ctrl  = 0.494

R^2 = SSB/SST = 0.2641483   (표에서 0.2641483)
모형 자유도 2,  잔차 자유도 27

C() 없이 숫자 코드를 그대로 넣으면
            df    sum_sq   mean_sq        F    PR(>F)
gcode      1.0   1.22018  1.220180  2.62037  0.116711
Residual  28.0  13.03825  0.465652      NaN       NaN
기울기 = 0.247  (집단당 '한 칸'씩 올라가는 직선)

from statsmodels.stats.multicomp import pairwise_tukeyhsd

# 분산분석이 유의했으니 이제 어느 쌍이 다른지를 본다. reject 열이 True 인
# 쌍이 유의한 쌍이고, 신뢰구간이 0 을 품지 않는 쌍과 정확히 일치한다.
tukey = pairwise_tukeyhsd(endog=df['weight'], groups=df['group'], alpha=0.05)
print(tukey)

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

MSW, n, k, nu = 0.388596, 10, 3, 27
q_crit = stats.studentized_range.ppf(0.95, k, nu)
HSD = q_crit * np.sqrt(MSW / n)
print(f"q(0.05, 3, 27) = {q_crit:.4f}")
print(f"HSD = q * sqrt(MSW/n) = {HSD:.4f}   (세 쌍 공통)")

means = df.groupby('group')['weight'].mean()
print(f"\n{'pair':<12}{'diff':>8}{'lower':>9}{'upper':>8}{'q_obs':>8}{'p-adj':>8}{'reject':>8}")
for g1, g2 in [('ctrl', 'trt1'), ('ctrl', 'trt2'), ('trt1', 'trt2')]:
    d = means[g2] - means[g1]
    q_obs = abs(d) / np.sqrt(MSW / n)
    p_adj = stats.studentized_range.sf(q_obs, k, nu)
    print(f"{g1 + ' vs ' + g2:<12}{d:>8.3f}{d - HSD:>9.4f}{d + HSD:>8.4f}"
          f"{q_obs:>8.4f}{p_adj:>8.4f}{str(abs(d) > HSD):>8}")

q(0.05, 3, 27) = 3.5064
HSD = q * sqrt(MSW/n) = 0.6912   (세 쌍 공통)

pair            diff    lower   upper   q_obs   p-adj  reject
ctrl vs trt1  -0.371  -1.0622  0.3202  1.8820  0.3909   False
ctrl vs trt2   0.494  -0.1972  1.1852  2.5060  0.1980   False
trt1 vs trt2   0.865   0.1738  1.5562  4.3880  0.0120    True

from itertools import combinations
from scipy.stats import ttest_ind
from statsmodels.stats.multitest import multipletests

groups = df['group'].unique()
p_raw, labels = [], []
for g1, g2 in combinations(groups, 2):
    x = df.loc[df['group'] == g1, 'weight'].values
    y = df.loc[df['group'] == g2, 'weight'].values
    stat, p = ttest_ind(x, y, equal_var=False)
    p_raw.append(p)
    labels.append(f"{g1} vs {g2}")

# Bonferroni는 문턱을 낮추는 대신 p-값에 m을 곱해 돌려준다.
# 그래서 보정 후에도 비교 대상은 여전히 alpha다.
_, p_bonf, _, _ = multipletests(p_raw, alpha=0.05, method='bonferroni')
for lbl, p, pb in zip(labels, p_raw, p_bonf):
    print(f"{lbl:<12}  p = {p:.4f}   p_bonf = {pb:.4f}")

ctrl vs trt1  p = 0.2504   p_bonf = 0.7511
ctrl vs trt2  p = 0.0479   p_bonf = 0.1437
trt1 vs trt2  p = 0.0093   p_bonf = 0.0279

MSW, n, k, nu, m = 0.388596, 10, 3, 27, 3

# 본페로니와 Tukey 의 문턱을 같은 자에 올려 견준다 (둘 다 합동분산 기준).
t_b = stats.t.ppf(1 - 0.05 / (2 * m), nu)
thr_bonf = t_b * np.sqrt(2 * MSW / n)
q_crit = stats.studentized_range.ppf(0.95, k, nu)
thr_tukey = q_crit * np.sqrt(MSW / n)
print(f"본페로니 문턱 = {t_b:.4f} * {np.sqrt(2 * MSW / n):.4f} = {thr_bonf:.4f}")
print(f"Tukey  문턱   = {q_crit:.4f} * {np.sqrt(MSW / n):.4f} = {thr_tukey:.4f}")
print(f"비 = {thr_bonf / thr_tukey:.5f}   (= sqrt(2) t / q = {np.sqrt(2) * t_b / q_crit:.5f})")

# 합동분산 t 로 본페로니를 다시 하면 Welch 와 얼마나 다른가.
print(f"\n{'pair':<14}{'Welch p':>10}{'Welch bonf':>12}"
      f"{'pooled p':>10}{'pooled bonf':>13}{'Tukey adj':>11}")
means = df.groupby('group')['weight'].mean()
for g1, g2 in combinations(['ctrl', 'trt1', 'trt2'], 2):
    x = df.loc[df['group'] == g1, 'weight'].values
    y = df.loc[df['group'] == g2, 'weight'].values
    p_w = stats.ttest_ind(x, y, equal_var=False).pvalue
    d = means[g2] - means[g1]
    t_pool = d / np.sqrt(2 * MSW / n)
    p_pool = 2 * stats.t.sf(abs(t_pool), nu)
    p_tuk = stats.studentized_range.sf(abs(d) / np.sqrt(MSW / n), k, nu)
    print(f"{g1 + ' vs ' + g2:<14}{p_w:>10.4f}{min(m * p_w, 1):>12.4f}"
          f"{p_pool:>10.4f}{min(m * p_pool, 1):>13.4f}{p_tuk:>11.4f}")

본페로니 문턱 = 2.5525 * 0.2788 = 0.7116
Tukey  문턱   = 3.5064 * 0.1971 = 0.6912
비 = 1.02946   (= sqrt(2) t / q = 1.02946)

pair             Welch p  Welch bonf  pooled p  pooled bonf  Tukey adj
ctrl vs trt1      0.2504      0.7511    0.1944       0.5832     0.3909
ctrl vs trt2      0.0479      0.1437    0.0877       0.2630     0.1980
trt1 vs trt2      0.0093      0.0279    0.0045       0.0134     0.0120

import matplotlib.pyplot as plt

# 상자그림으로 마무리한다. 검정 결과와 그림이 같은 이야기를 하는지 확인하는
# 것이 마지막 단계다. 순서를 못박아 두어야 그림이 자료 순서에 휘둘리지 않는다.
order = ['ctrl', 'trt1', 'trt2']
data = [df.loc[df['group'] == g, 'weight'].values for g in order]
plt.boxplot(data, labels=order)
plt.xlabel('Group')
plt.ylabel('Weight')
plt.title('PlantGrowth weights by group')
plt.tight_layout()
plt.show()

HSD = 0.6912
box = {g: np.percentile(df.loc[df['group'] == g, 'weight'], [25, 50, 75])
       for g in order}
means = df.groupby('group')['weight'].mean()
print(f"{'pair':<14}{'box overlap':>12}{'|diff|':>9}{'HSD':>8}{'Tukey':>8}")
for g1, g2 in [('ctrl', 'trt1'), ('ctrl', 'trt2'), ('trt1', 'trt2')]:
    lo = max(box[g1][0], box[g2][0])     # 두 상자 아랫변 중 높은 쪽
    hi = min(box[g1][2], box[g2][2])     # 두 상자 윗변 중 낮은 쪽
    overlap = max(hi - lo, 0.0)
    d = abs(means[g2] - means[g1])
    print(f"{g1 + ' vs ' + g2:<14}{overlap:>12.3f}{d:>9.3f}{HSD:>8.4f}"
          f"{str(d > HSD):>8}")

pair           box overlap   |diff|     HSD   Tukey
ctrl vs trt1         0.320    0.371  0.6912   False
ctrl vs trt2         0.025    0.494  0.6912   False
trt1 vs trt2         0.000    0.865  0.6912    True

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)

print("집단별 요약")
for nm in names:
    g = groups[nm]
    print(f"  {nm}: n={len(g)}  평균={g.mean():.4f}  "
          f"표준편차={g.std(ddof=1):.4f}  분산={g.var(ddof=1):.4f}")
v = [groups[nm].var(ddof=1) for nm in names]
print(f"  분산비 최대/최소 = {max(v) / min(v):.4f}\n")

print("본페로니 보정 Welch t")
pvals, labels = [], []
for a, b in combinations(range(3), 2):
    r = stats.ttest_ind(groups[names[a]], groups[names[b]], equal_var=False)
    pvals.append(r.pvalue)
    labels.append(f"{names[a]}-{names[b]}")
    print(f"  {names[a]}-{names[b]}:  t = {r.statistic:+.4f},  "
          f"df = {r.df:.2f},  p = {r.pvalue:.4f}")
adj = multipletests(pvals, method="bonferroni")[1]
print("  본페로니 조정 p: "
      + "   ".join(f"{labels[i]} {adj[i]:.4f}" for i in range(3)))

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

집단별 요약
  ctrl: n=10  평균=5.0320  표준편차=0.5831  분산=0.3400
  trt1: n=10  평균=4.6610  표준편차=0.7937  분산=0.6299
  trt2: n=10  평균=5.5260  표준편차=0.4426  분산=0.1959
  분산비 최대/최소 = 3.2160

본페로니 보정 Welch t
  ctrl-trt1:  t = +1.1913,  df = 16.52,  p = 0.2504
  ctrl-trt2:  t = -2.1340,  df = 16.79,  p = 0.0479
  trt1-trt2:  t = -3.0101,  df = 14.10,  p = 0.0093
  본페로니 조정 p: ctrl-trt1 0.7511   ctrl-trt2 0.1437   trt1-trt2 0.0279

투키 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 psturng, qsturng

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)
k = len(names)

print("게임스·하월")
for a, b in combinations(range(k), 2):
    x, y = groups[names[a]], groups[names[b]]
    n1, n2 = len(x), len(y)
    v1, v2 = x.var(ddof=1), y.var(ddof=1)
    se = np.sqrt((v1 / n1 + v2 / n2) / 2)
    df = ((v1 / n1 + v2 / n2)**2
          / ((v1 / n1)**2 / (n1 - 1) + (v2 / n2)**2 / (n2 - 1)))
    q = abs(x.mean() - y.mean()) / se
    p = float(np.atleast_1d(psturng(q, k, df))[0])
    half = qsturng(0.95, k, df) * se
    diff = x.mean() - y.mean()
    print(f"  {names[a]}-{names[b]}: 차이 {diff:+.4f}  q = {q:.4f}  "
          f"df = {df:.2f}  p = {p:.4f}")
    print(f"              95% CI ({diff - half:+.4f}, {diff + half:+.4f})")

게임스·하월
  ctrl-trt1: 차이 +0.3710  q = 1.6847  df = 16.52  p = 0.4761
              95% CI (-0.4299, +1.1719)
  ctrl-trt2: 차이 -0.4940  q = 3.0180  df = 16.79  p = 0.1128
              95% CI (-1.0884, +0.1004)
  trt1-trt2: 차이 -0.8650  q = 4.2569  df = 14.10  p = 0.0236
              95% CI (-1.6162, -0.1138)

rng = np.random.default_rng(1357)

def tukey_any(g, alpha=0.05):
    k = len(g)
    n = np.array([len(x) for x in g])
    N = n.sum()
    MSE = sum(((x - x.mean())**2).sum() for x in g) / (N - k)
    for i, j in combinations(range(k), 2):
        se = np.sqrt(MSE / 2 * (1 / n[i] + 1 / n[j]))
        if abs(g[i].mean() - g[j].mean()) / se > qsturng(1 - alpha, k, N - k):
            return True
    return False

def games_howell_any(g, alpha=0.05):
    k = len(g)
    for i, j in combinations(range(k), 2):
        x, y = g[i], g[j]
        n1, n2 = len(x), len(y)
        v1, v2 = x.var(ddof=1), y.var(ddof=1)
        se = np.sqrt((v1 / n1 + v2 / n2) / 2)
        df = ((v1 / n1 + v2 / n2)**2
              / ((v1 / n1)**2 / (n1 - 1) + (v2 / n2)**2 / (n2 - 1)))
        if abs(x.mean() - y.mean()) / se > qsturng(1 - alpha, k, df):
            return True
    return False

def bonf_welch_any(g, alpha=0.05):
    k = len(g)
    m = k * (k - 1) // 2
    return any(stats.ttest_ind(g[i], g[j], equal_var=False).pvalue < alpha / m
               for i, j in combinations(range(k), 2))

M = 5_000
print(f"{'상황':>24s} {'투키':>8s} {'게임스·하월':>12s} {'본페로니 Welch':>15s}")
for label, sig, ns in [("등분산 등n", (1, 1, 1), (10, 10, 10)),
                       ("이분산 등n", (1, 1, 3), (10, 10, 10)),
                       ("이분산 불균형(큰σ↔작은n)", (1, 1, 3), (15, 15, 5)),
                       ("이분산 불균형(큰σ↔큰n)", (1, 1, 3), (5, 15, 15))]:
    a = b = c = 0
    for _ in range(M):
        g = [rng.normal(0, s, n) for s, n in zip(sig, ns)]
        a += tukey_any(g)
        b += games_howell_any(g)
        c += bonf_welch_any(g)
    print(f"{label:>24s} {a / M:8.4f} {b / M:12.4f} {c / M:15.4f}")

                      상황       투키       게임스·하월      본페로니 Welch
                  등분산 등n   0.0464       0.0470          0.0400
                  이분산 등n   0.0702       0.0492          0.0414
         이분산 불균형(큰σ↔작은n)   0.2538       0.0506          0.0392
          이분산 불균형(큰σ↔큰n)   0.0326       0.0506          0.0436

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

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

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)
k, n = len(names), 10
N = k * n
values = np.concatenate([groups[x] for x in names])
grand = values.mean()
m = np.array([groups[x].mean() for x in names])
SST = n * ((m - grand)**2).sum()
SSE = sum(((x - x.mean())**2).sum() for x in groups.values())
SS_total = ((values - grand)**2).sum()
MSE = SSE / (N - k)
F = (SST / (k - 1)) / MSE

eta2 = SST / SS_total
omega2 = (SST - (k - 1) * MSE) / (SS_total + MSE)
print(f"F = {F:.4f},  p = {stats.f.sf(F, k - 1, N - k):.4f}")
print(f"η² = {eta2:.4f}   ω² = {omega2:.4f}   "
      f"Cohen f = {np.sqrt(eta2 / (1 - eta2)):.4f}")

def ncp_ci(F_obs, df1, df2, alpha=0.05, big=1e5):
    lo = (brentq(lambda l: stats.ncf.sf(F_obs, df1, df2, l) - alpha / 2, 0, big)
          if stats.ncf.sf(F_obs, df1, df2, 0) < alpha / 2 else 0.0)
    hi = (brentq(lambda l: stats.ncf.cdf(F_obs, df1, df2, l) - alpha / 2, 0, big)
          if stats.ncf.cdf(F_obs, df1, df2, 0) > alpha / 2 else 0.0)
    return lo, hi

lo, hi = ncp_ci(F, k - 1, N - k)
print(f"λ 의 95% CI ({lo:.4f}, {hi:.4f})")
print(f"η² 의 95% CI ({lo / (lo + N):.4f}, {hi / (hi + N):.4f})")

q = qsturng(0.95, k, N - k)
print(f"\n투키 신뢰구간 (q = {q:.4f})")
for a, b in combinations(range(k), 2):
    diff = m[a] - m[b]
    half = q * np.sqrt(MSE / n)
    print(f"  {names[a]}-{names[b]}: {diff:+.4f}  "
          f"95% CI ({diff - half:+.4f}, {diff + half:+.4f})")

F = 4.8461,  p = 0.0159
η² = 0.2641   ω² = 0.2041   Cohen f = 0.5991
λ 의 95% CI (0.3001, 25.9593)
η² 의 95% CI (0.0099, 0.4639)

투키 신뢰구간 (q = 3.5058)
  ctrl-trt1: +0.3710  95% CI (-0.3201, +1.0621)
  ctrl-trt2: -0.4940  95% CI (-1.1851, +0.1971)
  trt1-trt2: -0.8650  95% CI (-1.5561, -0.1739)

import numpy as np
from scipy import stats

def anova_from_summary(means, sds, ns):
    """집단 평균·표준편차·크기만으로 분산분석표를 만든다."""
    means = np.asarray(means, float)
    sds = np.asarray(sds, float)
    ns = np.asarray(ns, float)
    k = len(means)
    N = ns.sum()
    grand = (ns * means).sum() / N
    SST = (ns * (means - grand)**2).sum()
    SSE = ((ns - 1) * sds**2).sum()
    MST, MSE = SST / (k - 1), SSE / (N - k)
    F = MST / MSE
    return {"SST": SST, "SSE": SSE, "MST": MST, "MSE": MSE, "F": F,
            "df1": k - 1, "df2": int(N - k),
            "p": stats.f.sf(F, k - 1, N - k),
            "eta2": SST / (SST + SSE)}

out = anova_from_summary([5.0320, 4.6610, 5.5260],
                         [0.5831, 0.7937, 0.4426],
                         [10, 10, 10])
for key in ["SST", "SSE", "MST", "MSE", "F", "df1", "df2", "p", "eta2"]:
    val = out[key]
    print(f"  {key:5s} = {val:.4f}" if isinstance(val, float)
          else f"  {key:5s} = {val}")

# 균형 설계의 축약 공식으로도 같은 F 가 나오는지 확인
means = np.array([5.0320, 4.6610, 5.5260])
sds = np.array([0.5831, 0.7937, 0.4426])
n = 10
F_short = n * means.var(ddof=1) / (sds**2).mean()
print(f"\n축약 공식  F = n·s²(평균) / 평균(s²) = {F_short:.4f}")

  SST   = 3.7663
  SSE   = 10.4927
  MST   = 1.8832
  MSE   = 0.3886
  F     = 4.8458
  df1   = 2
  df2   = 27
  p     = 0.0159
  eta2  = 0.2641

축약 공식  F = n·s²(평균) / 평균(s²) = 4.8458

def welch_from_summary(means, sds, ns):
    means = np.asarray(means, float)
    v = np.asarray(sds, float)**2
    n = np.asarray(ns, float)
    k = len(means)
    w = n / v
    W = w.sum()
    m_tilde = (w * means).sum() / W
    tmp = np.sum((1 - w / W)**2 / (n - 1))
    F = ((w * (means - m_tilde)**2).sum() / (k - 1)) \
        / (1 + 2 * (k - 2) / (k * k - 1) * tmp)
    df2 = (k * k - 1) / (3 * tmp)
    return F, df2, stats.f.sf(F, k - 1, df2)

F_w, df2_w, p_w = welch_from_summary([5.0320, 4.6610, 5.5260],
                                      [0.5831, 0.7937, 0.4426],
                                      [10, 10, 10])
print(f"Welch:  F = {F_w:.4f},  df2 = {df2_w:.2f},  p = {p_w:.4f}")
print(f"고전 F: F = {out['F']:.4f},  df2 = {out['df2']},  p = {out['p']:.4f}")

Welch:  F = 5.1805,  df2 = 17.13,  p = 0.0174
고전 F: F = 4.8458,  df2 = 27,  p = 0.0159

① 자료 탐색
    집단별 n·평균·표준편차 표, 상자그림, 정규분위수그림
          ↓
② 가정 점검 (우선순위 순)
    독립성 → 등분산(분산비) → 정규성
          ↓
③ 검정 선택 (사전에 정한다)
    등분산·등n → 고전 F   /   그 밖 → Welch
          ↓
④ 옴니버스 검정
          ↓
⑤ 효과크기 + 신뢰구간
    η², ω², Cohen f, 그리고 η² 의 구간
          ↓
⑥ 사후비교 (유의할 때만)
    등분산 → 투키   /   이분산 → 게임스·하월
    대조군 대비 → 더넷
          ↓
⑦ 시각화
    집단별 상자그림 + 사후비교 신뢰구간 그림