분산분석의 실무 응용¶

import pandas as pd
import seaborn as sns
import statsmodels.api as sm
from statsmodels.formula.api import ols

# 붓꽃 자료에서 두 품종만 남긴다. 집단이 둘이면 분산분석과 이표본 t-검정이
# 같은 결론을 주며, F = t^2 이라는 관계도 확인할 수 있다.
data = sns.load_dataset("iris")
data = data[data["species"] != "setosa"]

model = ols('sepal_length ~ species', data=data).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)

           sum_sq    df          F        PR(>F)
species   10.6276   1.0  31.687502  1.724856e-07
Residual  32.8680  98.0        NaN           NaN

import numpy as np
import seaborn as sns
from scipy import stats
from statsmodels.formula.api import ols

data = sns.load_dataset("iris")
data = data[data["species"] != "setosa"]
g1 = data[data.species == 'versicolor'].sepal_length.values
g2 = data[data.species == 'virginica'].sepal_length.values
n1 = n2 = 50

print(f"versicolor: 평균 {g1.mean():.4f}, 표준편차 {g1.std(ddof=1):.4f}")
print(f"virginica : 평균 {g2.mean():.4f}, 표준편차 {g2.std(ddof=1):.4f}")
delta = g1.mean() - g2.mean()
SSB = n1 * n2 / (n1 + n2) * delta ** 2
print(f"\n차이 = {delta:.4f}")
print(f"공식 SSB = n1*n2/(n1+n2) * 차이^2 = {SSB:.6f}")

model = ols('sepal_length ~ species', data=data).fit()
MSE = model.mse_resid
t = stats.ttest_ind(g1, g2).statistic
print(f"MSE = s_p^2 = {MSE:.6f},  s_p = {np.sqrt(MSE):.6f}")
print(f"공식 F = SSB/MSE = {SSB / MSE:.6f}")
print(f"statsmodels F = {model.fvalue:.6f}")
print(f"이표본 t = {t:.6f},  t^2 = {t ** 2:.6f}")

print(f"\nCohen d = 차이 / s_p = {abs(delta) / np.sqrt(MSE):.4f}")
print(f"참고: 집단당 n = 25 였다면 F = {25 / 50 * SSB / MSE:.2f}, "
      f"p = {stats.f.sf(25 / 50 * SSB / MSE, 1, 48):.2e}")

versicolor: 평균 5.9360, 표준편차 0.5162
virginica : 평균 6.5880, 표준편차 0.6359

차이 = -0.6520
공식 SSB = n1*n2/(n1+n2) * 차이^2 = 10.627600
MSE = s_p^2 = 0.335388,  s_p = 0.579127
공식 F = SSB/MSE = 31.687502
statsmodels F = 31.687502
이표본 t = -5.629165,  t^2 = 31.687502

Cohen d = 차이 / s_p = 1.1258
참고: 집단당 n = 25 였다면 F = 15.84, p = 2.32e-04

import matplotlib.pyplot as plt
from scipy.stats import shapiro

# 정규성은 자료가 아니라 잔차에 요구된다.
sm.qqplot(model.resid, line='s')
plt.title("Q-Q Plot of Residuals")
plt.show()

# 표본이 크면 사소한 이탈에도 유의하게 나오므로 그림과 함께 읽는다.
stat, p_value = shapiro(model.resid)
print(f"Shapiro-Wilk Test: W = {stat:.4f}, p-value = {p_value:.4f}")

Shapiro-Wilk Test: W = 0.9831, p-value = 0.2285

import numpy as np
from scipy import stats

rng = np.random.default_rng(2024)
B = 20_000
dists = {
    "정규": lambda n: rng.normal(0, 1, n),
    "t(3)": lambda n: rng.standard_t(3, n),
    "지수": lambda n: rng.exponential(1, n),
    "로그정규": lambda n: np.exp(rng.normal(0, 1, n)),
}

print("집단당 n = 50, 두 집단, 평균이 실제로 같을 때 (명목 0.05)")
print(f"{'분포':>10}{'F 검정 오류율':>14}{'Shapiro 기각률':>16}")
for lab, f in dists.items():
    bad_f = bad_s = 0
    for _ in range(B):
        a, b = f(50), f(50)
        bad_f += stats.f_oneway(a, b).pvalue < 0.05
        resid = np.concatenate([a - a.mean(), b - b.mean()])
        bad_s += stats.shapiro(resid).pvalue < 0.05
    print(f"{lab:>10}{bad_f / B:>14.4f}{bad_s / B:>16.4f}")

집단당 n = 50, 두 집단, 평균이 실제로 같을 때 (명목 0.05)
        분포      F 검정 오류율     Shapiro 기각률
        정규        0.0510          0.0478
      t(3)        0.0493          0.8699
        지수        0.0501          1.0000
      로그정규        0.0394          1.0000

from scipy.stats import levene

# Levene 검정으로 두 집단의 분산이 같다고 볼 수 있는지 확인한다.
group1 = data[data['species'] == 'versicolor']['sepal_length']
group2 = data[data['species'] == 'virginica']['sepal_length']
stat, p_value = levene(group1, group2)
print(f"Levene's Test: F = {stat:.4f}, p-value = {p_value:.4f}")

Levene's Test: F = 1.0245, p-value = 0.3139

