scipy를 이용한 일원배치 분산분석과 그림¶

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

url = ('https://raw.githubusercontent.com/vincentarelbundock/'
       'Rdatasets/1dcc2bf5f955cc1224a3e1307256e1fe86b68dae/csv/datasets/PlantGrowth.csv')
df = pd.read_csv(url, usecols=[1, 2])
g = df.groupby('group')
# f_oneway는 집단을 **별도의 배열**로 받는다. 긴 형식 데이터프레임을 그대로
# 넘길 수 없어서 이렇게 쪼개야 한다. statsmodels의 ols 방식은 그 반대다.
ctrl = g.get_group('ctrl').weight.values
trt1 = g.get_group('trt1').weight.values
trt2 = g.get_group('trt2').weight.values

print(df.groupby('group').weight.agg(['count', 'mean', 'std']).round(4))

groups = [ctrl, trt1, trt2]
n = np.array([len(x) for x in groups])
mean_i = np.array([x.mean() for x in groups])
var_i = np.array([x.var(ddof=1) for x in groups])
N, k = n.sum(), len(groups)
grand = df.weight.mean()

# 제곱합 분해. 세 수를 따로 계산하고 항등식이 맞는지 본다.
SSB = (n * (mean_i - grand) ** 2).sum()
SSW = ((n - 1) * var_i).sum()
SST = ((df.weight - grand) ** 2).sum()
print(f"\nSSB = {SSB:.6f}   SSW = {SSW:.6f}   합 = {SSB + SSW:.6f}")
print(f"SST = {SST:.6f}   차이 = {abs(SST - SSB - SSW):.3e}")
print(f"자유도  {N - 1} = {k - 1} + {N - k}")

# 균형 설계이므로 MSW 는 표본분산의 단순평균이어야 한다.
MSW = SSW / (N - k)
print(f"\nMSW(=SSW/(N-k)) = {MSW:.8f}")
print(f"표본분산의 단순평균 = {var_i.mean():.8f}")

# 요약표의 반올림된 표준편차만 쓰면 얼마나 어긋나는가.
sd_rounded = np.array([0.5831, 0.7937, 0.4426])
print(f"반올림 표준편차로 = {(sd_rounded ** 2).mean():.8f}")
print(f"합동표준편차 sqrt(MSW) = {np.sqrt(MSW):.4f}")

       count   mean     std
group                      
ctrl      10  5.032  0.5831
trt1      10  4.661  0.7937
trt2      10  5.526  0.4426

SSB = 3.766340   SSW = 10.492090   합 = 14.258430
SST = 14.258430   차이 = 3.553e-15
자유도  29 = 2 + 27

MSW(=SSW/(N-k)) = 0.38859593
표본분산의 단순평균 = 0.38859593
반올림 표준편차로 = 0.38862002
합동표준편차 sqrt(MSW) = 0.6234

# 보기 1 의 ctrl, trt1, trt2, n, mean_i, MSW, k 를 이어 쓴다.
# F 는 집단 사이의 분산을 집단 안의 분산으로 나눈 값이다. 1 에 가까우면
# 집단을 나눈 것이 아무 설명도 하지 못한다는 뜻이다.
F, p = stats.f_oneway(ctrl, trt1, trt2)
print(f"F = {F:.4f}, p = {p:.4f}")

# 요약통계(개수와 평균)만으로 손계산한다.
grand = (n * mean_i).sum() / n.sum()
SSB_hand = (n * (mean_i - grand) ** 2).sum()
MSB_hand = SSB_hand / (k - 1)
print(f"\n전체평균 = {grand:.4f}")
print(f"SSB = {SSB_hand:.6f}   MSB = {MSB_hand:.6f}   MSW = {MSW:.6f}")
print(f"F(손계산) = {MSB_hand / MSW:.10f}")
print(f"F(scipy)  = {F:.10f}")

# k=2 이면 F = t^2 이어야 한다. ctrl 과 trt1 둘만 쓴다.
F2, p2 = stats.f_oneway(ctrl, trt1)
t2, pt = stats.ttest_ind(ctrl, trt1)
print(f"\nk=2:  F = {F2:.10f}   t = {t2:.10f}   t^2 = {t2 ** 2:.10f}")
print(f"      p(F) = {p2:.10f}   p(t) = {pt:.10f}")
print(f"\nE[F(2,27)] = {27 / 25:.4f}")

F = 4.8461, p = 0.0159

