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

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/ToothGrowth.csv')
df = pd.read_csv(url, usecols=[1, 2, 3])

model = ols('len ~ C(supp) + C(dose) + C(supp):C(dose)', data=df).fit()
# typ=2를 명시한다. statsmodels의 기본값은 typ=1(순차적 제곱합)이라
# 모형에 넣는 항의 **순서에 따라 결과가 달라진다**. 균형 설계에서는
# 세 유형이 모두 같지만, 불균형이면 갈린다.
aov2 = anova_lm(model, typ=2)
print(aov2)

                      sum_sq    df          F        PR(>F)
C(supp)           205.350000   1.0  15.571979  2.311828e-04
C(dose)          2426.434333   2.0  91.999965  4.046291e-18
C(supp):C(dose)   108.319000   2.0   4.106991  2.186027e-02
Residual          712.106000  54.0        NaN           NaN

# 균형 설계인가. 여섯 칸의 개수를 센다.
print("칸별 관측 수")
print(df.groupby(['supp', 'dose']).size().unstack())

a1 = anova_lm(model, typ=1)['sum_sq']
model_rev = ols('len ~ C(dose) + C(supp) + C(dose):C(supp)', data=df).fit()
a1r = anova_lm(model_rev, typ=1)['sum_sq']
a2 = anova_lm(model, typ=2)['sum_sq']
a3_bad = anova_lm(model, typ=3)['sum_sq']
# 제III형은 합이 0 이 되는 코딩(Sum)으로 적합해야 올바른 값이 나온다.
model_sum = ols('len ~ C(supp, Sum) + C(dose, Sum) + C(supp, Sum):C(dose, Sum)',
                data=df).fit()
a3 = anova_lm(model_sum, typ=3)['sum_sq']

print(f"\n{'방법':<22}{'supp':>12}{'dose':>14}{'supp:dose':>12}")
print(f"{'제I형 (supp 먼저)':<22}{a1['C(supp)']:>12.4f}{a1['C(dose)']:>14.4f}"
      f"{a1['C(supp):C(dose)']:>12.4f}")
print(f"{'제I형 (dose 먼저)':<22}{a1r['C(supp)']:>12.4f}{a1r['C(dose)']:>14.4f}"
      f"{a1r['C(dose):C(supp)']:>12.4f}")
print(f"{'제II형':<22}{a2['C(supp)']:>12.4f}{a2['C(dose)']:>14.4f}"
      f"{a2['C(supp):C(dose)']:>12.4f}")
print(f"{'제III형 (처리코딩)':<22}{a3_bad['C(supp)']:>12.4f}{a3_bad['C(dose)']:>14.4f}"
      f"{a3_bad['C(supp):C(dose)']:>12.4f}")
print(f"{'제III형 (합코딩)':<22}{a3['C(supp, Sum)']:>12.4f}{a3['C(dose, Sum)']:>14.4f}"
      f"{a3['C(supp, Sum):C(dose, Sum)']:>12.4f}")

print(f"\n네 제곱합의 합 = {a2.sum():.4f},  "
      f"SST = {((df['len'] - df['len'].mean()) ** 2).sum():.4f}")
print(f"자유도 59 = 1 + 2 + 2 + 54")

칸별 관측 수
dose  0.5  1.0  2.0
supp
OJ     10   10   10
VC     10   10   10

방법                            supp          dose   supp:dose
제I형 (supp 먼저)             205.3500     2426.4343    108.3190
제I형 (dose 먼저)             205.3500     2426.4343    108.3190
제II형                      205.3500     2426.4343    108.3190
제III형 (처리코딩)              137.8125      885.2647    108.3190
제III형 (합코딩)               205.3500     2426.4343    108.3190

네 제곱합의 합 = 3452.2093,  SST = 3452.2093
자유도 59 = 1 + 2 + 2 + 54

# OJ_0.5 에서 둘, VC_0.5 에서 넷을 지워 불균형을 만든다.
drop = df.index[(df['supp'] == 'OJ') & (df['dose'] == 0.5)][:2].tolist()
drop += df.index[(df['supp'] == 'VC') & (df['dose'] == 0.5)][:4].tolist()
dfu = df.drop(index=drop).reset_index(drop=True)
print("칸별 관측 수 (불균형)")
print(dfu.groupby(['supp', 'dose']).size().unstack())

