이원배치 Welch 분산분석 (로버스트 HC3)¶

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

# 칸마다 반복이 **둘 이상** 있어야 한다. 칸당 하나뿐이면 3x3 설계에서
# 모수 9개로 관측값 9개를 완전히 맞혀 버려 잔차 자유도가 0이 되고,
# 지렛값 h_ii가 1이 되어 HC3의 1/(1-h_ii)^2 이 발산한다.
data = {
    "Temperature": ["High"]*6 + ["Low"]*6 + ["Medium"]*6,
    "Fertilizer":  ["A", "A", "B", "B", "C", "C"] * 3,
    "Growth":      [12, 13, 15, 18, 14, 15,
                    10,  9, 13, 12, 11, 13,
                    14, 16, 16, 21, 15, 17],
}
df = pd.DataFrame(data)

model = ols("Growth ~ C(Temperature) * C(Fertilizer)", data=df).fit()
rob = model.get_robustcov_results(cov_type="HC3")

# Temperature의 주효과 검정.
# ":"가 든 이름을 빼야 한다. 교호작용 항의 이름도 "C(Temperature)["로 시작하므로
# 그냥 startswith만 쓰면 교호작용까지 함께 검정해 자유도가 2가 아니라 6이 된다.
pnames = model.params.index.tolist()
temp_params = [p for p in pnames
               if p.startswith("C(Temperature)[") and ":" not in p]
constraint = ", ".join([f"{t} = 0" for t in temp_params])
print("Main effect: Temperature")
print(rob.f_test(constraint))

# 교호작용 검정
inter_params = [p for p in pnames if ":" in p]
constraint_inter = ", ".join([f"{t} = 0" for t in inter_params])
print("Interaction: Temperature x Fertilizer")
print(rob.f_test(constraint_inter))

Main effect: Temperature
<F test: F=8.055555555555552, p=0.009878581991016048, df_denom=9, df_num=2>
Interaction: Temperature x Fertilizer
<F test: F=0.15370680044593157, p=0.956506962060376, df_denom=9, df_num=4>

import numpy as np

# 모든 지렛값이 정말 1/2 인가.
h = model.get_influence().hat_matrix_diag
print(f"지렛값 h_ii: 최소 {h.min():.6f}, 최대 {h.max():.6f}  "
      f"(모두 같은가: {np.allclose(h, 0.5)})")
print(f"HC3 가중 1/(1-h)^2 = {1 / (1 - h[0]) ** 2:.1f}")

# 칸평균 모수화로 다시 적합해 HC3 분산을 칸별 s^2/n 과 견준다.
cm = ols("Growth ~ C(Temperature):C(Fertilizer) - 1", data=df).fit()
cmr = cm.get_robustcov_results(cov_type="HC3")
hc3 = np.diag(cmr.cov_params())

cells = df.groupby(['Fertilizer', 'Temperature'])['Growth'].agg(['mean', 'var', 'count'])
classic = (cells['var'] / cells['count']).values
print(f"\n{'cell':<14}{'mean':>7}{'s^2':>8}{'s^2/n':>8}{'HC3':>8}{'비':>7}")
for (f, t), c, hv in zip(cells.index, classic, hc3):
    v = cells.loc[(f, t), 'var']
    print(f"{t + '-' + f:<14}{cells.loc[(f, t), 'mean']:>7.1f}{v:>8.2f}"
          f"{c:>8.2f}{hv:>8.2f}{hv / c:>7.2f}")
print(f"\n비가 모두 2 인가: {np.allclose(hc3 / classic, 2.0)}")

지렛값 h_ii: 최소 0.500000, 최대 0.500000  (모두 같은가: True)
HC3 가중 1/(1-h)^2 = 4.0

cell             mean     s^2   s^2/n     HC3      비
High-A           12.5    0.50    0.25    0.50   2.00
Low-A             9.5    0.50    0.25    0.50   2.00
Medium-A         15.0    2.00    1.00    2.00   2.00
High-B           16.5    4.50    2.25    4.50   2.00
Low-B            12.5    0.50    0.25    0.50   2.00
Medium-B         18.5   12.50    6.25   12.50   2.00
High-C           14.5    0.50    0.25    0.50   2.00
Low-C            12.0    2.00    1.00    2.00   2.00
Medium-C         16.0    2.00    1.00    2.00   2.00