import numpy as np
import seaborn as sns
from scipy import stats

data = sns.load_dataset("iris")
data = data[data["species"] != "setosa"]
g1 = data[data.species == 'versicolor'].sepal_length.values
g2 = data[data.species == 'virginica'].sepal_length.values
s1, s2 = g1.std(ddof=1), g2.std(ddof=1)
print(f"표준편차 {s1:.4f}, {s2:.4f}   비 {s2 / s1:.4f},  분산비 {(s2 / s1) ** 2:.4f}")
print(f"Levene  p = {stats.levene(g1, g2).pvalue:.4f}")
print(f"Bartlett p = {stats.bartlett(g1, g2).pvalue:.4f}")

t_pool = stats.ttest_ind(g1, g2)
t_welch = stats.ttest_ind(g1, g2, equal_var=False)
print(f"\n합동 t = {t_pool.statistic:.6f},  df = 98,       p = {t_pool.pvalue:.3e}")
print(f"Welch t = {t_welch.statistic:.6f},  df = {t_welch.df:.4f},  p = {t_welch.pvalue:.3e}")

rng = np.random.default_rng(55)
B = 10_000
print(f"\n집단당 n = 50, 평균이 같을 때. sigma = (1, r)")
print(f"{'r':>6}{'분산비':>8}{'Levene 적발':>12}{'고전 F 오류율':>14}{'Welch 오류율':>13}")
for r in [1.0, 1.2319, 1.5, 2.0, 3.0]:
    lv = ff = wl = 0
    for _ in range(B):
        a, b = rng.normal(0, 1, 50), rng.normal(0, r, 50)
        lv += stats.levene(a, b).pvalue < 0.05
        ff += stats.f_oneway(a, b).pvalue < 0.05
        wl += stats.ttest_ind(a, b, equal_var=False).pvalue < 0.05
    print(f"{r:>6.4f}{r ** 2:>8.3f}{lv / B:>12.4f}{ff / B:>14.4f}{wl / B:>13.4f}")

표준편차 0.5162, 0.6359   비 1.2319,  분산비 1.5176
Levene  p = 0.3139
Bartlett p = 0.1478

합동 t = -5.629165,  df = 98,       p = 1.725e-07
Welch t = -5.629165,  df = 94.0255,  p = 1.866e-07

집단당 n = 50, 평균이 같을 때. sigma = (1, r)
     r     분산비   Levene 적발      고전 F 오류율    Welch 오류율
1.0000   1.000      0.0460        0.0476       0.0476
1.2319   1.518      0.2552        0.0516       0.0516
1.5000   2.250      0.7217        0.0507       0.0503
2.0000   4.000      0.9924        0.0509       0.0502
3.0000   9.000      1.0000        0.0488       0.0469

# 일원배치에서 적합값은 집단평균뿐이므로 세로줄이 집단 수만큼만 생긴다.
# 각 줄의 퍼짐이 비슷한지를 본다.
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()

import numpy as np
import seaborn as sns
from statsmodels.formula.api import ols

data = sns.load_dataset("iris")
data = data[data["species"] != "setosa"]
model = ols('sepal_length ~ species', data=data).fit()
e = model.resid.values
fit = model.fittedvalues.values
sp = np.array(data.species)

print(f"서로 다른 적합값 = {np.unique(np.round(fit, 6))}")
print(f"\n{'띠':>12}{'합':>11}{'최소':>9}{'최대':>9}{'폭':>8}{'sd(ddof=0)':>12}")
for s in ['versicolor', 'virginica']:
    ei = e[sp == s]
    print(f"{s:>12}{ei.sum():>11.1e}{ei.min():>9.3f}{ei.max():>9.3f}"
          f"{ei.max() - ei.min():>8.3f}{ei.std(ddof=0):>12.4f}")

r = model.get_influence().resid_studentized_internal
print(f"\n표준화 잔차: 최대 |r| = {np.abs(r).max():.4f},  |r|>2 인 개수 = {(np.abs(r) > 2).sum()}")
print(f"관측값 100 개가 갖는 서로 다른 값의 수 = {len(set(data.sepal_length))}")
print(f"  (0.1 눈금으로 기록되어 점이 겹친다)")
print(f"두 종이 겹치는 구간 = [{data[data.species == 'virginica'].sepal_length.min():.1f}, "
      f"{data[data.species == 'versicolor'].sepal_length.max():.1f}]")

서로 다른 적합값 = [5.936 6.588]

           띠          합       최소       최대       폭  sd(ddof=0)
  versicolor    4.6e-14   -1.036    1.064   2.100      0.5110
   virginica    4.4e-14   -1.688    1.312   3.000      0.6295

표준화 잔차: 최대 |r| = 2.9443,  |r|>2 인 개수 = 2
관측값 100 개가 갖는 서로 다른 값의 수 = 28
  (0.1 눈금으로 기록되어 점이 겹친다)
두 종이 겹치는 구간 = [4.9, 7.0]

import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols

# 근무 형태 세 가지에 따른 생산성. 집단마다 다섯 명씩이다.
data = pd.DataFrame({
    'productivity': [68, 75, 80, 65, 85, 78, 70, 82, 90, 88, 72, 95, 67, 85, 79],
    'environment': ['remote']*5 + ['office']*5 + ['hybrid']*5
})