mu = ols('len ~ C(supp) + C(dose) + C(supp):C(dose)', data=dfu).fit()
mur = ols('len ~ C(dose) + C(supp) + C(dose):C(supp)', data=dfu).fit()
u1, u1r = anova_lm(mu, typ=1)['sum_sq'], anova_lm(mur, typ=1)['sum_sq']
u2 = anova_lm(mu, typ=2)['sum_sq']
print(f"\n{'방법':<22}{'supp':>12}{'dose':>14}")
print(f"{'제I형 (supp 먼저)':<22}{u1['C(supp)']:>12.4f}{u1['C(dose)']:>14.4f}")
print(f"{'제I형 (dose 먼저)':<22}{u1r['C(supp)']:>12.4f}{u1r['C(dose)']:>14.4f}")
print(f"{'제II형':<22}{u2['C(supp)']:>12.4f}{u2['C(dose)']:>14.4f}")
print(f"\n제II형의 supp 가 '제I형에서 supp 를 나중에 넣은 값'과 같은가: "
      f"{abs(u2['C(supp)'] - u1r['C(supp)']) < 1e-9}")
print(f"제II형의 dose 가 '제I형에서 dose 를 나중에 넣은 값'과 같은가: "
      f"{abs(u2['C(dose)'] - u1['C(dose)']) < 1e-9}")

칸별 관측 수 (불균형)
dose  0.5  1.0  2.0
supp
OJ      8   10   10
VC      6   10   10

방법                            supp          dose
제I형 (supp 먼저)              75.2894     2061.9230
제I형 (dose 먼저)             125.6611     2011.5513
제II형                      125.6611     2061.9230

제II형의 supp 가 '제I형에서 supp 를 나중에 넣은 값'과 같은가: True
제II형의 dose 가 '제I형에서 dose 를 나중에 넣은 값'과 같은가: True

from statsmodels.stats.multicomp import pairwise_tukeyhsd

# 두 요인의 주효과를 각각 사후비교한다. 다른 요인은 잠시 무시하는 셈이다.
print(pairwise_tukeyhsd(endog=df['len'], groups=df['dose'], alpha=0.05))
print(pairwise_tukeyhsd(endog=df['len'], groups=df['supp'], alpha=0.05))

Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj  lower   upper  reject
---------------------------------------------------
   0.5    1.0     9.13   0.0  5.9018 12.3582   True
   0.5    2.0   15.495   0.0 12.2668 18.7232   True
   1.0    2.0    6.365   0.0  3.1368  9.5932   True
---------------------------------------------------
Multiple Comparison of Means - Tukey HSD, FWER=0.05
=================================================
group1 group2 meandiff p-adj  lower  upper reject
-------------------------------------------------
    OJ     VC     -3.7 0.0604 -7.567 0.167  False
-------------------------------------------------

import numpy as np
from scipy import stats

SS_supp, SS_dose, SS_int, SSE2 = 205.350, 2426.434333, 108.319, 712.106

# supp 만 보는 일원배치의 분모를 손으로 만든다.
SSW1 = SS_dose + SS_int + SSE2
MSE1, MSE2 = SSW1 / 58, SSE2 / 54
print(f"SSW(1way) = {SS_dose:.3f} + {SS_int:.3f} + {SSE2:.3f} = {SSW1:.3f}  (df 58)")
print(f"MSE(1way) = {MSE1:.4f}   MSE(2way) = {MSE2:.4f}   비 = {MSE1 / MSE2:.4f}")
print(f"F(1way) = {SS_supp / MSE1:.4f}   F(2way) = {SS_supp / MSE2:.4f}")

oj = df.loc[df['supp'] == 'OJ', 'len'].values
vc = df.loc[df['supp'] == 'VC', 'len'].values
F1, p1 = stats.f_oneway(oj, vc)
print(f"\nf_oneway(OJ, VC):  F = {F1:.4f},  p = {p1:.6f}")

d = vc.mean() - oj.mean()
q_obs = abs(d) / np.sqrt(MSE1 / 30)
p_tukey = stats.studentized_range.sf(q_obs, 2, 58)
HSD = stats.studentized_range.ppf(0.95, 2, 58) * np.sqrt(MSE1 / 30)
print(f"Tukey(supp):  meandiff = {d:.2f},  p-adj = {p_tukey:.6f}")
print(f"              HSD = {HSD:.4f}  ->  CI = ({d - HSD:.3f}, {d + HSD:.3f})")

