이원배치 분산분석: 교호작용 효과¶

import matplotlib.pyplot as plt
import pandas as pd
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
from statsmodels.graphics.factorplots import interaction_plot

def load_data():
    """ToothGrowth 자료를 읽는다. 보충제 종류와 투여량, 그리고 치아 길이다."""
    url = 'https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/1dcc2bf5f955cc1224a3e1307256e1fe86b68dae/csv/datasets/ToothGrowth.csv'
    df = pd.read_csv(url, usecols=[1, 2, 3])
    return df

def plot_interaction(df):
    """교호작용 그림. 두 선이 나란하면 교호작용이 없다는 뜻이다.

    선이 벌어지거나 엇갈리면 한 요인의 효과가 다른 요인의 수준에 따라
    달라진다는 것이고, 그때는 주효과만 말해서는 안 된다.
    """
    fig, ax = plt.subplots(figsize=(12, 3))
    interaction_plot(df.dose, df.supp, df.len,
                     colors=['red', 'blue'],
                     markers=['*', 'P'],
                     markersize=7, ax=ax,
                     legendloc='lower right',
                     linestyles=["--", "--"])
    ax.set_title("Interaction Plot: Dose vs. Supplement on Tooth Length")
    ax.set_xlabel("Dose Level")
    ax.set_ylabel("Tooth Length")
    plt.show()

def perform_two_way_anova(df):
    """이원배치 분산분석. 그림에서 본 것이 통계적으로도 유의한지 확인한다."""
    model = ols('len ~ C(supp) + C(dose) + C(supp):C(dose)', data=df).fit()
    anova_results = anova_lm(model)
    print("\nTwo-Way ANOVA Results:")
    print(anova_results)

df = load_data()
plot_interaction(df)
perform_two_way_anova(df)

Two-Way ANOVA Results:
                   df       sum_sq      mean_sq          F        PR(>F)
C(supp)           1.0   205.350000   205.350000  15.571979  2.311828e-04
C(dose)           2.0  2426.434333  1213.217167  91.999965  4.046291e-18
C(supp):C(dose)   2.0   108.319000    54.159500   4.106991  2.186027e-02
Residual         54.0   712.106000    13.187148        NaN           NaN

import numpy as np

c = 10
M = df.pivot_table(index='supp', columns='dose', values='len', aggfunc='mean')
print("칸평균")
print(M)

row, col = M.mean(axis=1), M.mean(axis=0)
grand = M.values.mean()
SSA = 3 * c * ((row - grand) ** 2).sum()
SSB = 2 * c * ((col - grand) ** 2).sum()
D = M.values - row.values[:, None] - col.values[None, :] + grand
SSAB = c * (D ** 2).sum()
print(f"\n행평균 {row.values.round(4)}  열평균 {col.values.round(4)}  전체 {grand:.4f}")
print(f"SSA  = {SSA:.4f}   (표: 205.3500)")
print(f"SSB  = {SSB:.4f}   (표: 2426.4343)")
print(f"SSAB = {SSAB:.4f}   (표: 108.3190)")

# a=2 의 지름길: 격차의 흩어짐
delta = M.loc['OJ'] - M.loc['VC']
print(f"\n격차 delta_j = {delta.values.round(2)},  평균 {delta.mean():.2f}")
print(f"(c/2) * sum (delta - mean)^2 = {c / 2 * ((delta - delta.mean()) ** 2).sum():.4f}")

SSE = sum(((v.values - v.values.mean()) ** 2).sum()
          for _, v in df.groupby(['supp', 'dose'])['len'])
SST = ((df['len'] - df['len'].mean()) ** 2).sum()
print(f"\nSSE = {SSE:.4f},  네 제곱합의 합 = {SSA + SSB + SSAB + SSE:.4f},  SST = {SST:.4f}")

칸평균
dose    0.5    1.0    2.0
supp
OJ    13.23  22.70  26.06
VC     7.98  16.77  26.14

행평균 [20.6633 16.9633]  열평균 [10.605 19.735 26.1  ]  전체 18.8133
SSA  = 205.3500   (표: 205.3500)
SSB  = 2426.4343   (표: 2426.4343)
SSAB = 108.3190   (표: 108.3190)

