교호작용 효과 그림¶

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

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

# 가로축이 한 요인, 선의 색이 다른 요인, 세로축이 반응의 평균이다.
# 두 선이 평행이면 교호작용이 없고, 벌어지거나 엇갈리면 있는 것이다.
fig, ax = plt.subplots(figsize=(10, 4))
interaction_plot(df['dose'], df['supp'], df['len'],
                 ax=ax, markers=['o', 's'], linestyles=['--', '-.'])
ax.set_title("Interaction: dose x supp on tooth length")
ax.set_xlabel("dose")
ax.set_ylabel("len")
plt.tight_layout()
plt.show()

import numpy as np
from scipy import stats

MSE, c = 13.187148, 10          # 이원배치 잔차 평균제곱과 칸 크기
M = df.pivot_table(index='supp', columns='dose', values='len', aggfunc='mean')
print("칸평균")
print(M)

se_gap = np.sqrt(2 * MSE / c)
t_crit = stats.t.ppf(0.975, 54)
print(f"\n격차의 표준오차 = sqrt(2*MSE/n) = {se_gap:.4f},  t(0.975,54) = {t_crit:.4f}")
print(f"{'dose':>6}{'OJ':>8}{'VC':>8}{'OJ-VC':>9}{'95% CI':>22}")
for d in (0.5, 1.0, 2.0):
    g = M.loc['OJ', d] - M.loc['VC', d]
    lo, hi = g - t_crit * se_gap, g + t_crit * se_gap
    print(f"{d:>6}{M.loc['OJ', d]:>8.2f}{M.loc['VC', d]:>8.2f}{g:>9.2f}"
          f"   ({lo:>6.2f}, {hi:>6.2f})")

# 표본에서 격차의 부호가 바뀌는가.
gaps = (M.loc['OJ'] - M.loc['VC']).values
print(f"\n격차의 부호: {np.sign(gaps)}  ->  표본 평균선은 실제로 교차한다")
print(f"다만 마지막 격차 -0.08 은 표준오차 {se_gap:.2f} 의 {abs(gaps[-1]) / se_gap:.3f} 배다")

for s in ('OJ', 'VC'):
    print(f"{s}: 0.5->1.0 {M.loc[s, 1.0] - M.loc[s, 0.5]:+.2f}, "
          f"1.0->2.0 {M.loc[s, 2.0] - M.loc[s, 1.0]:+.2f}")

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

격차의 표준오차 = sqrt(2*MSE/n) = 1.6240,  t(0.975,54) = 2.0049
  dose      OJ      VC    OJ-VC                95% CI
   0.5   13.23    7.98     5.25   (  1.99,   8.51)
   1.0   22.70   16.77     5.93   (  2.67,   9.19)
   2.0   26.06   26.14    -0.08   ( -3.34,   3.18)

격차의 부호: [ 1.  1. -1.]  ->  표본 평균선은 실제로 교차한다
다만 마지막 격차 -0.08 은 표준오차 1.62 의 0.049 배다
OJ: 0.5->1.0 +9.47, 1.0->2.0 +3.36
VC: 0.5->1.0 +8.79, 1.0->2.0 +9.37

import numpy as np

rng = np.random.default_rng(2468)
sigma = 1.0

print(f"{'n':>4s} {'SE(ψ̂)':>9s} {'평균|ψ̂|':>9s} {'|ψ̂|>1 비율':>11s}")
for n in [5, 10, 20, 40]:
    M = 4_000
    se = sigma * np.sqrt(4 / n)
    psi = []
    for _ in range(M):
        m = rng.normal(0, sigma / np.sqrt(n), 4)   # 참 교호작용은 0
        psi.append(abs((m[0] - m[1]) - (m[2] - m[3])))
    psi = np.array(psi)
    print(f"{n:4d} {se:9.4f} {psi.mean():9.4f} {(psi > 1).mean():11.4f}")

   n    SE(ψ̂)    평균|ψ̂|   |ψ̂|>1 비율
   5    0.8944    0.7197      0.2660
  10    0.6325    0.5019      0.1138
  20    0.4472    0.3541      0.0248
  40    0.3162    0.2543      0.0020

import numpy as np
import pandas as pd

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,
})

m = df.groupby(["supp", "dose"]).len.mean().unstack()
print(m.round(3))

print("\n배치 A (가로축 dose, 선 supp) — 선의 구간별 기울기")
for s in ["OJ", "VC"]:
    r = m.loc[s]
    print(f"  {s}: 0.5→1.0 {(r[1.0] - r[0.5]) / 0.5:+7.3f},  "
          f"1.0→2.0 {(r[2.0] - r[1.0]) / 1.0:+7.3f}")

print("\n배치 B (가로축 supp, 선 dose) — VC→OJ 의 증가량")
for d in [0.5, 1.0, 2.0]:
    print(f"  dose {d}: {m.loc['OJ', d] - m.loc['VC', d]:+7.3f}")