전체평균 = 5.0730
SSB = 3.766340   MSB = 1.883170   MSW = 0.388596
F(손계산) = 4.8460878624
F(scipy)  = 4.8460878624

k=2:  F = 1.4191012974   t = 1.1912603818   t^2 = 1.4191012974
      p(F) = 0.2490231660   p(t) = 0.2490231660

E[F(2,27)] = 1.0800

import matplotlib.pyplot as plt

# 상자그림이 실제로 그리는 수들을 먼저 찍어 둔다.
names = ['ctrl', 'trt1', 'trt2']
groups = [ctrl, trt1, trt2]
print(f"{'group':>6}{'mean':>8}{'median':>8}{'Q1':>8}{'Q3':>8}{'IQR':>8}{'sd':>8}")
for name, x in zip(names, groups):
    q1, q2, q3 = np.percentile(x, [25, 50, 75])
    print(f"{name:>6}{x.mean():>8.3f}{q2:>8.3f}{q1:>8.3f}{q3:>8.3f}"
          f"{q3 - q1:>8.3f}{x.std(ddof=1):>8.3f}")

# 상자 바깥으로 따로 찍히는 점이 어느 집단에 있는가.
for name, x in zip(names, groups):
    q1, q3 = np.percentile(x, [25, 75])
    lo, hi = q1 - 1.5 * (q3 - q1), q3 + 1.5 * (q3 - q1)
    flagged = np.sort(x[(x < lo) | (x > hi)])
    print(f"{name}: 울타리 밖 {flagged if len(flagged) else '없음'}")

# 그림에서는 읽히지 않는 것 — 쌍별 t-검정의 p-값.
for i in range(3):
    for j in range(i + 1, 3):
        t, p = stats.ttest_ind(groups[i], groups[j])
        print(f"{names[i]} vs {names[j]}: t = {t:+.4f}, p = {p:.4f}")

# 검정을 하기 전이든 뒤든 그림을 본다. 상자가 겹치는 정도가 F 값과
# 어떻게 맞물리는지 눈에 익혀 두면 좋다.
fig, ax = plt.subplots(figsize=(6, 4))
ax.boxplot([ctrl, trt1, trt2], labels=['ctrl', 'trt1', 'trt2'])
ax.set_xlabel('Group')
ax.set_ylabel('Weight')
ax.set_title('Plant weights by group')
plt.tight_layout()
plt.show()

 group    mean  median      Q1      Q3     IQR      sd
  ctrl   5.032   5.155   4.550   5.293   0.742   0.583
  trt1   4.661   4.550   4.207   4.870   0.662   0.794
  trt2   5.526   5.435   5.268   5.735   0.467   0.443
ctrl: 울타리 밖 없음
trt1: 울타리 밖 [5.87 6.03]
trt2: 울타리 밖 없음
ctrl vs trt1: t = +1.1913, p = 0.2490
ctrl vs trt2: t = -2.1340, p = 0.0469
trt1 vs trt2: t = -3.0101, p = 0.0075

import numpy as np

# 자유도는 (집단 수 - 1, 전체 수 - 집단 수) = (2, 27) 이다.
m = 27

# 닫힌 꼴 꼬리확률과 scipy 를 맞춰 본다.
p_closed = (1 + 2 * F / m) ** (-m / 2)
p_scipy = stats.f(2, m).sf(F)
print(f"닫힌 꼴 p = {p_closed:.10f}")
print(f"scipy   p = {p_scipy:.10f}")
print(f"차이      = {abs(p_closed - p_scipy):.3e}")

# 밀도도 닫힌 꼴이다. x=0 에서 1 이고 단조감소한다.
x = np.array([0.0, 0.5, 1.0, 2.0, 3.0, F, 8.0])
pdf_closed = (1 + 2 * x / m) ** (-(m + 2) / 2)
pdf_scipy = stats.f(2, m).pdf(x)
print(f"\n{'x':>8}{'closed':>14}{'scipy':>14}")
for xi, a, b in zip(x, pdf_closed, pdf_scipy):
    print(f"{xi:>8.3f}{a:>14.8f}{b:>14.8f}")

# 임계값도 뒤집어 풀린다.
crit_closed = m / 2 * (0.05 ** (-2 / m) - 1)
print(f"\n임계값 닫힌 꼴 = {crit_closed:.8f}")
print(f"임계값 scipy   = {stats.f(2, m).ppf(0.95):.8f}")
print(f"E[F(2,27)] = {m / (m - 2):.4f},  최빈값 = 0 (밀도가 단조감소)")