격차 delta_j = [ 5.25  5.93 -0.08],  평균 3.70
(c/2) * sum (delta - mean)^2 = 108.3190

SSE = 712.1060,  네 제곱합의 합 = 3452.2093,  SST = 3452.2093

from statsmodels.stats.multicomp import pairwise_tukeyhsd

# 1단계: 이원배치 분산분석
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()
anova_results = anova_lm(model)
print(anova_results, end="\n\n")

# 2단계: 주효과에 대한 Tukey HSD
tukey_dose = pairwise_tukeyhsd(endog=df['len'], groups=df['dose'], alpha=0.05)
print(tukey_dose, end="\n\n")

tukey_supp = pairwise_tukeyhsd(endog=df['len'], groups=df['supp'], alpha=0.05)
print(tukey_supp, end="\n\n")

# Step 3: 교호작용에 대한 사후검정
# 두 요인을 붙여 하나의 요인으로 만든다. 수준이 2 x 3 = 6개가 되고
# 비교는 15쌍으로 늘어난다. 그만큼 Tukey의 보정도 커진다.
df['supp_dose'] = df['supp'].astype(str) + "_" + df['dose'].astype(str)
tukey_interaction = pairwise_tukeyhsd(endog=df['len'], groups=df['supp_dose'], alpha=0.05)
print(tukey_interaction, end="\n\n")

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

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
-------------------------------------------------

 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
------------------------------------------------------

import numpy as np
from scipy import stats

# dose 주효과 Tukey 표를 스튜던트화 범위 분포로 직접 만든다.
g = df.groupby('dose')['len']
m_i = g.mean().values
n_i = g.count().values
k, N = len(n_i), n_i.sum()
dfe = N - k
MSE = ((n_i - 1) * g.var(ddof=1).values).sum() / dfe
se = np.sqrt(MSE / n_i[0])                      # 집단평균 하나의 표준오차
q_crit = stats.studentized_range.ppf(0.95, k, dfe)

print(f"MSE = {MSE:.6f},  오차 자유도 = {dfe}")
print(f"집단평균 하나의 표준오차 sqrt(MSE/m) = {se:.6f}")
print(f"q_crit(0.05, k={k}, df={dfe}) = {q_crit:.4f}")
print(f"HSD = q_crit * sqrt(MSE/m)   = {q_crit * se:.4f}   <- 모든 쌍에 같은 폭")
print(f"대조: 보정 없는 t 임계값 {stats.t.ppf(0.975, dfe):.4f},  "
      f"Tukey 는 t 척도로 {q_crit / np.sqrt(2):.4f}")

lv = g.mean().index.values
print(f"\n{'group1':>6s} {'group2':>6s} {'meandiff':>9s} {'q':>8s} "
      f"{'p-adj':>9s} {'lower':>9s} {'upper':>9s}")
for a in range(k):
    for b in range(a + 1, k):
        d = m_i[b] - m_i[a]
        q = abs(d) / se
        print(f"{lv[a]:>6.1f} {lv[b]:>6.1f} {d:9.4f} {q:8.4f} "
              f"{stats.studentized_range.sf(q, k, dfe):9.6f} "
              f"{d - q_crit * se:9.4f} {d + q_crit * se:9.4f}")

MSE = 17.996053,  오차 자유도 = 57
집단평균 하나의 표준오차 sqrt(MSE/m) = 0.948579
q_crit(0.05, k=3, df=57) = 3.4032
HSD = q_crit * sqrt(MSE/m)   = 3.2282   <- 모든 쌍에 같은 폭
대조: 보정 없는 t 임계값 2.0025,  Tukey 는 t 척도로 2.4064

group1 group2  meandiff        q     p-adj     lower     upper
   0.5    1.0    9.1300   9.6249  0.000000    5.9018   12.3582
   0.5    2.0   15.4950  16.3350  0.000000   12.2668   18.7232
   1.0    2.0    6.3650   6.7100  0.000042    3.1368    9.5932

import numpy as np
from itertools import combinations
from scipy import stats