model = ols('productivity ~ environment', data=data).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)

             sum_sq    df         F    PR(>F)
environment   130.0   2.0  0.768019  0.485443
Residual     1015.6  12.0       NaN       NaN

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
from scipy import optimize, stats

data = pd.DataFrame({
    'productivity': [68, 75, 80, 65, 85, 78, 70, 82, 90, 88, 72, 95, 67, 85, 79],
    'environment': ['remote'] * 5 + ['office'] * 5 + ['hybrid'] * 5,
})
n, k = 5, 3
N = n * k
print(data.groupby('environment', sort=False).productivity.agg(['mean', 'std']).round(3))

SSB, SSE = 130.0, 1015.6
MSE = SSE / (N - k)
sigma = np.sqrt(MSE)
print(f"\nMSE = {MSE:.4f},  합동 표준편차 = {sigma:.4f}")
print(f"관측된 평균의 최대-최소 = {81.6 - 74.6:.1f}")
print(f"eta^2 = SSB/SST = {SSB / (SSB + SSE):.4f}")
print(f"omega^2 = (SSB - (k-1)MSE)/(SST + MSE) = "
      f"{(SSB - (k - 1) * MSE) / (SSB + SSE + MSE):.4f}")

crit = stats.f.ppf(0.95, k - 1, N - k)
power = lambda f: stats.ncf.sf(crit, k - 1, N - k, N * f ** 2)
print(f"\n임계값 F_0.95(2,12) = {crit:.4f}")
print(f"평균이 (mu-d/2, mu, mu+d/2) 꼴이면 f = d/(sigma*sqrt(6))")
print(f"{'Cohen f':>9}{'비중심모수':>12}{'검정력':>9}{'끝 두 집단의 차 d':>18}")
for f in [0.25, 0.40, 0.60, 0.80, 1.00]:
    print(f"{f:>9.4f}{N * f ** 2:>12.3f}{power(f):>9.4f}{f * sigma * np.sqrt(6):>18.2f}")

f80 = optimize.brentq(lambda f: power(f) - 0.80, 0.1, 1.5)
print(f"\n검정력 0.80 이 되는 f = {f80:.4f},  그때 d = {f80 * sigma * np.sqrt(6):.2f}")
f_obs = np.sqrt(SSB / SSE)
print(f"관측된 f = sqrt(SSB/SSE) = {f_obs:.4f},  그 크기에서의 검정력 = {power(f_obs):.4f}")

for nn in [5, 10, 20, 40, 80]:
    NN = nn * k
    c = stats.f.ppf(0.95, k - 1, NN - k)
    pw = stats.ncf.sf(c, k - 1, NN - k, NN * f_obs ** 2)
    print(f"  집단당 n = {nn:>3}: 관측 효과크기 f = {f_obs:.4f} 의 검정력 = {pw:.4f}")

             mean     std
environment              
remote       74.6   8.264
office       81.6   8.050
hybrid       79.6  10.991

MSE = 84.6333,  합동 표준편차 = 9.1996
관측된 평균의 최대-최소 = 7.0
eta^2 = SSB/SST = 0.1135
omega^2 = (SSB - (k-1)MSE)/(SST + MSE) = -0.0319

임계값 F_0.95(2,12) = 3.8853
평균이 (mu-d/2, mu, mu+d/2) 꼴이면 f = d/(sigma*sqrt(6))
  Cohen f       비중심모수      검정력       끝 두 집단의 차 d
   0.2500       0.938   0.1095              5.63
   0.4000       2.400   0.2137              9.01
   0.6000       5.400   0.4334             13.52
   0.8000       9.600   0.6827             18.03
   1.0000      15.000   0.8696             22.53

검정력 0.80 이 되는 f = 0.9130,  그때 d = 20.57
관측된 f = sqrt(SSB/SSE) = 0.3578,  그 크기에서의 검정력 = 0.1786
  집단당 n =   5: 관측 효과크기 f = 0.3578 의 검정력 = 0.1786
  집단당 n =  10: 관측 효과크기 f = 0.3578 의 검정력 = 0.3630
  집단당 n =  20: 관측 효과크기 f = 0.3578 의 검정력 = 0.6741
  집단당 n =  40: 관측 효과크기 f = 0.3578 의 검정력 = 0.9436
  집단당 n =  80: 관측 효과크기 f = 0.3578 의 검정력 = 0.9993

import matplotlib.pyplot as plt
from scipy.stats import shapiro, levene

# 정규성 — 잔차의 Q-Q 그림과 Shapiro-Wilk 검정
sm.qqplot(model.resid, line='s')
plt.title("Q-Q Plot of Residuals")
plt.show()

stat, p_value = shapiro(model.resid)
print(f"Shapiro-Wilk Test: p-value = {p_value:.4f}")