# dose 주효과 쪽도 같은 방식으로 되살려 본다.
SSW1d = SS_supp + SS_int + SSE2
HSD_dose = stats.studentized_range.ppf(0.95, 3, 57) * np.sqrt((SSW1d / 57) / 20)
print(f"\nTukey(dose):  MSE(1way) = {SSW1d / 57:.4f},  HSD = {HSD_dose:.4f}")

SSW(1way) = 2426.434 + 108.319 + 712.106 = 3246.859  (df 58)
MSE(1way) = 55.9803   MSE(2way) = 13.1871   비 = 4.2451
F(1way) = 3.6683   F(2way) = 15.5720

f_oneway(OJ, VC):  F = 3.6683,  p = 0.060393
Tukey(supp):  meandiff = -3.70,  p-adj = 0.060393
              HSD = 3.8670  ->  CI = (-7.567, 0.167)

Tukey(dose):  MSE(1way) = 17.9961,  HSD = 3.2282

# 교호작용이 유의하면 주효과만으로는 부족하다. 두 요인을 붙여 만든 여섯 칸을
# 서로 견주어야 "어느 조합이 어느 조합과 다른가"를 말할 수 있다.
df['supp_dose'] = df['supp'].astype(str) + "_" + df['dose'].astype(str)
print(pairwise_tukeyhsd(endog=df['len'], groups=df['supp_dose'], alpha=0.05))

 Multiple Comparison of Means - Tukey HSD, FWER=0.05
======================================================
group1 group2 meandiff p-adj   lower    upper   reject
------------------------------------------------------
OJ_0.5 OJ_1.0     9.47    0.0   4.6719  14.2681   True
OJ_0.5 OJ_2.0    12.83    0.0   8.0319  17.6281   True
OJ_0.5 VC_0.5    -5.25 0.0243 -10.0481  -0.4519   True
OJ_0.5 VC_1.0     3.54  0.264  -1.2581   8.3381  False
OJ_0.5 VC_2.0    12.91    0.0   8.1119  17.7081   True
OJ_1.0 OJ_2.0     3.36 0.3187  -1.4381   8.1581  False
OJ_1.0 VC_0.5   -14.72    0.0 -19.5181  -9.9219   True
OJ_1.0 VC_1.0    -5.93 0.0074 -10.7281  -1.1319   True
OJ_1.0 VC_2.0     3.44 0.2936  -1.3581   8.2381  False
OJ_2.0 VC_0.5   -18.08    0.0 -22.8781 -13.2819   True
OJ_2.0 VC_1.0    -9.29    0.0 -14.0881  -4.4919   True
OJ_2.0 VC_2.0     0.08    1.0  -4.7181   4.8781  False
VC_0.5 VC_1.0     8.79    0.0   3.9919  13.5881   True
VC_0.5 VC_2.0    18.16    0.0  13.3619  22.9581   True
VC_1.0 VC_2.0     9.37    0.0   4.5719  14.1681   True
------------------------------------------------------

# 여섯 칸 일원배치의 MSW 가 이원배치 잔차와 같은 수인가.
cells = [v.values for _, v in df.groupby('supp_dose')['len']]
SSW_cells = sum(((x - x.mean()) ** 2).sum() for x in cells)
MSE = SSW_cells / 54
print(f"6칸 일원배치 SSW = {SSW_cells:.3f} (df 54),  MSE = {MSE:.6f}")
print(f"이원배치 잔차     = 712.106   (df 54),  MS  = {712.106 / 54:.6f}")

q6 = stats.studentized_range.ppf(0.95, 6, 54)
HSD = q6 * np.sqrt(MSE / 10)
print(f"\nq(0.05, 6, 54) = {q6:.4f}   HSD = {HSD:.4f}  (15개 비교 공통)")