# 평균이 서로 같은 집단들 사이에서 "한 번이라도 거짓 기각이 나는" 확률을 센다.
# 이것이 집단별오류율(FWER)이고 이론값은 0.05 다.
rng_sim = np.random.default_rng(2026)
n, B = 10, 2000
print(f"집단당 n = {n}, 반복 {B}회, 명목 FWER 0.05, "
      f"몬테카를로 오차 {np.sqrt(0.05 * 0.95 / B):.4f}")
print(f"\n{'k':>3s} {'delta':>6s} {'같은쌍':>7s} {'F 유의':>7s} "
      f"{'Tukey':>7s} {'보호 LSD':>9s}")
for k in [3, 4, 5, 6]:
    dfe = k * n - k
    q_crit = stats.studentized_range.ppf(0.95, k, dfe)
    t_crit = stats.t.ppf(0.975, dfe)
    f_crit = stats.f.ppf(0.95, k - 1, dfe)
    for delta in [0.0, 8.0]:
        # delta = 0 이면 완전 영가설, delta = 8 이면 한 집단만 멀리 떨어진 부분 영가설
        mu = np.array([0.0] * (k - 1) + [delta])
        pairs = (list(combinations(range(k), 2)) if delta == 0.0
                 else list(combinations(range(k - 1), 2)))
        Y = rng_sim.normal(mu[None, :, None], 1.0, size=(B, k, n))
        gm = Y.mean(axis=2)
        MSW = ((Y - gm[:, :, None]) ** 2).sum(axis=(1, 2)) / dfe
        SSB = (n * (gm - gm.mean(axis=1, keepdims=True)) ** 2).sum(axis=1)
        sig_F = (SSB / (k - 1)) / MSW > f_crit

        any_tukey = np.zeros(B, bool)
        any_lsd = np.zeros(B, bool)
        for a, b in pairs:
            d = np.abs(gm[:, a] - gm[:, b])
            any_tukey |= d / np.sqrt(MSW / n) > q_crit
            any_lsd |= d / np.sqrt(MSW * 2 / n) > t_crit
        print(f"{k:3d} {delta:6g} {len(pairs):7d} {sig_F.mean():7.4f} "
              f"{any_tukey.mean():7.4f} {(sig_F & any_lsd).mean():9.4f}")

집단당 n = 10, 반복 2000회, 명목 FWER 0.05, 몬테카를로 오차 0.0049

  k  delta     같은쌍    F 유의   Tukey    보호 LSD
  3      0       3  0.0595  0.0560    0.0595
  3      8       1  1.0000  0.0195    0.0520
  4      0       6  0.0480  0.0425    0.0480
  4      8       3  1.0000  0.0320    0.1185
  5      0      10  0.0500  0.0495    0.0500
  5      8       6  1.0000  0.0355    0.1980
  6      0      15  0.0525  0.0545    0.0525
  6      8      10  1.0000  0.0330    0.2810

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

# 교수법 둘과 학습시간 둘, 칸마다 두 명씩인 2x2 설계다.
# 칸당 반복이 둘뿐이라 자유도가 매우 적다. 교호작용을 재려면 칸마다
# 반복이 적어도 둘은 있어야 한다는 요구를 겨우 맞춘 셈이다.
data = {
    'Teaching_Method': ['Traditional', 'Traditional', 'Traditional', 'Traditional',
                        'Online', 'Online', 'Online', 'Online'],
    'Study_Time': ['1 Hour', '1 Hour', '2 Hours', '2 Hours',
                   '1 Hour', '1 Hour', '2 Hours', '2 Hours'],
    'Score': [60, 62, 68, 70, 65, 63, 72, 74]
}
df = pd.DataFrame(data)

model = ols('Score ~ C(Teaching_Method) + C(Study_Time) + C(Teaching_Method):C(Study_Time)', data=df).fit()
anova_results = anova_lm(model)
print("Two-Way ANOVA Results:")
print(anova_results)

Two-Way ANOVA Results:
                                   df  sum_sq  mean_sq      F    PR(>F)
C(Teaching_Method)                1.0    24.5     24.5  12.25  0.024896
C(Study_Time)                     1.0   144.5    144.5  72.25  0.001051
C(Teaching_Method):C(Study_Time)  1.0     0.5      0.5   0.25  0.643330
Residual                          4.0     8.0      2.0    NaN       NaN

import warnings

import numpy as np
from scipy import stats