# 칠해진 오른쪽 꼬리의 넓이가 곧 p-값이다.
x = np.linspace(0, 8, 400)
pdf = stats.f(2, 27).pdf(x)

fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(x, pdf, label='F(2, 27) PDF')
mask = x >= F
ax.fill_between(x[mask], pdf[mask], alpha=0.3, label='Observed tail')
ax.set_title('F-distribution and observed tail')
ax.legend()
plt.tight_layout()
plt.show()

닫힌 꼴 p = 0.0159099583
scipy   p = 0.0159099583
차이      = 2.776e-17

       x        closed         scipy
   0.000    1.00000000    1.00000000
   0.500    0.59017815    0.59017815
   1.000    0.35481634    0.35481634
   2.000    0.13490561    0.13490561
   3.000    0.05449071    0.05449071
   4.846    0.01170737    0.01170737
   8.000    0.00117350    0.00117350

임계값 닫힌 꼴 = 3.35413083
임계값 scipy   = 3.35413083
E[F(2,27)] = 1.0800,  최빈값 = 0 (밀도가 단조감소)

import numpy as np
from scipy import stats

print(f"{'df1':>4s} {'df2':>5s} {'평균':>8s} {'최빈값':>9s} "
      f"{'왜도':>8s} {'F(0.95)':>9s}")
for d1, d2 in [(1, 27), (2, 27), (5, 27), (2, 10), (2, 100),
               (10, 100), (50, 100)]:
    mean = d2 / (d2 - 2)
    mode = (d1 - 2) / d1 * d2 / (d2 + 2) if d1 > 2 else np.nan
    skew = float(stats.f.stats(d1, d2, moments='s'))
    print(f"{d1:4d} {d2:5d} {mean:8.4f} {mode:9.4f} {skew:8.4f} "
          f"{stats.f.ppf(0.95, d1, d2):9.4f}")

 df1   df2       평균       최빈값       왜도   F(0.95)
   1    27   1.0800       nan   3.4203    4.2100
   2    27   1.0800       nan   2.5491    3.3541
   5    27   1.0800    0.5586   1.8459    2.5719
   2    10   1.2500       nan   4.6476    4.1028
   2   100   1.0204       nan   2.1264    3.0873
  10   100   1.0204    0.7843   1.0586    1.9267
  50   100   1.0204    0.9412   0.6786    1.4772

for d1 in [1, 2, 5, 10]:
    print(f"d1={d1:2d}:  F(0.95, d1, ∞) = "
          f"{stats.chi2.ppf(0.95, d1) / d1:.4f}   "
          f"F(0.95, d1, 1000) = {stats.f.ppf(0.95, d1, 1000):.4f}")

d1= 1:  F(0.95, d1, ∞) = 3.8415   F(0.95, d1, 1000) = 3.8508
d1= 2:  F(0.95, d1, ∞) = 2.9957   F(0.95, d1, 1000) = 3.0047
d1= 5:  F(0.95, d1, ∞) = 2.2141   F(0.95, d1, 1000) = 2.2231
d1=10:  F(0.95, d1, ∞) = 1.8307   F(0.95, d1, 1000) = 1.8402

import numpy as np
from scipy import stats

rng = np.random.default_rng(555)
for n1, n2 in [(10, 10), (15, 12)]:
    x = rng.normal(0, 1, n1)
    y = rng.normal(0.6, 1, n2)
    F, p_F = stats.f_oneway(x, y)
    t, p_t = stats.ttest_ind(x, y)
    print(f"n = ({n1},{n2}):  F = {F:.6f},  t² = {t**2:.6f}")
    print(f"              p_F = {p_F:.8f},  p_t = {p_t:.8f}")
    print(f"              t = {t:+.6f}  ← 부호가 방향을 알려준다")

n = (10,10):  F = 2.610177,  t² = 2.610177
              p_F = 0.12357147,  p_t = 0.12357147
              t = -1.615604  ← 부호가 방향을 알려준다
n = (15,12):  F = 2.912193,  t² = 2.912193
              p_F = 0.10030624,  p_t = 0.10030624
              t = -1.706515  ← 부호가 방향을 알려준다

x = rng.normal(0, 1, 20)
y = rng.normal(0.7, 1, 20)
print(f"양측      p = {stats.ttest_ind(x, y).pvalue:.4f}")
print(f"단측(x<y) p = {stats.ttest_ind(x, y, alternative='less').pvalue:.4f}")
print(f"F 검정    p = {stats.f_oneway(x, y).pvalue:.4f}  ← 언제나 양측")