cell = df.groupby(['supp', 'dose'])['len'].mean()
print(f"\n{'dose':>6}{'OJ':>8}{'VC':>8}{'OJ-VC':>9}{'p-adj':>9}")
for d in (0.5, 1.0, 2.0):
    diff = cell[('OJ', d)] - cell[('VC', d)]
    p = stats.studentized_range.sf(abs(diff) / np.sqrt(MSE / 10), 6, 54)
    print(f"{d:>6}{cell[('OJ', d)]:>8.2f}{cell[('VC', d)]:>8.2f}{diff:>9.2f}{p:>9.4f}")

contrast = ((cell[('OJ', 0.5)] - cell[('VC', 0.5)])
            - (cell[('OJ', 2.0)] - cell[('VC', 2.0)]))
print(f"\n교호작용 대비 (0.5 의 격차) - (2.0 의 격차) = {contrast:.2f}")
print(f"이 대비의 표준오차 = {np.sqrt(4 * MSE / 10):.4f},  "
      f"t = {contrast / np.sqrt(4 * MSE / 10):.4f}")

6칸 일원배치 SSW = 712.106 (df 54),  MSE = 13.187148
이원배치 잔차     = 712.106   (df 54),  MS  = 13.187148

q(0.05, 6, 54) = 4.1783   HSD = 4.7981  (15개 비교 공통)

  dose      OJ      VC    OJ-VC    p-adj
   0.5   13.23    7.98     5.25   0.0243
   1.0   22.70   16.77     5.93   0.0074
   2.0   26.06   26.14    -0.08   1.0000

교호작용 대비 (0.5 의 격차) - (2.0 의 격차) = 5.33
이 대비의 표준오차 = 2.2967,  t = 2.3207

import matplotlib.pyplot as plt
from statsmodels.graphics.factorplots import interaction_plot

# 마지막으로 그림을 다시 본다. 사후비교로 갈린 쌍이 그림에서 어디인지
# 짚어 보면 결과가 몸에 붙는다.
fig, ax = plt.subplots(figsize=(8, 4))
interaction_plot(df['dose'], df['supp'], df['len'], ax=ax,
                 markers=['o', 's'], linestyles=['--', '-.'])
ax.set_title("Interaction: dose x supp")
ax.set_xlabel("dose")
ax.set_ylabel("len")
plt.tight_layout()
plt.show()

tab = df.groupby(['supp', 'dose'])['len'].agg(['mean', 'std', 'count'])
print(tab.round(3))

cell = df.groupby(['supp', 'dose'])['len'].mean()
MSE = 13.187148
print(f"\n칸평균의 표준오차 = sqrt(MSE/n) = {np.sqrt(MSE / 10):.4f}"
      f"  (그림에 그려지지 않는다)")
print(f"{'dose':>6}{'OJ-VC':>9}{'+-1.96SE_diff':>16}")
for d in (0.5, 1.0, 2.0):
    diff = cell[('OJ', d)] - cell[('VC', d)]
    sed = np.sqrt(2 * MSE / 10)
    print(f"{d:>6}{diff:>9.2f}   ({diff - 1.96 * sed:>6.2f}, {diff + 1.96 * sed:>6.2f})")

# 가로축을 '칸 번호'로 바꾸면 선의 모양이 달라 보인다.
print(f"\n선분의 기울기 (가로축이 실제 용량일 때)")
for s in ('OJ', 'VC'):
    print(f"  {s}: 0.5->1.0 {(cell[(s, 1.0)] - cell[(s, 0.5)]) / 0.5:>7.2f}"
          f"   1.0->2.0 {(cell[(s, 2.0)] - cell[(s, 1.0)]) / 1.0:>7.2f}")
print(f"선분이 올라간 양 (가로축이 순서 1,2,3 일 때)")
for s in ('OJ', 'VC'):
    print(f"  {s}: 1->2 {cell[(s, 1.0)] - cell[(s, 0.5)]:>7.2f}"
          f"   2->3 {cell[(s, 2.0)] - cell[(s, 1.0)]:>7.2f}")

            mean    std  count
supp dose
OJ   0.5   13.23  4.460     10
     1.0   22.70  3.911     10
     2.0   26.06  2.655     10
VC   0.5    7.98  2.747     10
     1.0   16.77  2.515     10
     2.0   26.14  4.798     10