c = 2
M = df.pivot_table(index='Teaching_Method', columns='Study_Time',
                   values='Score', aggfunc='mean')
print("칸평균")
print(M)

row, col = M.mean(axis=1), M.mean(axis=0)
grand = M.values.mean()
SSA = 2 * c * ((row - grand) ** 2).sum()
SSB = 2 * c * ((col - grand) ** 2).sum()
L = (M.iloc[0, 0] - M.iloc[0, 1]) - (M.iloc[1, 0] - M.iloc[1, 1])
SSAB = c * L ** 2 / 4
SSE = sum(((v.values - v.values.mean()) ** 2).sum()
          for _, v in df.groupby(['Teaching_Method', 'Study_Time'])['Score'])
print(f"\n행평균 {row.values}  열평균 {col.values}  전체 {grand}")
print(f"SSA = {SSA}   SSB = {SSB}   L = {L}   SSAB = c L^2/4 = {SSAB}")
print(f"SSE = {SSE}   합 = {SSA + SSB + SSAB + SSE}   "
      f"SST = {((df['Score'] - df['Score'].mean()) ** 2).sum()}")

MSE = SSE / 4
se_L = np.sqrt(4 * MSE / c)
t_L = L / se_L
print(f"\nMSE = {MSE},  SE(L) = {se_L},  t = {t_L},  t^2 = {t_L ** 2}")

# 자유도 4 에서 잡아낼 수 있는 교호작용의 크기는 얼마인가.
t_crit = stats.t.ppf(0.975, 4)
print(f"\nt(0.975, 4) = {t_crit:.4f}  ->  |L| > {t_crit * se_L:.4f} 이어야 유의하다")
print(f"견줄 것: 학습시간 주효과의 크기 = {col.iloc[1] - col.iloc[0]:.1f}")

F_crit = stats.f(1, 4).ppf(0.95)
print(f"\n{'참 L':>8}{'lambda':>10}{'검정력':>10}")
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    for L_true in (1, 2, 4, 6, 8):
        lam = c * L_true ** 2 / (4 * MSE)
        print(f"{L_true:>8}{lam:>10.3f}{stats.ncf(1, 4, lam).sf(F_crit):>10.4f}")

칸평균
Study_Time       1 Hour  2 Hours
Teaching_Method
Online             64.0     73.0
Traditional        61.0     69.0

행평균 [68.5 65. ]  열평균 [62.5 71. ]  전체 66.75
SSA = 24.5   SSB = 144.5   L = -1.0   SSAB = c L^2/4 = 0.5
SSE = 8.0   합 = 177.5   SST = 177.5

MSE = 2.0,  SE(L) = 2.0,  t = -0.5,  t^2 = 0.25

t(0.975, 4) = 2.7764  ->  |L| > 5.5529 이어야 유의하다
견줄 것: 학습시간 주효과의 크기 = 8.5

     참 L    lambda       검정력
       1     0.250    0.0679
       2     1.000    0.1222
       4     4.000    0.3360
       6     9.000    0.6193
       8    16.000    0.8433

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

import numpy as np
from scipy import stats

def n_needed(delta, var_coef, sigma=1.0, power=0.80, alpha=0.05):
    """대비의 분산이 var_coef·σ²/n 일 때 크기 delta 를 잡는 데 필요한 n."""
    for n in range(3, 100_000):
        se = sigma * np.sqrt(var_coef / n)
        df = 4 * (n - 1)
        ncp = (delta / se)**2
        if stats.ncf.sf(stats.f.ppf(1 - alpha, 1, df), 1, df, ncp) >= power:
            return n

print(f"{'δ':>5s} {'주효과 n':>10s} {'교호작용 n':>12s} {'비':>6s}")
for d in [0.3, 0.5, 0.8, 1.0]:
    a = n_needed(d, 1.0)
    b = n_needed(d, 4.0)
    print(f"{d:5.1f} {a:10d} {b:12d} {b / a:6.2f}")

    δ      주효과 n       교호작용 n      비
  0.3         88          350   3.98
  0.5         32          127   3.97
  0.8         13           50   3.85
  1.0          9           32   3.56

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