비가 모두 2 인가: True

import statsmodels.api as sm

# 같은 자료를 보통의 분산분석으로 돌려 견준다. 등분산이 깨진 설계에서
# 두 방법의 p-값이 얼마나 갈리는지가 요점이다.
print(sm.stats.anova_lm(model, typ=2))

                                 sum_sq   df      F    PR(>F)
C(Temperature)                81.444444  2.0  14.66  0.001475
C(Fertilizer)                 36.777778  2.0   6.62  0.017060
C(Temperature):C(Fertilizer)   2.555556  4.0   0.23  0.914666
Residual                      25.000000  9.0    NaN       NaN

from scipy import stats

MSE = 25.0 / 9
cells = df.groupby(['Fertilizer', 'Temperature'])['Growth'].agg(['mean', 'var'])
print(f"합동 MSE = {MSE:.4f}   (칸 분산은 {cells['var'].min():.2f} 에서 "
      f"{cells['var'].max():.2f} 까지)")

print(f"\n{'비료':>5}{'평균':>8}{'Var(합동)':>12}{'Var(칸별)':>12}{'Var(HC3)':>11}")
for f in ('A', 'B', 'C'):
    v_cells = cells.loc[f, 'var'].values
    var_pooled = MSE / 6
    var_cell = (v_cells / 2).sum() / 9
    print(f"{f:>5}{cells.loc[f, 'mean'].mean():>8.3f}{var_pooled:>12.4f}"
          f"{var_cell:>12.4f}{2 * var_cell:>11.4f}")

# A 와 B 의 차에 대한 t 를 두 방식으로
mA, mB = cells.loc['A', 'mean'].mean(), cells.loc['B', 'mean'].mean()
vA = (cells.loc['A', 'var'].values / 2).sum() / 9
vB = (cells.loc['B', 'var'].values / 2).sum() / 9
se_pooled = np.sqrt(2 * MSE / 6)
se_hc3 = np.sqrt(2 * (vA + vB))
print(f"\nB - A = {mB - mA:.4f}")
print(f"  합동 SE = {se_pooled:.4f}  ->  t = {(mB - mA) / se_pooled:.4f}")
print(f"  HC3  SE = {se_hc3:.4f}  ->  t = {(mB - mA) / se_hc3:.4f}")

# 세 항의 F 를 나란히
pnames = model.params.index.tolist()
terms = {
    "Temperature": [p for p in pnames
                    if p.startswith("C(Temperature)[") and ":" not in p],
    "Fertilizer": [p for p in pnames
                   if p.startswith("C(Fertilizer)[") and ":" not in p],
    "Interaction": [p for p in pnames if ":" in p]}
aov = sm.stats.anova_lm(model, typ=2)
print(f"\n{'항':>12}{'고전 F':>10}{'고전 p':>10}{'HC3 F':>10}{'HC3 p':>10}")
for name, ps in terms.items():
    ft = rob.f_test(", ".join(f"{t} = 0" for t in ps))
    row = aov.index[list(terms).index(name)]
    print(f"{name:>12}{aov.loc[row, 'F']:>10.4f}{aov.loc[row, 'PR(>F)']:>10.4f}"
          f"{float(ft.fvalue):>10.4f}{float(ft.pvalue):>10.4f}")
print(f"\n5% 임계값 F(2,9) = {stats.f(2, 9).ppf(0.95):.4f}")

합동 MSE = 2.7778   (칸 분산은 0.50 에서 12.50 까지)

   비료      평균     Var(합동)     Var(칸별)   Var(HC3)
    A  12.333      0.4630      0.1667     0.3333
    B  15.833      0.4630      0.9722     1.9444
    C  14.167      0.4630      0.2500     0.5000

B - A = 3.5000
  합동 SE = 0.9623  ->  t = 3.6373
  HC3  SE = 1.5092  ->  t = 2.3191

           항      고전 F      고전 p     HC3 F     HC3 p
 Temperature   14.6600    0.0015    8.0556    0.0099
  Fertilizer    6.6200    0.0171    2.9474    0.1036
 Interaction    0.2300    0.9147    0.1537    0.9565

5% 임계값 F(2,9) = 4.2565

import warnings
warnings.filterwarnings("ignore")

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