칸평균의 표준오차 = sqrt(MSE/n) = 1.1484  (그림에 그려지지 않는다)
  dose    OJ-VC   +-1.96SE_diff
   0.5     5.25   (  2.07,   8.43)
   1.0     5.93   (  2.75,   9.11)
   2.0    -0.08   ( -3.26,   3.10)

선분의 기울기 (가로축이 실제 용량일 때)
  OJ: 0.5->1.0   18.94   1.0->2.0    3.36
  VC: 0.5->1.0   17.58   1.0->2.0    9.37
선분이 올라간 양 (가로축이 순서 1,2,3 일 때)
  OJ: 1->2    9.47   2->3    3.36
  VC: 1->2    8.79   2->3    9.37

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

def twoway_report(data, y, A, B, alpha=0.05):
    """이원배치 분산분석의 전 과정을 한 번에 수행한다."""
    print(f"=== {y} ~ {A} × {B} ===")

    cell = data.groupby([A, B])[y].agg(["size", "mean", "std"]).round(3)
    print("\n[1] 칸 요약")
    print(cell.to_string())
    print(f"  균형 설계: {cell['size'].nunique() == 1}")

    groups = [g[y].values for _, g in data.groupby([A, B])]
    lev = stats.levene(*groups, center="median")
    print(f"\n[2] 등분산(브라운-포사이드): "
          f"W = {lev.statistic:.4f}, p = {lev.pvalue:.4f}")

    fit = ols(f"{y} ~ C({A})*C({B})", data=data).fit()
    tab = sm.stats.anova_lm(fit, typ=2)
    SSE = tab.loc["Residual", "sum_sq"]
    MSE = SSE / tab.loc["Residual", "df"]
    tab["부분 eta2"] = tab["sum_sq"] / (tab["sum_sq"] + SSE)
    tab.loc["Residual", "부분 eta2"] = np.nan     # 잔차에는 의미가 없다
    print("\n[3] 분산분석표 (제II형)")
    print(tab.round(4).to_string())

    sw = stats.shapiro(fit.resid)
    print(f"\n[4] 잔차 정규성(샤피로): "
          f"W = {sw.statistic:.4f}, p = {sw.pvalue:.4f}")

    key = f"C({A}):C({B})"
    if tab.loc[key, "PR(>F)"] < alpha:
        print(f"\n[5] 교호작용 유의 → {A} 의 단순효과 ({B} 수준별, 본페로니)")
        lv = sorted(data[B].unique())
        m = data.groupby([A, B])[y].mean()
        a0, a1 = sorted(data[A].unique())[:2]
        n = cell["size"].min()
        dfe = int(tab.loc["Residual", "df"])
        se = np.sqrt(2 * MSE / n)
        for b in lv:
            d = m[(a1, b)] - m[(a0, b)]
            p = 2 * stats.t.sf(abs(d / se), dfe)
            print(f"  {B}={b}: {a1}-{a0} = {d:+7.3f},  "
                  f"p = {min(1, len(lv) * p):.4f}")
    else:
        print("\n[5] 교호작용 유의하지 않음 → 주효과를 해석한다")

tooth = [4.2, 11.5, 7.3, 5.8, 6.4, 10, 11.2, 11.2, 5.2, 7,
         16.5, 16.5, 15.2, 17.3, 22.5, 17.3, 13.6, 14.5, 18.8, 15.5,
         23.6, 18.5, 33.9, 25.5, 26.4, 32.5, 26.7, 21.5, 23.3, 29.5,
         15.2, 21.5, 17.6, 9.7, 14.5, 10, 8.2, 9.4, 16.5, 9.7,
         19.7, 23.3, 23.6, 26.4, 20, 25.2, 25.8, 21.2, 14.5, 27.3,
         25.5, 26.4, 22.4, 24.5, 24.8, 30.9, 26.4, 27.3, 29.4, 23]
df = pd.DataFrame({
    "len": tooth,
    "supp": ["VC"] * 30 + ["OJ"] * 30,
    "dose": ([0.5] * 10 + [1.0] * 10 + [2.0] * 10) * 2,
})
twoway_report(df, "len", "supp", "dose")

=== len ~ supp × dose ===

[1] 칸 요약
           size   mean    std
supp dose
OJ   0.5     10  13.23  4.460
     1.0     10  22.70  3.911
     2.0     10  26.06  2.655