양측      p = 0.1230
단측(x<y) p = 0.0615
F 검정    p = 0.1230  ← 언제나 양측

diff = x.mean() - y.mean()
n1 = n2 = 20
sp = np.sqrt(((n1 - 1) * x.var(ddof=1) + (n2 - 1) * y.var(ddof=1))
             / (n1 + n2 - 2))
se = sp * np.sqrt(1 / n1 + 1 / n2)
tc = stats.t.ppf(0.975, n1 + n2 - 2)
print(f"평균 차이 {diff:+.4f},  95% CI "
      f"({diff - tc * se:+.4f}, {diff + tc * se:+.4f})")

평균 차이 -0.5079,  95% CI (-1.1596, +0.1439)

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

df = pd.DataFrame({
    "group": ["ctrl"] * 10 + ["trt1"] * 10 + ["trt2"] * 10,
    "weight": [4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14,
               4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69,
               6.31, 5.12, 5.54, 5.50, 5.37, 5.29, 4.92, 6.15, 5.80, 5.26],
})

# ① 올바른 사용: 집단을 별도의 배열로 쪼개 넘긴다
groups = [v.values for _, v in df.groupby("group")["weight"]]
F, p = stats.f_oneway(*groups)
print(f"① 올바름:              F = {F:.6f},  p = {p:.6f}")

# ② 흔한 실수 1: 두 열을 그대로 넘긴다
try:
    stats.f_oneway(df["weight"], df["group"])
except Exception as e:
    print(f"② 두 열을 넘김:        {type(e).__name__}: {str(e)[:60]}")

# ③ 흔한 실수 2: 리스트 하나로 감싼다 (별표를 빠뜨림)
try:
    stats.f_oneway(groups)
except Exception as e:
    print(f"③ 별표를 빠뜨림:       {type(e).__name__}: {str(e)[:60]}")

# ④ 결측값
with_nan = [np.r_[groups[0], np.nan], groups[1], groups[2]]
print(f"④ NaN 포함(기본):      {stats.f_oneway(*with_nan)}")
F2, p2 = stats.f_oneway(*with_nan, nan_policy="omit")
print(f"   nan_policy='omit':  F = {F2:.6f},  p = {p2:.6f}")

① 올바름:              F = 4.846088,  p = 0.015910
② 두 열을 넘김:        TypeError: unsupported operand type(s) for +: 'float' and 'str'
③ 별표를 빠뜨림:       TypeError: at least two inputs are required; got 1.
④ NaN 포함(기본):      F_onewayResult(statistic=nan, pvalue=nan)
   nan_policy='omit':  F = 4.846088,  p = 0.015910

# 위험: 집단 순서가 정렬 순서에 의존한다
groups_a = [v.values for _, v in df.groupby("group")["weight"]]
# 안전: 순서를 명시한다
order = ["ctrl", "trt1", "trt2"]
groups_b = [df.loc[df["group"] == g, "weight"].values for g in order]
print(f"같은 결과인가: "
      f"{np.isclose(stats.f_oneway(*groups_a).statistic, stats.f_oneway(*groups_b).statistic)}")
print(f"집단 크기: {[len(g) for g in groups_b]}")
print(f"집단 이름: {order}")

같은 결과인가: True
집단 크기: [10, 10, 10]
집단 이름: ['ctrl', 'trt1', 'trt2']

import numpy as np
from scipy import stats

groups = {
    "ctrl": np.array([4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14]),
    "trt1": np.array([4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69]),
    "trt2": np.array([6.31, 5.12, 5.54, 5.50, 5.37, 5.29, 4.92, 6.15, 5.80, 5.26]),
}
g = list(groups.values())

def welch_anova(*grp):
    k = len(grp)
    n = np.array([len(x) for x in grp], float)
    m = np.array([x.mean() for x in grp])
    v = np.array([x.var(ddof=1) for x in grp])
    w = n / v
    W = w.sum()
    tmp = np.sum((1 - w / W)**2 / (n - 1))
    m_tilde = (w * m).sum() / W
    F = ((w * (m - m_tilde)**2).sum() / (k - 1)) \
        / (1 + 2 * (k - 2) / (k * k - 1) * tmp)
    return F, stats.f.sf(F, k - 1, (k * k - 1) / (3 * tmp))

rng = np.random.default_rng(2468)