rng = np.random.default_rng(1122)
M, sigma = 2_000, 1.0
print(f"\n2×2, 칸당 n, 칸 평균 (0, d, d, 3d) — 교호작용 크기 = d")
print(f"{'n':>4s} {'d':>5s} {'A 주효과':>10s} {'B 주효과':>10s} {'교호작용':>10s}")
for n in [10, 20, 40]:
    for d in [0.5, 1.0]:
        mu = {(0, 0): 0.0, (0, 1): d, (1, 0): d, (1, 1): 3 * d}
        a = b = c = 0
        for _ in range(M):
            frames = []
            for i in range(2):
                for j in range(2):
                    frames.append(pd.DataFrame(
                        {"y": rng.normal(mu[(i, j)], sigma, n),
                         "A": i, "B": j}))
            df = pd.concat(frames)
            tab = sm.stats.anova_lm(ols("y ~ C(A)*C(B)", data=df).fit(), typ=2)
            a += tab.loc["C(A)", "PR(>F)"] < 0.05
            b += tab.loc["C(B)", "PR(>F)"] < 0.05
            c += tab.loc["C(A):C(B)", "PR(>F)"] < 0.05
        print(f"{n:4d} {d:5.1f} {a / M:10.4f} {b / M:10.4f} {c / M:10.4f}")

2×2, 칸당 n, 칸 평균 (0, d, d, 3d) — 교호작용 크기 = d
   n     d      A 주효과      B 주효과       교호작용
  10   0.5     0.6425     0.6515     0.1105
  10   1.0     0.9965     0.9945     0.3260
  20   0.5     0.9070     0.9045     0.1995
  20   1.0     1.0000     1.0000     0.6050
  40   0.5     0.9960     0.9980     0.3225
  40   1.0     1.0000     1.0000     0.8725

import numpy as np
from scipy import stats

SSA, SSB, SSAB, SSE = 30.0, 80.0, 24.0, 60.0
a, b, c = 2, 3, 4                      # 2×3 요인, 칸당 반복 4
N = a * b * c
df_a, df_b, df_ab, df_e = a - 1, b - 1, (a - 1) * (b - 1), a * b * (c - 1)
SS_total = SSA + SSB + SSAB + SSE

print(f"{'원천':>8s} {'SS':>8s} {'df':>4s} {'MS':>9s} {'F':>9s} {'p':>9s}")
MSE = SSE / df_e
for name, ss, dfx in [("A", SSA, df_a), ("B", SSB, df_b),
                      ("A×B", SSAB, df_ab)]:
    ms = ss / dfx
    F = ms / MSE
    print(f"{name:>8s} {ss:8.2f} {dfx:4d} {ms:9.4f} {F:9.4f} "
          f"{stats.f.sf(F, dfx, df_e):9.4f}")
print(f"{'오차':>8s} {SSE:8.2f} {df_e:4d} {MSE:9.4f}")
print(f"{'전체':>8s} {SS_total:8.2f} {N - 1:4d}")

print(f"\n효과크기")
for name, ss, dfx in [("A", SSA, df_a), ("B", SSB, df_b),
                      ("A×B", SSAB, df_ab)]:
    eta2 = ss / SS_total
    partial = ss / (ss + SSE)
    omega2 = (ss - dfx * MSE) / (SS_total + MSE)
    print(f"  {name:>4s}:  η² = {eta2:.4f},  부분 η² = {partial:.4f},  "
          f"ω² = {omega2:.4f}")

      원천       SS   df        MS         F         p
       A    30.00    1   30.0000    9.0000    0.0077
       B    80.00    2   40.0000   12.0000    0.0005
     A×B    24.00    2   12.0000    3.6000    0.0484
      오차    60.00   18    3.3333
      전체   194.00   23

효과크기
     A:  η² = 0.1546,  부분 η² = 0.3333,  ω² = 0.1351
     B:  η² = 0.4124,  부분 η² = 0.5714,  ω² = 0.3716
   A×B:  η² = 0.1237,  부분 η² = 0.2857,  ω² = 0.0878

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

rng = np.random.default_rng(7788)
a, b = 4, 5
base = np.add.outer(np.array([0., 1, 2, 3]), np.array([0., 0.5, 1, 1.5, 2]))