VC   0.5     10   7.98  2.747
     1.0     10  16.77  2.515
     2.0     10  26.14  4.798
  균형 설계: True

[2] 등분산(브라운-포사이드): W = 1.7086, p = 0.1484

[3] 분산분석표 (제II형)
                    sum_sq    df       F  PR(>F)  부분 eta2
C(supp)           205.3500   1.0  15.572  0.0002   0.2238
C(dose)          2426.4343   2.0  92.000  0.0000   0.7731
C(supp):C(dose)   108.3190   2.0   4.107  0.0219   0.1320
Residual          712.1060  54.0     NaN     NaN      NaN

[4] 잔차 정규성(샤피로): W = 0.9850, p = 0.6694

[5] 교호작용 유의 → supp 의 단순효과 (dose 수준별, 본페로니)
  dose=0.5: VC-OJ =  -5.250,  p = 0.0063
  dose=1.0: VC-OJ =  -5.930,  p = 0.0018
  dose=2.0: VC-OJ =  +0.080,  p = 1.0000

[1] 칸 요약 ──→ 균형인가? 칸이 빈 곳은 없는가?
     ↓
[2] 등분산 ──→ 깨지면 웰치 이원배치나 로그 변환으로
     ↓
[3] 분산분석 ──→ 교호작용의 p 와 효과크기
     ↓
[4] 잔차 정규성 ──→ 깨지면 순열검정·변환
     ↓
[5] 분기 ──→ 유의: 단순효과 / 비유의: 주효과

import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.libqsturng import qsturng

tooth = [4.2, 11.5, 7.3, 5.8, 6.4, 10, 11.2, 11.2, 5.2, 7,
         16.5, 16.5, 15.2, 17.3, 22.5, 17.3, 13.6, 14.5, 18.8, 15.5,
         23.6, 18.5, 33.9, 25.5, 26.4, 32.5, 26.7, 21.5, 23.3, 29.5,
         15.2, 21.5, 17.6, 9.7, 14.5, 10, 8.2, 9.4, 16.5, 9.7,
         19.7, 23.3, 23.6, 26.4, 20, 25.2, 25.8, 21.2, 14.5, 27.3,
         25.5, 26.4, 22.4, 24.5, 24.8, 30.9, 26.4, 27.3, 29.4, 23]
df = pd.DataFrame({
    "len": tooth,
    "supp": ["VC"] * 30 + ["OJ"] * 30,
    "dose": ([0.5] * 10 + [1.0] * 10 + [2.0] * 10) * 2,
})

tab = sm.stats.anova_lm(ols("len ~ C(supp)*C(dose)", data=df).fit(), typ=2)
MSE = tab.loc["Residual", "sum_sq"] / tab.loc["Residual", "df"]
dfe = int(tab.loc["Residual", "df"])
se = np.sqrt(2 * MSE / 10)

q = qsturng(0.95, 6, dfe)                                # 6칸 Tukey
t_bon = stats.t.ppf(1 - 0.05 / (2 * 3), dfe)             # 관심 3개 본페로니
t_sid = stats.t.ppf(1 - (1 - 0.95**(1 / 3)) / 2, dfe)    # 관심 3개 시닥
t_raw = stats.t.ppf(0.975, dfe)                          # 보정 없음

print(f"공통 SE(차이) = {se:.4f},  df = {dfe}")
print(f"\n{'방법':>18s} {'배수':>8s} {'임계 차이':>10s}")
print(f"{'6칸 Tukey (15개)':>18s} {q / np.sqrt(2):8.4f} "
      f"{q / np.sqrt(2) * se:10.4f}")
print(f"{'본페로니 (3개)':>18s} {t_bon:8.4f} {t_bon * se:10.4f}")
print(f"{'시닥 (3개)':>18s} {t_sid:8.4f} {t_sid * se:10.4f}")
print(f"{'보정 없음':>18s} {t_raw:8.4f} {t_raw * se:10.4f}")

print("\n같은 용량에서의 OJ-VC 차이")
for d in [0.5, 1.0, 2.0]:
    diff = (df[(df.supp == "OJ") & (df.dose == d)].len.mean()
            - df[(df.supp == "VC") & (df.dose == d)].len.mean())
    print(f"  dose {d}: {diff:+7.3f}   "
          f"Tukey15 {'통과' if abs(diff) > q / np.sqrt(2) * se else ' - '}   "
          f"본페로니3 {'통과' if abs(diff) > t_bon * se else ' - '}")