def wald(fit, pick):
    idx = list(fit.params.index)
    names = [x for x in idx if pick(x)]
    R = np.zeros((len(names), len(idx)))
    for r, x in enumerate(names):
        R[r, idx.index(x)] = 1
    w = fit.f_test(R)
    return float(w.fvalue), int(w.df_num), int(w.df_denom), float(w.pvalue)

data = {
    "Temperature": ["High"] * 6 + ["Low"] * 6 + ["Medium"] * 6,
    "Fertilizer": ["A", "A", "B", "B", "C", "C"] * 3,
    "Growth": [12, 13, 15, 18, 14, 15,
               10, 9, 13, 12, 11, 13,
               14, 16, 16, 21, 15, 17],
}
df = pd.DataFrame(data)

print("=== 본문 자료 (3×3, 칸당 2, 교호작용 거의 없음) ===")
for lab, form, pre in [
    ("처리 대비 (기본)", "Growth ~ C(Temperature)*C(Fertilizer)",
     "C(Temperature)["),
    ("합 대비", "Growth ~ C(Temperature, Sum)*C(Fertilizer, Sum)",
     "C(Temperature, Sum)["),
]:
    fit = ols(form, data=df).fit(cov_type="HC3")
    F, d1, d2, p = wald(fit, lambda s, pre=pre: s.startswith(pre) and ":" not in s)
    Fi, i1, i2, pi = wald(fit, lambda s: ":" in s)
    print(f"  {lab:16s} Temp 주효과 F({d1},{d2}) = {F:7.4f}, p = {p:.4f}"
          f"   |  교호작용 F({i1},{i2}) = {Fi:6.4f}, p = {pi:.4f}")

# 교호작용이 뚜렷한 자료에서 다시
rng = np.random.default_rng(515)
mu = np.array([[10.0, 12.0, 14.0], [10.0, 13.0, 22.0]])
sd = np.array([[1.0, 1.5, 4.0], [1.0, 2.0, 5.0]])
rows = []
for i in range(2):
    for j in range(3):
        rows.append(pd.DataFrame({"y": rng.normal(mu[i, j], sd[i, j], 10),
                                  "A": f"A{i}", "B": f"B{j}"}))
d2f = pd.concat(rows, ignore_index=True)

print("\n=== 교호작용이 뚜렷한 2×3 (칸당 10, 균형) ===")
print(d2f.groupby(["A", "B"]).y.agg(["mean", "var"]).round(3).to_string())
for lab, form, pre in [("처리 대비 (기본)", "y ~ C(A)*C(B)", "C(A)["),
                       ("합 대비", "y ~ C(A, Sum)*C(B, Sum)", "C(A, Sum)[")]:
    fit = ols(form, data=d2f).fit(cov_type="HC3")
    F, d1, dd, p = wald(fit, lambda s, pre=pre: s.startswith(pre) and ":" not in s)
    print(f"  {lab:16s} A 주효과 F({d1},{dd}) = {F:8.4f}, p = {p:.4f}")

m = d2f.groupby(["A", "B"]).y.mean()
print(f"\n  A 의 단순효과: B0 {m[('A1', 'B0')] - m[('A0', 'B0')]:+.3f}, "
      f"B1 {m[('A1', 'B1')] - m[('A0', 'B1')]:+.3f}, "
      f"B2 {m[('A1', 'B2')] - m[('A0', 'B2')]:+.3f}")
print(f"  비가중 주변평균 차이: {m['A1'].mean() - m['A0'].mean():+.3f}")

=== 본문 자료 (3×3, 칸당 2, 교호작용 거의 없음) ===
  처리 대비 (기본)       Temp 주효과 F(2,9) =  8.0556, p = 0.0099   |  교호작용 F(4,9) = 0.1537, p = 0.9565
  합 대비             Temp 주효과 F(2,9) =  9.3094, p = 0.0064   |  교호작용 F(4,9) = 0.1537, p = 0.9565

=== 교호작용이 뚜렷한 2×3 (칸당 10, 균형) ===
         mean     var
A  B
A0 B0   9.178   1.742
   B1  12.621   2.614
   B2  12.117  13.669
A1 B0  10.180   1.113
   B1  13.729   5.214
   B2  24.821  32.266
  처리 대비 (기본)       A 주효과 F(1,54) =   3.1623, p = 0.0810
  합 대비             A 주효과 F(1,54) =  34.8784, p = 0.0000

  A 의 단순효과: B0 +1.002, B1 +1.107, B2 +12.704
  비가중 주변평균 차이: +4.938