def tukey_nonadditivity(Y):
    """반복 없는 이원배치에서 터키의 1자유도 비가법성 검정."""
    a, b = Y.shape
    grand = Y.mean()
    alpha = Y.mean(1) - grand
    beta = Y.mean(0) - grand
    fitted = grand + alpha[:, None] + beta[None, :]
    resid = Y - fitted
    SSE = (resid**2).sum()
    SS_na = ((resid * np.outer(alpha, beta)).sum()**2
             / ((alpha**2).sum() * (beta**2).sum()))
    df_e = (a - 1) * (b - 1)
    F = SS_na / ((SSE - SS_na) / (df_e - 1))
    return SS_na, SSE, F, stats.f.sf(F, 1, df_e - 1)

# ① 진짜 가법 모형
Y1 = base + rng.normal(0, 0.8, (a, b))
df1 = pd.DataFrame([{"y": Y1[i, j], "A": i, "B": j}
                    for i in range(a) for j in range(b)])
print("① 가법 자료 — 반복 없는 이원배치")
print(sm.stats.anova_lm(ols("y ~ C(A)+C(B)", data=df1).fit(), typ=2).round(4))
ss, sse, F, p = tukey_nonadditivity(Y1)
print(f"  터키 비가법성: SS = {ss:.4f} / SSE = {sse:.4f},  "
      f"F(1,{(a - 1) * (b - 1) - 1}) = {F:.4f},  p = {p:.4f}")

# ② 곱셈형(비가법) 자료
Y2 = np.exp(base * 0.6) + rng.normal(0, 0.8, (a, b))
ss2, sse2, F2, p2 = tukey_nonadditivity(Y2)
print(f"\n② 곱셈형 자료")
print(f"  터키 비가법성: SS = {ss2:.4f} / SSE = {sse2:.4f},  "
      f"F = {F2:.4f},  p = {p2:.4f}")

① 가법 자료 — 반복 없는 이원배치
           sum_sq    df        F  PR(>F)
C(A)      45.9658   3.0  27.2948  0.0000
C(B)      19.6479   4.0   8.7503  0.0015
Residual   6.7362  12.0      NaN     NaN
  터키 비가법성: SS = 0.2574 / SSE = 6.7362,  F(1,11) = 0.4371,  p = 0.5221

② 곱셈형 자료
  터키 비가법성: SS = 44.3733 / SSE = 58.1807,  F = 35.3512,  p = 0.0001

Y2_log = np.log(np.clip(Y2, 0.1, None))
ss3, sse3, F3, p3 = tukey_nonadditivity(Y2_log)
print(f"로그 변환 후: SS = {ss3:.4f} / SSE = {sse3:.4f},  "
      f"F = {F3:.4f},  p = {p3:.4f}")

로그 변환 후: SS = 1.7980 / SSE = 4.5126,  F = 7.2859,  p = 0.0207

import numpy as np
import pandas as pd
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 = ols("len ~ C(supp) * C(dose)", data=df).fit()
add = ols("len ~ C(supp) + C(dose)", data=df).fit()
print("교호작용 포함 모형")
print(sm.stats.anova_lm(full, typ=2).round(4))
print("\n주효과만 있는 모형")
print(sm.stats.anova_lm(add, typ=2).round(4))
print(f"\nR²:  포함 {full.rsquared:.4f}   주효과만 {add.rsquared:.4f}")
print(f"조정 R²:  포함 {full.rsquared_adj:.4f}   주효과만 {add.rsquared_adj:.4f}")
print(f"AIC:  포함 {full.aic:.4f}   주효과만 {add.aic:.4f}")

교호작용 포함 모형
                    sum_sq    df       F  PR(>F)
C(supp)           205.3500   1.0  15.572  0.0002
C(dose)          2426.4343   2.0  92.000  0.0000
C(supp):C(dose)   108.3190   2.0   4.107  0.0219
Residual          712.1060  54.0     NaN     NaN

주효과만 있는 모형
             sum_sq    df        F  PR(>F)
C(supp)    205.3500   1.0  14.0166  0.0004
C(dose)   2426.4343   2.0  82.8109  0.0000
Residual   820.4250  56.0      NaN     NaN

R²:  포함 0.7937   주효과만 0.7623
조정 R²:  포함 0.7746   주효과만 0.7496
AIC:  포함 330.7056   주효과만 335.2013