공통 SE(차이) = 1.6240,  df = 54

                방법       배수      임계 차이
    6칸 Tukey (15개)   2.9545     4.7981
         본페로니 (3개)   2.4708     4.0127
           시닥 (3개)   2.4641     4.0017
             보정 없음   2.0049     3.2560

같은 용량에서의 OJ-VC 차이
  dose 0.5:  +5.250   Tukey15 통과   본페로니3 통과
  dose 1.0:  +5.930   Tukey15 통과   본페로니3 통과
  dose 2.0:  -0.080   Tukey15  -    본페로니3  -

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

tooth = [4.2, 11.5, 7.3, 5.8, 6.4, 10, 11.2, 11.2, 5.2, 7,
         16.5, 16.5, 15.2, 17.3, 22.5, 17.3, 13.6, 14.5, 18.8, 15.5,
         23.6, 18.5, 33.9, 25.5, 26.4, 32.5, 26.7, 21.5, 23.3, 29.5,
         15.2, 21.5, 17.6, 9.7, 14.5, 10, 8.2, 9.4, 16.5, 9.7,
         19.7, 23.3, 23.6, 26.4, 20, 25.2, 25.8, 21.2, 14.5, 27.3,
         25.5, 26.4, 22.4, 24.5, 24.8, 30.9, 26.4, 27.3, 29.4, 23]
df = pd.DataFrame({
    "len": tooth,
    "supp": ["VC"] * 30 + ["OJ"] * 30,
    "dose": ([0.5] * 10 + [1.0] * 10 + [2.0] * 10) * 2,
})

full = sm.stats.anova_lm(ols("len ~ C(supp)*C(dose)", data=df).fit(), typ=2)
one = sm.stats.anova_lm(ols("len ~ C(supp)", data=df).fit(), typ=2)

d = df[df.supp == "OJ"].len.mean() - df[df.supp == "VC"].len.mean()
print(f"OJ - VC = {d:.4f}  (두 검정 모두 같다)")

for lab, tab in [("supp 만 있는 모형(주변 Tukey)", one),
                 ("교호작용 포함 모형(분산분석)", full)]:
    ms = tab.loc["Residual", "sum_sq"] / tab.loc["Residual", "df"]
    dd = int(tab.loc["Residual", "df"])
    se = np.sqrt(2 * ms / 30)
    t = d / se
    print(f"  {lab}\n    MSE = {ms:8.4f}  df = {dd:3d}  "
          f"SE = {se:6.4f}  t = {t:7.4f}  p = {2 * stats.t.sf(abs(t), dd):.6f}")

sd = full.loc["C(dose)", "sum_sq"]
sr = one.loc["Residual", "sum_sq"]
print(f"\n용량이 설명하는 제곱합 = {sd:.2f}")
print(f"주변 모형의 잔차 제곱합 = {sr:.2f}  →  {100 * sd / sr:.1f}% 가 용량 탓")

OJ - VC = 3.7000  (두 검정 모두 같다)
  supp 만 있는 모형(주변 Tukey)
    MSE =  55.9803  df =  58  SE = 1.9318  t =  1.9153  p = 0.060393
  교호작용 포함 모형(분산분석)
    MSE =  13.1871  df =  54  SE = 0.9376  t =  3.9461  p = 0.000231

용량이 설명하는 제곱합 = 2426.43
주변 모형의 잔차 제곱합 = 3246.86  →  74.7% 가 용량 탓

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols
from patsy.contrasts import Sum

tooth = [4.2, 11.5, 7.3, 5.8, 6.4, 10, 11.2, 11.2, 5.2, 7,
         16.5, 16.5, 15.2, 17.3, 22.5, 17.3, 13.6, 14.5, 18.8, 15.5,
         23.6, 18.5, 33.9, 25.5, 26.4, 32.5, 26.7, 21.5, 23.3, 29.5,
         15.2, 21.5, 17.6, 9.7, 14.5, 10, 8.2, 9.4, 16.5, 9.7,
         19.7, 23.3, 23.6, 26.4, 20, 25.2, 25.8, 21.2, 14.5, 27.3,
         25.5, 26.4, 22.4, 24.5, 24.8, 30.9, 26.4, 27.3, 29.4, 23]