print("\n교호작용 대비 (dose 0.5 vs 2.0) × (OJ vs VC)")
psi_B = (m.loc['OJ', 0.5] - m.loc['VC', 0.5]) - (m.loc['OJ', 2.0] - m.loc['VC', 2.0])
psi_A = (m.loc['OJ', 0.5] - m.loc['OJ', 2.0]) - (m.loc['VC', 0.5] - m.loc['VC', 2.0])
print(f"  배치 A 로 읽기: {psi_A:+.3f}")
print(f"  배치 B 로 읽기: {psi_B:+.3f}")

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

배치 A (가로축 dose, 선 supp) — 선의 구간별 기울기
  OJ: 0.5→1.0 +18.940,  1.0→2.0  +3.360
  VC: 0.5→1.0 +17.580,  1.0→2.0  +9.370

배치 B (가로축 supp, 선 dose) — VC→OJ 의 증가량
  dose 0.5:  +5.250
  dose 1.0:  +5.930
  dose 2.0:  -0.080

교호작용 대비 (dose 0.5 vs 2.0) × (OJ vs VC)
  배치 A 로 읽기: +5.330
  배치 B 로 읽기: +5.330

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

rng = np.random.default_rng(1357)
A = np.array([0.0, 0.8])
B = np.array([0.0, 0.5, 1.0])
n = 12

rows = []
for i, a in enumerate(A):
    for j, b in enumerate(B):
        mu = np.exp(a + b)
        # 오차도 곱셈형: y = μ · exp(ε)
        rows.append(pd.DataFrame({"y": mu * np.exp(rng.normal(0, 0.25, n)),
                                  "A": i, "B": j}))
df = pd.concat(rows, ignore_index=True)

print("원 척도 칸 평균")
print(df.groupby(["A", "B"]).y.mean().unstack().round(3))
print(sm.stats.anova_lm(ols("y ~ C(A)*C(B)", data=df).fit(), typ=2).round(4))

df["ly"] = np.log(df.y)
print("\n로그 척도 칸 평균")
print(df.groupby(["A", "B"]).ly.mean().unstack().round(3))
print(sm.stats.anova_lm(ols("ly ~ C(A)*C(B)", data=df).fit(), typ=2).round(4))

m = df.groupby(["A", "B"]).y.mean().unstack()
ml = df.groupby(["A", "B"]).ly.mean().unstack()
print("\n원 척도 A 의 단순효과:", (m.loc[1] - m.loc[0]).round(3).tolist())
print("로그 척도 A 의 단순효과:", (ml.loc[1] - ml.loc[0]).round(3).tolist())

원 척도 칸 평균
B      0      1      2
A
0  1.212  1.728  2.852
1  2.033  3.396  5.337
            sum_sq    df        F  PR(>F)
C(A)       49.4801   1.0  73.6164  0.0000
C(B)       74.7443   2.0  55.6022  0.0000
C(A):C(B)   8.2963   2.0   6.1716  0.0035
Residual   44.3609  66.0      NaN     NaN

로그 척도 칸 평균
B      0      1      2
A
0  0.146  0.529  0.976
1  0.684  1.204  1.649
           sum_sq    df        F  PR(>F)
C(A)       7.1135   1.0  97.6776  0.0000
C(B)       9.6745   2.0  66.4218  0.0000
C(A):C(B)  0.0740   2.0   0.5079  0.6041
Residual   4.8066  66.0      NaN     NaN

원 척도 A 의 단순효과: [0.822, 1.668, 2.485]
로그 척도 A 의 단순효과: [0.538, 0.675, 0.673]

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,
})

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, n = int(tab.loc["Residual", "df"]), 10
se = np.sqrt(2 * MSE / n)
print(f"MSE = {MSE:.4f},  df = {dfe},  SE(차이) = {se:.4f}")

print(f"\n{'dose':>5s} {'OJ':>7s} {'VC':>7s} {'OJ-VC':>8s} "
      f"{'95% CI':>20s} {'p':>9s} {'Bonf p':>9s}")
for d in [0.5, 1.0, 2.0]:
    a = df[(df.supp == "OJ") & (df.dose == d)].len.mean()
    b = df[(df.supp == "VC") & (df.dose == d)].len.mean()
    diff = a - b
    t = diff / se
    p = 2 * stats.t.sf(abs(t), dfe)
    h = stats.t.ppf(0.975, dfe) * se
    print(f"{d:5.1f} {a:7.2f} {b:7.2f} {diff:8.2f}  "
          f"[{diff - h:7.2f},{diff + h:7.2f}] {p:9.4f} {min(1, 3 * p):9.4f}")

half = stats.t.ppf(0.975, dfe) * np.sqrt(MSE / n)
print(f"\n칸 평균의 SE = {np.sqrt(MSE / n):.4f}  (오차막대 ±{half:.4f})")

MSE = 13.1871,  df = 54,  SE(차이) = 1.6240

 dose      OJ      VC    OJ-VC               95% CI         p    Bonf p
  0.5   13.23    7.98     5.25  [   1.99,   8.51]    0.0021    0.0063
  1.0   22.70   16.77     5.93  [   2.67,   9.19]    0.0006    0.0018
  2.0   26.06   26.14    -0.08  [  -3.34,   3.18]    0.9609    1.0000

칸 평균의 SE = 1.1484  (오차막대 ±2.3023)