print("주효과만 있는 모형의 예측값 (칸별)")
grid = pd.DataFrame([{"supp": s, "dose": d}
                     for s in ["OJ", "VC"] for d in [0.5, 1.0, 2.0]])
grid["예측"] = add.predict(grid).round(4)
grid["실제"] = [df[(df.supp == r.supp) & (df.dose == r.dose)].len.mean()
               for r in grid.itertuples()]
grid["잔차"] = (grid["실제"] - grid["예측"]).round(4)
print(grid.to_string(index=False))

주효과만 있는 모형의 예측값 (칸별)
supp  dose     예측    실제     잔차
  OJ   0.5 12.455 13.23  0.775
  OJ   1.0 21.585 22.70  1.115
  OJ   2.0 27.950 26.06 -1.890
  VC   0.5  8.755  7.98 -0.775
  VC   1.0 17.885 16.77 -1.115
  VC   2.0 24.250 26.14  1.890

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

rng = np.random.default_rng(4321)
n, sigma = 20, 1.0

def build(cells, rng, n=20, sigma=1.0):
    frames = []
    for (i, j), mu in cells.items():
        frames.append(pd.DataFrame({"y": rng.normal(mu, sigma, n),
                                    "A": i, "B": j}))
    return pd.concat(frames, ignore_index=True)

scenarios = {
    "교호작용 없음": {(0, 0): 0.0, (0, 1): 1.0, (1, 0): 2.0, (1, 1): 3.0},
    "순서형 (크기만 다름)": {(0, 0): 0.0, (0, 1): 1.0, (1, 0): 2.0, (1, 1): 4.0},
    "비순서형 (교차)": {(0, 0): 0.0, (0, 1): 2.0, (1, 0): 2.0, (1, 1): 0.0},
}
for label, cells in scenarios.items():
    df = build(cells, rng, n, sigma)
    tab = sm.stats.anova_lm(ols("y ~ C(A)*C(B)", data=df).fit(), typ=2)
    m = df.groupby(["A", "B"]).y.mean()
    simple0 = m[(1, 0)] - m[(0, 0)]      # B=0 에서 A 의 효과
    simple1 = m[(1, 1)] - m[(0, 1)]      # B=1 에서 A 의 효과
    marg = df.groupby("A").y.mean()
    print(f"[{label}]")
    print(f"  A 의 단순효과:  B=0 에서 {simple0:+.3f},  B=1 에서 {simple1:+.3f}")
    print(f"  A 의 주효과(주변):  {marg[1] - marg[0]:+.3f}")
    print(f"  A: p={tab.loc['C(A)', 'PR(>F)']:.4f}   "
          f"B: p={tab.loc['C(B)', 'PR(>F)']:.4f}   "
          f"A×B: p={tab.loc['C(A):C(B)', 'PR(>F)']:.4f}\n")

[교호작용 없음]
  A 의 단순효과:  B=0 에서 +2.143,  B=1 에서 +1.848
  A 의 주효과(주변):  +1.996
  A: p=0.0000   B: p=0.0000   A×B: p=0.5066

[순서형 (크기만 다름)]
  A 의 단순효과:  B=0 에서 +1.800,  B=1 에서 +3.089
  A 의 주효과(주변):  +2.444
  A: p=0.0000   B: p=0.0000   A×B: p=0.0018

[비순서형 (교차)]
  A 의 단순효과:  B=0 에서 +1.505,  B=1 에서 -1.700
  A 의 주효과(주변):  -0.097
  A: p=0.6632   B: p=0.5486   A×B: p=0.0000

교호작용 F 검정
    │
    ├─ 유의하지 않음
    │     ├─ 검정력이 충분했는가? (효과크기·구간 확인)
    │     └─ 주효과를 해석하되 "교호작용을 찾지 못했다"고 서술
    │
    └─ 유의함
          ├─ 단순효과를 계산한다
          │     ├─ 부호가 같음(순서형)
          │     │     └─→ 주효과 + "크기가 다름" 명시
          │     └─ 부호가 다름(비순서형)
          │           └─→ 주효과 보고 금지, 단순효과만
          ├─ 다중비교 보정
          └─ 순서형 요인이면 교호작용을 선형·이차로 분해