# 등분산성 — Levene 검정
group1 = data[data['environment'] == 'remote']['productivity']
group2 = data[data['environment'] == 'office']['productivity']
group3 = data[data['environment'] == 'hybrid']['productivity']
stat, p_value = levene(group1, group2, group3)
print(f"Levene's Test: p-value = {p_value:.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()

Shapiro-Wilk Test: p-value = 0.7449
Levene's Test: p-value = 0.7631

import numpy as np
from scipy import stats

rng = np.random.default_rng(7)
B = 20_000
n, k = 5, 3

print(f"집단당 n = {n}, k = {k} (N = {n * k}).  명목 0.05 에서 기각률")
print(f"\nShapiro-Wilk 를 잔차 {n * k} 개에 적용")
for lab, f in [("정규 (귀무가 참)", lambda m: rng.normal(0, 1, m)),
               ("t(3)", lambda m: rng.standard_t(3, m)),
               ("지수", lambda m: rng.exponential(1, m)),
               ("로그정규", lambda m: np.exp(rng.normal(0, 1, m)))]:
    c = 0
    for _ in range(B):
        gs = [f(n) for _ in range(k)]
        r = np.concatenate([g - g.mean() for g in gs])
        c += stats.shapiro(r).pvalue < 0.05
    print(f"  {lab:>16}: {c / B:.4f}")

print(f"\nLevene(중앙값) 을 정규자료에 적용, sigma = (1, 1, r)")
for r in [1.0, 1.5, 2.0, 3.0, 5.0]:
    c = sum(stats.levene(rng.normal(0, 1, n), rng.normal(0, 1, n),
                         rng.normal(0, r, n)).pvalue < 0.05 for _ in range(B))
    print(f"  sigma 비 r = {r:>4}: {c / B:.4f}")

집단당 n = 5, k = 3 (N = 15).  명목 0.05 에서 기각률

Shapiro-Wilk 를 잔차 15 개에 적용
        정규 (귀무가 참): 0.0410
              t(3): 0.1759
                지수: 0.3719
              로그정규: 0.5527

Levene(중앙값) 을 정규자료에 적용, sigma = (1, 1, r)
  sigma 비 r =  1.0: 0.0052
  sigma 비 r =  1.5: 0.0132
  sigma 비 r =  2.0: 0.0408
  sigma 비 r =  3.0: 0.1338
  sigma 비 r =  5.0: 0.3220

import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols

# 지점 네 곳의 고객만족도. 지점마다 다섯 건씩이다.
data = pd.DataFrame({
    'satisfaction': [4.5, 3.8, 4.7, 4.2, 4.9, 4.1, 3.5, 4.3, 4.8, 3.9,
                     4.4, 4.0, 3.7, 4.2, 4.6, 4.8, 3.6, 4.3, 4.1, 4.7],
    'location': ['A']*5 + ['B']*5 + ['C']*5 + ['D']*5
})

model = ols('satisfaction ~ location', data=data).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)

          sum_sq    df         F    PR(>F)
location  0.2655   3.0  0.456186  0.716615
Residual  3.1040  16.0       NaN       NaN

import numpy as np
import pandas as pd
from scipy import stats

data = pd.DataFrame({
    'satisfaction': [4.5, 3.8, 4.7, 4.2, 4.9, 4.1, 3.5, 4.3, 4.8, 3.9,
                     4.4, 4.0, 3.7, 4.2, 4.6, 4.8, 3.6, 4.3, 4.1, 4.7],
    'location': ['A'] * 5 + ['B'] * 5 + ['C'] * 5 + ['D'] * 5,
})
n, k = 5, 4
N = n * k
g = data.groupby('location').satisfaction
print(g.agg(['mean', 'std']).round(4))

SSB = n * ((g.mean().values - data.satisfaction.mean()) ** 2).sum()
SSE = sum(((data[data.location == l].satisfaction
            - data[data.location == l].satisfaction.mean()) ** 2).sum()
          for l in 'ABCD')
MSE = SSE / (N - k)
F = (SSB / (k - 1)) / MSE
print(f"\nSSB = {SSB:.4f},  SSE = {SSE:.4f},  MSE = {MSE:.6f},  F = {F:.6f}")
print(f"귀무가설 아래 SSB 의 기댓값 = (k-1)*sigma^2 = {(k - 1) * MSE:.4f}")
print(f"관측된 SSB 는 그 {SSB / ((k - 1) * MSE):.2f} 배")
print(f"\neta^2  = SSB/SST                 = {SSB / (SSB + SSE):.4f}")
print(f"omega^2 = (SSB-(k-1)MSE)/(SST+MSE) = {(SSB - (k - 1) * MSE) / (SSB + SSE + MSE):.4f}")
print(f"\nP(F(3,16) <= {F:.4f}) = {stats.f.cdf(F, k - 1, N - k):.4f}")
print(f"곧 귀무가설이 참이어도 이만큼 작은 F 가 나올 확률이 "
      f"{stats.f.cdf(F, k - 1, N - k):.1%} 다")
print(f"평균의 최대-최소 = {g.mean().max() - g.mean().min():.2f},  합동 표준편차 = {np.sqrt(MSE):.4f}")
print(f"관측 효과크기 f = sqrt(SSB/SSE) = {np.sqrt(SSB / SSE):.4f}")
crit = stats.f.ppf(0.95, k - 1, N - k)
print(f"검정력 (관측 f 가 참일 때) = {stats.ncf.sf(crit, k - 1, N - k, N * SSB / SSE):.4f}")

          mean     std
location              
A         4.42  0.4324
B         4.12  0.4817
C         4.18  0.3493
D         4.30  0.4848

SSB = 0.2655,  SSE = 3.1040,  MSE = 0.194000,  F = 0.456186
귀무가설 아래 SSB 의 기댓값 = (k-1)*sigma^2 = 0.5820
관측된 SSB 는 그 0.46 배