def perm_f(grp, B, rng):
    obs = stats.f_oneway(*grp).statistic
    sizes = [len(x) for x in grp]
    pooled = np.concatenate(grp)
    cnt = 1
    for _ in range(B):
        rng.shuffle(pooled)
        parts = np.split(pooled, np.cumsum(sizes)[:-1])
        cnt += stats.f_oneway(*parts).statistic >= obs - 1e-12
    return cnt / (B + 1)

print(f"집단 분산: "
      f"{[round(float(x.var(ddof=1)), 4) for x in g]}")
print(f"분산비 = {max(x.var(ddof=1) for x in g) / min(x.var(ddof=1) for x in g):.4f}\n")
F1, p1 = stats.f_oneway(*g)
F2, p2 = welch_anova(*g)
ag = stats.alexandergovern(*g)
H, p4 = stats.kruskal(*g)
print(f"고전 F            F = {F1:.4f},  p = {p1:.4f}")
print(f"Welch 분산분석     F = {F2:.4f},  p = {p2:.4f}")
print(f"알렉산더·고번      A = {ag.statistic:.4f},  p = {ag.pvalue:.4f}")
print(f"크러스컬·월리스    H = {H:.4f},  p = {p4:.4f}")
print(f"순열 F                            p = {perm_f(g, 9_999, rng):.4f}")
print(f"\n등분산 검정 (참고용, 검정 선택에는 쓰지 말 것)")
print(f"  레빈(중앙값)  p = {stats.levene(*g, center='median').pvalue:.4f}")
print(f"  바틀렛        p = {stats.bartlett(*g).pvalue:.4f}")

집단 분산: [0.34, 0.6299, 0.1959]
분산비 = 3.2160

고전 F            F = 4.8461,  p = 0.0159
Welch 분산분석     F = 5.1810,  p = 0.0174
알렉산더·고번      A = 8.3285,  p = 0.0155
크러스컬·월리스    H = 7.9882,  p = 0.0184
순열 F                            p = 0.0169

등분산 검정 (참고용, 검정 선택에는 쓰지 말 것)
  레빈(중앙값)  p = 0.3412
  바틀렛        p = 0.2371

import numpy as np
from scipy import stats

def oneway_report(groups, labels=None, alpha=0.05):
    """f_oneway 에 진단과 효과크기를 붙인 보고 함수."""
    g = [np.asarray(x, float) for x in groups]
    k = len(g)
    n = np.array([len(x) for x in g])
    N = n.sum()
    labels = labels or [f"집단{i + 1}" for i in range(k)]

    for lab, x in zip(labels, g):
        print(f"  {lab:>8s}: n={len(x):3d}  평균={x.mean():8.4f}  "
              f"표준편차={x.std(ddof=1):7.4f}")
    v = [x.var(ddof=1) for x in g]
    ratio = max(v) / min(v)
    print(f"  분산비 = {ratio:.4f}"
          + ("   ⚠ 4 초과 — Welch 를 고려하라" if ratio > 4 else ""))
    if n.max() / n.min() > 1.5:
        print(f"  ⚠ 표본크기 불균형 ({n.tolist()})")

    F, p = stats.f_oneway(*g)
    grand = np.concatenate(g).mean()
    SST = sum(len(x) * (x.mean() - grand)**2 for x in g)
    SSE = sum(((x - x.mean())**2).sum() for x in g)
    MSE = SSE / (N - k)
    print(f"\n  F({k - 1}, {N - k}) = {F:.4f},  p = {p:.6f}"
          f"   → {'기각' if p < alpha else '기각 못 함'}")
    print(f"  η² = {SST / (SST + SSE):.4f},  "
          f"ω² = {(SST - (k - 1) * MSE) / (SST + SSE + MSE):.4f}")
    return {"F": F, "p": p, "df1": k - 1, "df2": int(N - k),
            "eta2": SST / (SST + SSE)}

_ = oneway_report(
    [[4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14],
     [4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69],
     [6.31, 5.12, 5.54, 5.50, 5.37, 5.29, 4.92, 6.15, 5.80, 5.26]],
    labels=["ctrl", "trt1", "trt2"])

      ctrl: n= 10  평균=  5.0320  표준편차= 0.5831
      trt1: n= 10  평균=  4.6610  표준편차= 0.7937
      trt2: n= 10  평균=  5.5260  표준편차= 0.4426
  분산비 = 3.2160

  F(2, 27) = 4.8461,  p = 0.015910   → 기각
  η² = 0.2641,  ω² = 0.2041