df = pd.DataFrame({
    "len": tooth,
    "supp": ["VC"] * 30 + ["OJ"] * 30,
    "dose": ([0.5] * 10 + [1.0] * 10 + [2.0] * 10) * 2,
})

# 두 칸에서만 관측을 떨어뜨려 불균형을 만든다
rng = np.random.default_rng(99)
drop = []
drop += list(rng.choice(df[(df.supp == "VC") & (df.dose == 0.5)].index,
                        5, replace=False))
drop += list(rng.choice(df[(df.supp == "OJ") & (df.dose == 2.0)].index,
                        3, replace=False))
un = df.drop(index=drop).reset_index(drop=True)
print("칸별 n")
print(un.groupby(["supp", "dose"]).size().unstack().to_string())

print("\n제I형 (supp 를 먼저 넣음)")
print(sm.stats.anova_lm(
    ols("len ~ C(supp)+C(dose)+C(supp):C(dose)", data=un).fit(),
    typ=1).round(4).to_string())

print("\n제I형 (dose 를 먼저 넣음)")
print(sm.stats.anova_lm(
    ols("len ~ C(dose)+C(supp)+C(dose):C(supp)", data=un).fit(),
    typ=1).round(4).to_string())

print("\n제II형")
print(sm.stats.anova_lm(ols("len ~ C(supp)*C(dose)", data=un).fit(),
                        typ=2).round(4).to_string())

print("\n제III형 (합 대비를 반드시 쓴다)")
print(sm.stats.anova_lm(
    ols("len ~ C(supp, Sum)*C(dose, Sum)", data=un).fit(),
    typ=3).round(4).to_string())

칸별 n
dose  0.5  1.0  2.0
supp
OJ     10   10    7
VC      5   10   10

제I형 (supp 를 먼저 넣음)
                   df     sum_sq   mean_sq        F  PR(>F)
C(supp)           1.0    27.6193   27.6193   2.0119  0.1628
C(dose)           2.0  1938.4311  969.2156  70.6013  0.0000
C(supp):C(dose)   2.0   100.2390   50.1195   3.6509  0.0338
Residual         46.0   631.4883   13.7280      NaN     NaN

제I형 (dose 를 먼저 넣음)
                   df     sum_sq   mean_sq        F  PR(>F)
C(dose)           2.0  1770.0452  885.0226  64.4684  0.0000
C(supp)           1.0   196.0052  196.0052  14.2778  0.0005
C(dose):C(supp)   2.0   100.2390   50.1195   3.6509  0.0338
Residual         46.0   631.4883   13.7280      NaN     NaN

제II형
                    sum_sq    df        F  PR(>F)
C(supp)           196.0052   1.0  14.2778  0.0005
C(dose)          1938.4311   2.0  70.6013  0.0000
C(supp):C(dose)   100.2390   2.0   3.6509  0.0338
Residual          631.4883  46.0      NaN     NaN

제III형 (합 대비를 반드시 쓴다)
                               sum_sq    df          F  PR(>F)
Intercept                  16925.8079   1.0  1232.9400  0.0000
C(supp, Sum)                 189.7146   1.0    13.8195  0.0005
C(dose, Sum)                1864.7392   2.0    67.9173  0.0000
C(supp, Sum):C(dose, Sum)    100.2390   2.0     3.6509  0.0338
Residual                     631.4883  46.0        NaN     NaN

0. 설계 확인 ── 칸별 n, 빈 칸, 균형 여부
   ↓
1. 기술통계 ── 칸 평균·표준편차, 교호작용 그림(오차막대 포함)
   ↓
2. 가정 점검 ── 등분산(레빈), 잔차 정규성, 독립성
   ↓
3. 분산분석 ── 균형이면 아무 유형, 불균형이면 II/III 선택
   ↓
4. 교호작용 판정
   ├─ 유의 ──→ 단순효과 (보정 필수)
   └─ 비유의 ──→ 주효과 사후검정 + 검정력 언급
   ↓
5. 효과크기·구간 ── 부분 η², ω², 대비의 신뢰구간
   ↓
6. 보고 ── 표 + 그림 + 코드