eta^2  = SSB/SST                 = 0.0788
omega^2 = (SSB-(k-1)MSE)/(SST+MSE) = -0.0888

P(F(3,16) <= 0.4562) = 0.2834
곧 귀무가설이 참이어도 이만큼 작은 F 가 나올 확률이 28.3% 다
평균의 최대-최소 = 0.30,  합동 표준편차 = 0.4405
관측 효과크기 f = sqrt(SSB/SSE) = 0.2925
검정력 (관측 f 가 참일 때) = 0.1419

import matplotlib.pyplot as plt
from scipy.stats import shapiro, levene

# 정규성 — 잔차의 Q-Q 그림과 Shapiro-Wilk 검정
sm.qqplot(model.resid, line='s')
plt.title("Q-Q Plot of Residuals")
plt.show()

stat, p_value = shapiro(model.resid)
print(f"Shapiro-Wilk Test: p-value = {p_value:.4f}")

# 등분산성 — Levene 검정
groups = [data[data['location'] == loc]['satisfaction'] for loc in ['A', 'B', 'C', 'D']]
stat, p_value = levene(*groups)
print(f"Levene's Test: p-value = {p_value:.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()

Shapiro-Wilk Test: p-value = 0.5488
Levene's Test: p-value = 0.9343

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(99)
B = 10_000
n, k = 5, 4

def nrm(s):
    return lambda: rng.normal(0, s, n)

def lgn(s):
    return lambda: s * np.exp(rng.normal(0, 1, n))

scen = {
    "가정 모두 성립": [nrm(1)] * 4,
    "비정규만 (로그정규)": [lgn(1)] * 4,
    "이분산만 (1,1,1,3)": [nrm(1), nrm(1), nrm(1), nrm(3)],
    "둘 다": [lgn(1)] * 3 + [lgn(3)],
}

print(f"집단당 n = {n}, k = {k}.  평균은 언제나 같다.  명목 0.05")
print(f"{'상황':>20}{'Shapiro':>9}{'Levene':>8}{'둘 중 하나':>11}"
      f"{'고전 F':>8}{'Welch 류':>10}")
for lab, fl in scen.items():
    a = b = c = d = e = 0
    for _ in range(B):
        gs = [f() for f in fl]
        r = np.concatenate([g - g.mean() for g in gs])
        s = stats.shapiro(r).pvalue < 0.05
        l = stats.levene(*gs).pvalue < 0.05
        a += s
        b += l
        c += (s or l)
        d += stats.f_oneway(*gs).pvalue < 0.05
        e += stats.alexandergovern(*gs).pvalue < 0.05
    print(f"{lab:>20}{a / B:>9.4f}{b / B:>8.4f}{c / B:>11.4f}{d / B:>8.4f}{e / B:>10.4f}")

집단당 n = 5, k = 4.  평균은 언제나 같다.  명목 0.05
                  상황  Shapiro  Levene     둘 중 하나    고전 F   Welch 류
            가정 모두 성립   0.0437  0.0036     0.0472  0.0474    0.0463
         비정규만 (로그정규)   0.6747  0.0205     0.6878  0.0320    0.0361
      이분산만 (1,1,1,3)   0.3043  0.1749     0.4092  0.0923    0.0468
                 둘 다   0.7061  0.0831     0.7352  0.2468    0.0953

from statsmodels.stats.multicomp import pairwise_tukeyhsd

# 분산분석이 유의했으므로 어느 지점 쌍이 다른지 사후비교로 좁힌다.
tukey = pairwise_tukeyhsd(data['satisfaction'], data['location'], alpha=0.05)
print(tukey)

Multiple Comparison of Means - Tukey HSD, FWER=0.05
=================================================
group1 group2 meandiff p-adj  lower  upper reject
-------------------------------------------------
     A      B     -0.3  0.708 -1.097 0.497  False
     A      C    -0.24 0.8243 -1.037 0.557  False
     A      D    -0.12 0.9723 -0.917 0.677  False
     B      C     0.06 0.9963 -0.737 0.857  False
     B      D     0.18 0.9154 -0.617 0.977  False
     C      D     0.12 0.9723 -0.677 0.917  False
-------------------------------------------------

import numpy as np
from scipy import stats

rng = np.random.default_rng(314)
B = 20_000
n, k = 5, 4
N, nu, m = n * k, n * k - k, k * (k - 1) // 2
q_c = stats.studentized_range.ppf(0.95, k, nu)
t_c = stats.t.ppf(0.975, nu)
F_c = stats.f.ppf(0.95, k - 1, nu)
print(f"k = {k}, 집단당 n = {n}, nu = {nu}, 쌍 수 m = {m}")
print(f"임계값: q(0.95,4,16)/sqrt2 = {q_c / np.sqrt(2):.4f},  t(0.975,16) = {t_c:.4f},  "
      f"F(0.95,3,16) = {F_c:.4f}")

cnt = dict(tukey=0, tukey_gated=0, lsd=0, lsd_gated=0, anova=0, tukey_only=0)
for _ in range(B):
    gs = [rng.normal(0, 1, n) for _ in range(k)]
    mi = np.array([g.mean() for g in gs])
    MSE = np.mean([g.var(ddof=1) for g in gs])
    SE = np.sqrt(2 * MSE / n)
    tmax = max(abs(mi[i] - mi[j]) for i in range(k) for j in range(i + 1, k)) / SE
    Fv = (n * ((mi - mi.mean()) ** 2).sum() / (k - 1)) / MSE
    sig_F = Fv > F_c
    sig_T = tmax > q_c / np.sqrt(2)
    sig_L = tmax > t_c
    cnt['anova'] += sig_F
    cnt['tukey'] += sig_T
    cnt['tukey_gated'] += (sig_F and sig_T)
    cnt['lsd'] += sig_L
    cnt['lsd_gated'] += (sig_F and sig_L)
    cnt['tukey_only'] += (sig_T and not sig_F)

print(f"\n평균이 모두 같을 때 (귀무가 참), 명목 0.05")
for lab, key in [("전역 F 만", 'anova'), ("투키 (전역검정 없이)", 'tukey'),
                 ("투키 (전역검정 통과 뒤에만)", 'tukey_gated'),
                 ("LSD (전역검정 없이)", 'lsd'), ("LSD (전역검정 통과 뒤에만)", 'lsd_gated')]:
    print(f"  {lab:>26}: {cnt[key] / B:.4f}")
print(f"\n전역 F 는 기각하지 못했는데 투키가 어떤 쌍을 기각한 비율 = {cnt['tukey_only'] / B:.4f}")

k = 4, 집단당 n = 5, nu = 16, 쌍 수 m = 6
임계값: q(0.95,4,16)/sqrt2 = 2.8610,  t(0.975,16) = 2.1199,  F(0.95,3,16) = 3.2389

평균이 모두 같을 때 (귀무가 참), 명목 0.05
                      전역 F 만: 0.0501
                투키 (전역검정 없이): 0.0510
            투키 (전역검정 통과 뒤에만): 0.0437
               LSD (전역검정 없이): 0.1916
           LSD (전역검정 통과 뒤에만): 0.0501

전역 F 는 기각하지 못했는데 투키가 어떤 쌍을 기각한 비율 = 0.0073

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.libqsturng import qsturng
from sklearn.datasets import load_iris

d = load_iris()
df = pd.DataFrame(d.data, columns=["sl", "sw", "pl", "pw"])
df["sp"] = [d.target_names[i] for i in d.target]

print("사례 1: 붓꽃 꽃받침 너비(sepal width)")
g = df.groupby("sp").sw.agg(["size", "mean", "std"]).round(4)
print(g.to_string())

gs = [v.sw.values for _, v in df.groupby("sp")]
lev = stats.levene(*gs, center="median")
print(f"\n브라운-포사이드: W = {lev.statistic:.4f}, p = {lev.pvalue:.4f}")
r = stats.f_oneway(*gs)
print(f"표준 F: {r.statistic:.4f}, p = {r.pvalue:.3e}")

def welch(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))
    df2 = (3 / (k**2 - 1) * lam)**-1
    return F, df2, stats.f.sf(F, k - 1, df2)