B0 에서  +1.00   ← 처리 대비가 검정하는 것
B1 에서  +1.11
B2 에서 +12.70
------------------
평균     +4.94   ← 합 대비가 검정하는 것 (주효과)

import warnings
warnings.filterwarnings("ignore")

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

def wald_ps(df, cov):
    fit = (ols("y ~ C(A, Sum)*C(B, Sum)", data=df).fit(cov_type=cov)
           if cov else ols("y ~ C(A, Sum)*C(B, Sum)", data=df).fit())
    idx = list(fit.params.index)
    ps = []
    for pick in [lambda s: s.startswith("C(A, Sum)[") and ":" not in s,
                 lambda s: s.startswith("C(B, Sum)[") and ":" not in s,
                 lambda s: ":" in s]:
        nm = [x for x in idx if pick(x)]
        R = np.zeros((len(nm), len(idx)))
        for r, x in enumerate(nm):
            R[r, idx.index(x)] = 1
        ps.append(float(fit.f_test(R).pvalue))
    return np.array(ps)

NS = np.array([[8, 6, 4], [7, 5, 4]])
SD = np.array([[1.0, 2.0, 5.0], [1.0, 2.0, 5.0]])
rng = np.random.default_rng(4242)
B = 1_200
meths = [("고전(OLS)", None), ("HC0", "HC0"), ("HC1", "HC1"),
         ("HC2", "HC2"), ("HC3", "HC3")]
cnt = {k: np.zeros(3) for k, _ in meths}
for _ in range(B):
    rows = []
    for i in range(2):
        for j in range(3):
            rows.append(pd.DataFrame({"y": rng.normal(0, SD[i, j], NS[i, j]),
                                      "A": f"A{i}", "B": f"B{j}"}))
    df = pd.concat(rows, ignore_index=True)
    for k, c in meths:
        cnt[k] += wald_ps(df, c) < 0.05

print("2×3 완전 귀무, 칸 n=[[8,6,4],[7,5,4]], σ=[[1,2,5],[1,2,5]], B=1200")
print(f"{'방법':>10s} {'A 주효과':>9s} {'B 주효과':>9s} {'A×B':>9s}")
for k, _ in meths:
    x = cnt[k]
    print(f"{k:>10s} {x[0] / B:9.4f} {x[1] / B:9.4f} {x[2] / B:9.4f}")

2×3 완전 귀무, 칸 n=[[8,6,4],[7,5,4]], σ=[[1,2,5],[1,2,5]], B=1200
        방법     A 주효과     B 주효과       A×B
   고전(OLS)    0.1458    0.1767    0.1442
       HC0    0.1075    0.1375    0.1208
       HC1    0.0725    0.1108    0.0892
       HC2    0.0683    0.1033    0.0850
       HC3    0.0475    0.0583    0.0550

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

for n in [1, 2, 3, 5, 10]:
    rows = []
    for i in range(2):
        for j in range(3):
            rows.append(pd.DataFrame({"y": np.arange(n, dtype=float),
                                      "A": f"A{i}", "B": f"B{j}"}))
    df = pd.concat(rows, ignore_index=True)
    X = ols("y ~ C(A)*C(B)", data=df).fit().model.exog
    H = X @ np.linalg.pinv(X.T @ X) @ X.T
    h = np.diag(H)
    hc3 = "∞" if n == 1 else f"{1 / (1 - h[0])**2:.4f}"
    hc2 = "∞" if n == 1 else f"{1 / (1 - h[0]):.4f}"
    print(f"  칸당 n={n:2d}: h_ii = {h[0]:.4f} (이론 1/n = {1 / n:.4f}), "
          f"HC3 배율 1/(1-h)² = {hc3}, HC2 배율 = {hc2}")

  칸당 n= 1: h_ii = 1.0000 (이론 1/n = 1.0000), HC3 배율 1/(1-h)² = ∞, HC2 배율 = ∞
  칸당 n= 2: h_ii = 0.5000 (이론 1/n = 0.5000), HC3 배율 1/(1-h)² = 4.0000, HC2 배율 = 2.0000
  칸당 n= 3: h_ii = 0.3333 (이론 1/n = 0.3333), HC3 배율 1/(1-h)² = 2.2500, HC2 배율 = 1.5000
  칸당 n= 5: h_ii = 0.2000 (이론 1/n = 0.2000), HC3 배율 1/(1-h)² = 1.5625, HC2 배율 = 1.2500
  칸당 n=10: h_ii = 0.1000 (이론 1/n = 0.1000), HC3 배율 1/(1-h)² = 1.2346, HC2 배율 = 1.1111