F, d2, p = welch(gs)
print(f"웰치 F = {F:.4f}, df2 = {d2:.3f}, p = {p:.3e}")

N, k = len(df), 3
print(f"η² = {2 * r.statistic / (2 * r.statistic + N - k):.4f},  "
      f"ω² = {2 * (r.statistic - 1) / (2 * r.statistic + N - k + 1):.4f}")

print("\n게임스-하웰")
names = list(g.index)
n = np.array([len(x) for x in gs], float)
m = np.array([x.mean() for x in gs])
v = np.array([x.var(ddof=1) for x in gs])
for i in range(3):
    for j in range(i + 1, 3):
        s = v[i] / n[i] + v[j] / n[j]
        dfw = s**2 / (v[i]**2 / (n[i]**2 * (n[i] - 1))
                      + v[j]**2 / (n[j]**2 * (n[j] - 1)))
        diff = m[j] - m[i]
        h = qsturng(0.95, 3, dfw) * np.sqrt(s / 2)
        print(f"  {names[j]:>12s} - {names[i]:<12s} {diff:+7.4f}  "
              f"95% [{diff - h:+7.4f},{diff + h:+7.4f}]  df={dfw:6.2f}")

사례 1: 붓꽃 꽃받침 너비(sepal width)
            size   mean     std
sp                             
setosa        50  3.428  0.3791
versicolor    50  2.770  0.3138
virginica     50  2.974  0.3225

브라운-포사이드: W = 0.5902, p = 0.5555
표준 F: 49.1600, p = 4.492e-17
웰치 F = 45.0120, df2 = 97.402, p = 1.433e-14
η² = 0.4008,  ω² = 0.3910

게임스-하웰
    versicolor - setosa       -0.6580  95% [-0.8237,-0.4923]  df= 94.70
     virginica - setosa       -0.4540  95% [-0.6216,-0.2864]  df= 95.55
     virginica - versicolor   +0.2040  95% [+0.0525,+0.3555]  df= 97.93

붓꽃 세 종(각 n = 50)의 꽃받침 너비를 비교했다.

기술통계
  setosa      3.428 ± 0.379
  versicolor  2.770 ± 0.314
  virginica   2.974 ± 0.323

가정 점검
  브라운-포사이드 등분산 검정  W = 0.59, p = 0.56  (문제 없음)
  잔차 Q-Q 그림에서 뚜렷한 이탈 없음

분석
  Welch F(2, 97.4) = 45.01,  p < 0.001,  ω² = 0.391
  사후검정: 게임스-하웰 (FWER = 0.05)

  versicolor − setosa     −0.658  [−0.824, −0.492]
  virginica  − setosa     −0.454  [−0.622, −0.286]
  virginica  − versicolor +0.204  [ 0.053,  0.356]

결론
  세 종의 꽃받침 너비가 모두 서로 다르다. setosa 가 가장 넓고,
  versicolor 가 가장 좁다. 종이 변동의 39% 를 설명한다.

관찰된 것:  재택 근무자의 생산성이 사무실 근무자보다 높다

가능한 설명 넷
  (가) 재택 근무가 생산성을 높인다            ← 우리가 원하는 해석
  (나) 생산성 높은 사람이 재택을 선택했다      ← 자기 선택
  (다) 생산성 높은 직무가 재택 가능하다        ← 직무 교란
  (라) 재택 근무자의 생산성 측정이 다르다      ← 측정 편향

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(19001)
B = 3_000
print("근무 형태의 참 효과는 0. 다만 '능력'이 선택과 생산성에 모두 영향")
print(f"{'선택 강도 (능력→재택 확률)':>26s} {'관측 차이 평균':>13s} {'유의 비율':>9s}")
for gamma in [0.0, 0.5, 1.0, 2.0]:
    diffs, hit = [], 0
    for _ in range(B):
        N = 150
        ability = rng.normal(0, 1, N)
        p = 1 / (1 + np.exp(-gamma * ability))       # 능력이 높을수록 재택
        home = rng.random(N) < p
        prod = 50 + 5 * ability + rng.normal(0, 3, N)  # 근무 형태의 효과는 0
        if home.sum() < 5 or (~home).sum() < 5:
            continue
        d = prod[home].mean() - prod[~home].mean()
        diffs.append(d)
        hit += stats.ttest_ind(prod[home], prod[~home],
                               equal_var=False).pvalue < 0.05
    print(f"{gamma:26.1f} {np.mean(diffs):13.3f} {hit / len(diffs):9.4f}")

근무 형태의 참 효과는 0. 다만 '능력'이 선택과 생산성에 모두 영향
          선택 강도 (능력→재택 확률)      관측 차이 평균     유의 비율
                       0.0         0.023    0.0410
                       0.5         2.354    0.7060
                       1.0         4.124    0.9943
                       2.0         6.052    1.0000

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(18001)
B = 6_000
doses = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
k, n = len(doses), 10
N, dfe = k * n, k * n - k

lin = doses - doses.mean()
lin /= np.sqrt((lin**2).sum())
quad = (doses - doses.mean())**2
quad -= quad.mean()
quad /= np.sqrt((quad**2).sum())

print("용량당 n=10, 오차 SD=1.0, 명목 0.05")
print(f"{'참 반응 모양':>22s} {'전체 F':>8s} {'선형 대비':>9s} {'이차 대비':>9s}")
for lab, mus in [("평평 (귀무)", np.zeros(k)),
                 ("선형 증가 0.3/단위", 0.3 * doses),
                 ("포화 (로그형)", 1.2 * np.log(doses + 1)),
                 ("역U자", -0.35 * (doses - 2)**2 + 1.4),
                 ("마지막만 점프", np.array([0, 0, 0, 0, 1.5]))]:
    a = b = c = 0
    for _ in range(B):
        gs = [rng.normal(mu, 1.0, n) for mu in mus]
        m = np.array([g.mean() for g in gs])
        MSE = np.array([g.var(ddof=1) for g in gs]).mean()
        a += stats.f_oneway(*gs).pvalue < 0.05
        for cvec, which in [(lin, "b"), (quad, "c")]:
            se = np.sqrt(MSE * (cvec**2).sum() / n)
            p = 2 * stats.t.sf(abs(cvec @ m / se), dfe)
            if which == "b":
                b += p < 0.05
            else:
                c += p < 0.05
    print(f"{lab:>22s} {a / B:8.4f} {b / B:9.4f} {c / B:9.4f}")

용량당 n=10, 오차 SD=1.0, 명목 0.05
               참 반응 모양     전체 F     선형 대비     이차 대비
               평평 (귀무)   0.0492    0.0498    0.0443
          선형 증가 0.3/단위   0.6067    0.8287    0.0503
              포화 (로그형)   0.9717    0.9973    0.1843
                   역U자   0.8987    0.0498    0.9815
               마지막만 점프   0.9188    0.8345    0.6983

학교 20곳 → 프로그램 A(10곳) / 프로그램 B(10곳)
각 학교에서 학생 30명 측정 → 총 600명

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(19002)
B = 3_000
n_school, n_stud = 10, 30      # 프로그램당 학교 수, 학교당 학생 수

print("프로그램의 참 효과는 0. 학교 효과(ICC)만 존재. 명목 0.05")
print(f"{'ICC':>6s} {'학생 단위 t 검정':>14s} {'학교 평균 t 검정':>15s} "
      f"{'설계효과':>9s}")