import warnings
warnings.filterwarnings("ignore")

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

def welch_james(m, v, n, Cm):
    m, v, n = np.asarray(m, float), np.asarray(v, float), np.asarray(n, float)
    V = np.diag(v / n)
    Q = Cm.T @ np.linalg.inv(Cm @ V @ Cm.T) @ Cm
    T = m @ Q @ m
    q = np.linalg.matrix_rank(Cm)
    VQ = V @ Q
    A = 0.0
    for i in range(len(n)):
        E = np.zeros_like(V); E[i, i] = 1
        M = VQ @ E
        A += 0.5 * (np.trace(M)**2 + np.trace(M @ M)) / (n[i] - 1)
    c = q + 2 * A - 6 * A / (q + 2)
    return stats.f.sf(T / c, q, q * (q + 2) / (3 * A))

CAB = np.kron(np.array([[1.0, -1]]),
              np.array([[1.0, -1, 0], [0, 1, -1]]))     # 교호작용 대비

NS = np.array([[10, 8, 6], [9, 8, 6]])
SD = np.array([[1.0, 2.0, 4.0], [1.0, 2.0, 4.0]])
MU0 = np.zeros((2, 3))
MU1 = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 4.0]])      # 교호작용 있음

def run(MU, B, seed):
    rng = np.random.default_rng(seed)
    wj = hc = 0
    for _ in range(B):
        rows = []
        for i in range(2):
            for j in range(3):
                rows.append(pd.DataFrame({"y": rng.normal(MU[i, j], SD[i, j],
                                                          NS[i, j]),
                                          "A": f"A{i}", "B": f"B{j}"}))
        df = pd.concat(rows, ignore_index=True)
        cell = df.groupby(["A", "B"]).y.agg(["size", "mean", "var"])
        wj += welch_james(cell["mean"].values, cell["var"].values,
                          cell["size"].values, CAB) < 0.05
        fit = ols("y ~ C(A, Sum)*C(B, Sum)", data=df).fit(cov_type="HC3")
        idx = list(fit.params.index)
        nm = [x for x in idx if ":" in x]
        R = np.zeros((len(nm), len(idx)))
        for r, x in enumerate(nm):
            R[r, idx.index(x)] = 1
        hc += float(fit.f_test(R).pvalue) < 0.05
    return wj / B, hc / B

B = 1_500
s_wj, s_hc = run(MU0, B, 808)
p_wj, p_hc = run(MU1, B, 909)
print(f"교호작용 검정, 2×3, n=[[10,8,6],[9,8,6]], "
      f"σ=[[1,2,4],[1,2,4]], B={B}")
print(f"{'':12s} {'크기':>8s} {'검정력':>8s}")
print(f"{'웰치-제임스':>12s} {s_wj:8.4f} {p_wj:8.4f}")
print(f"{'HC3':>12s} {s_hc:8.4f} {p_hc:8.4f}")

교호작용 검정, 2×3, n=[[10,8,6],[9,8,6]], σ=[[1,2,4],[1,2,4]], B=1500
                   크기      검정력
      웰치-제임스   0.0420   0.2607
         HC3   0.0427   0.2593

1. 합 대비로 완전 요인 모형 (필수)
2. HC3 로버스트 공분산 (필수)
     fit = ols("y ~ C(A, Sum)*C(B, Sum)", data=df).fit(cov_type="HC3")

3. 항별 왈드 F — 이름 선택에 ":" 조건을 반드시 넣는다
     A 주효과 : startswith("C(A, Sum)[") and ":" not in name
     교호작용 : ":" in name

2×3 요인설계, 칸별 n = 6~10
칸 분산이 1.1~32.3 (29배) → 등분산 가정 불가
합 대비 완전 요인 OLS + HC3 로버스트 공분산, 왈드 F 검정

  A 주효과   F(1, 54) = 34.88,  p < 0.001
  B 주효과   F(2, 54) = ...
  A×B       F(2, 54) = ...

교호작용이 유의하므로 단순효과를 본페로니 보정하여 보고한다.