for icc in [0.0, 0.05, 0.10, 0.20]:
    a = b = 0
    for _ in range(B):
        sa, sb = [], []
        for _ in range(n_school):
            u = rng.normal(0, np.sqrt(icc))
            sa.append(u + rng.normal(0, np.sqrt(1 - icc), n_stud))
            u = rng.normal(0, np.sqrt(icc))
            sb.append(u + rng.normal(0, np.sqrt(1 - icc), n_stud))
        a += stats.ttest_ind(np.concatenate(sa),
                             np.concatenate(sb)).pvalue < 0.05
        b += stats.ttest_ind([x.mean() for x in sa],
                             [x.mean() for x in sb]).pvalue < 0.05
    print(f"{icc:6.2f} {a / B:14.4f} {b / B:15.4f} "
          f"{1 + (n_stud - 1) * icc:9.2f}")

프로그램의 참 효과는 0. 학교 효과(ICC)만 존재. 명목 0.05
   ICC     학생 단위 t 검정      학교 평균 t 검정      설계효과
  0.00         0.0500          0.0537      1.00
  0.05         0.2143          0.0523      2.45
  0.10         0.3173          0.0497      3.90
  0.20         0.4580          0.0477      6.80

설계: 군집 무작위 배정 (학교 20곳, 학생 600명)
분석 단위: 학교 (배정 단위와 일치)

  프로그램 A  학교 10곳, 학교 평균 점수 72.4 ± 4.1
  프로그램 B  학교 10곳, 학교 평균 점수 76.8 ± 3.9

  t(18) = 2.46,  p = 0.024,  차이 +4.4 [0.6, 8.2]

급내상관 ICC = 0.12 (설계효과 4.5). 학생 600명은 실질적으로
독립 관측 133명에 해당한다.

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.libqsturng import qsturng

def anova_case_study(groups, names=None, alpha=0.05):
    """일원배치 분산분석의 전 과정을 수행하고 보고서를 출력한다."""
    k = len(groups)
    names = names or [f"G{i}" for i in range(k)]
    n = np.array([len(g) for g in groups], float)
    m = np.array([g.mean() for g in groups])
    v = np.array([g.var(ddof=1) for g in groups])
    res = np.concatenate([g - g.mean() for g in groups])

    print("[1] 기술통계")
    for i in range(k):
        print(f"    {names[i]:>12s}  n={int(n[i]):4d}  "
              f"평균={m[i]:9.4f}  s={np.sqrt(v[i]):8.4f}")
    s = np.sqrt(v)
    print(f"    SD 비 = {s.max() / s.min():.2f},  "
          f"n 비 = {n.max() / n.min():.2f}")

    print("\n[2] 가정 점검")
    lev = stats.levene(*groups, center="median")
    sw = stats.shapiro(res)
    print(f"    브라운-포사이드  W = {lev.statistic:8.4f}  p = {lev.pvalue:.4f}")
    print(f"    샤피로(잔차)     W = {sw.statistic:8.4f}  p = {sw.pvalue:.4f}")
    print(f"    잔차 치우침 = {stats.skew(res):+.3f}, "
          f"초과첨도 = {stats.kurtosis(res):+.3f}")

    print("\n[3] 전체 검정")
    r = stats.f_oneway(*groups)
    N = int(n.sum())
    w = n / v
    W = w.sum()
    mt = (w * m).sum() / W
    lam = ((1 - w / W)**2 / (n - 1)).sum()
    Fw = ((w * (m - mt)**2).sum() / (k - 1)
          / (1 + 2 * (k - 2) / (k**2 - 1) * lam))
    df2 = (3 / (k**2 - 1) * lam)**-1
    print(f"    표준 F({k - 1}, {N - k}) = {r.statistic:.4f}, p = {r.pvalue:.3e}")
    print(f"    웰치  F({k - 1}, {df2:.1f}) = {Fw:.4f}, "
          f"p = {stats.f.sf(Fw, k - 1, df2):.3e}")

    print("\n[4] 효과크기")
    e2 = 2 * r.statistic / (2 * r.statistic + N - k) if k == 3 else \
        (k - 1) * r.statistic / ((k - 1) * r.statistic + N - k)
    w2 = ((k - 1) * (r.statistic - 1)
          / ((k - 1) * r.statistic + N - k + 1))
    print(f"    η² = {e2:.4f},  ω² = {w2:.4f},  "
          f"Cohen f = {np.sqrt(max(w2, 0) / max(1 - w2, 1e-9)):.4f}")

    print("\n[5] 사후비교 (게임스-하웰, 등분산 가정 없음)")
    for i in range(k):
        for j in range(i + 1, k):
            sij = v[i] / n[i] + v[j] / n[j]
            dfw = sij**2 / (v[i]**2 / (n[i]**2 * (n[i] - 1))
                            + v[j]**2 / (n[j]**2 * (n[j] - 1)))
            diff = m[j] - m[i]
            h = qsturng(1 - alpha, k, dfw) * np.sqrt(sij / 2)
            mark = "*" if abs(diff) > h else " "
            print(f"    {names[j]:>12s} - {names[i]:<12s} {diff:+8.4f}  "
                  f"95% [{diff - h:+8.4f},{diff + h:+8.4f}] {mark}")

rng = np.random.default_rng(19003)
gs = [rng.normal(mu, sd, n) for n, mu, sd in
      [(24, 72.0, 6.0), (30, 76.5, 8.5), (18, 74.0, 5.0)]]
anova_case_study(gs, ["대조", "처치A", "처치B"])