사후비교: Tukey HSD¶

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

# R 의 PlantGrowth 자료. 대조군과 처리군 둘, 모두 세 집단이다.
url = 'https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/1dcc2bf5f955cc1224a3e1307256e1fe86b68dae/csv/datasets/PlantGrowth.csv'
df = pd.read_csv(url, usecols=[1, 2])

# 분산분석이 먼저다. 여기서 유의하지 않으면 사후비교로 넘어갈 까닭이 없다.
model = ols('weight ~ C(group)', data=df).fit()
anova_results = anova_lm(model)
print("One-Way ANOVA Results:")
print(anova_results)

One-Way ANOVA Results:
            df    sum_sq   mean_sq         F   PR(>F)
C(group)   2.0   3.76634  1.883170  4.846088  0.01591
Residual  27.0  10.49209  0.388596       NaN      NaN

import numpy as np
import pandas as pd

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

g = df.groupby('group').weight
m, s = g.mean().values, g.std().values
print(f"집단평균 = {np.round(m, 4)},  집단표준편차 = {np.round(s, 6)}")
print(f"전체평균 = {df.weight.mean():.6f}")

SSB_dev = n * ((m - m.mean()) ** 2).sum()
SSB_pair = (n / k) * sum((m[i] - m[j]) ** 2 for i in range(k) for j in range(i + 1, k))
MSE = (s ** 2).mean()
print(f"\nSSB (편차 꼴) = {SSB_dev:.6f}")
print(f"SSB (쌍  꼴) = {SSB_pair:.6f}")
print(f"MSE = 집단분산의 평균 = {MSE:.6f}")
print(f"F = (SSB/{k - 1})/MSE = {(SSB_dev / (k - 1)) / MSE:.6f}")

print(f"\n쌍별 제곱차가 SSB 에 기여하는 몫")
labels = ['ctrl', 'trt1', 'trt2']
for i in range(k):
    for j in range(i + 1, k):
        d2 = (m[i] - m[j]) ** 2
        print(f"  {labels[i]}-{labels[j]}: (차)^2 = {d2:.6f},  몫 = {(n / k) * d2 / SSB_dev:.1%}")

집단평균 = [5.032 4.661 5.526],  집단표준편차 = [0.583091 0.793676 0.442573]
전체평균 = 5.073000

SSB (편차 꼴) = 3.766340
SSB (쌍  꼴) = 3.766340
MSE = 집단분산의 평균 = 0.388596
F = (SSB/2)/MSE = 4.846088

쌍별 제곱차가 SSB 에 기여하는 몫
  ctrl-trt1: (차)^2 = 0.137641,  몫 = 12.2%
  ctrl-trt2: (차)^2 = 0.244036,  몫 = 21.6%
  trt1-trt2: (차)^2 = 0.748225,  몫 = 66.2%

from statsmodels.stats.multicomp import pairwise_tukeyhsd

# 분산분석은 "어딘가 다르다"까지만 말한다. 어느 쌍이 다른지는 사후비교의 몫이다.
tukey_result = pairwise_tukeyhsd(endog=df['weight'], groups=df['group'], alpha=0.05)
print("Tukey's HSD Test Results:")
print(tukey_result)

Tukey's HSD Test Results:
Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj   lower  upper  reject
---------------------------------------------------
  ctrl   trt1   -0.371 0.3909 -1.0622 0.3202  False
  ctrl   trt2    0.494  0.198 -0.1972 1.1852  False
  trt1   trt2    0.865  0.012  0.1738 1.5562   True
---------------------------------------------------

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])
n, k, nu = 10, 3, 27

g = df.groupby('group').weight
m, s = g.mean().values, g.std().values
MSE = (s ** 2).mean()
SE = np.sqrt(MSE * (1 / n + 1 / n))

q = stats.studentized_range.ppf(0.95, k, nu)
HSD = q * np.sqrt(MSE / n)
print(f"MSE = {MSE:.6f},  SE = sqrt(MSE*(1/n+1/n)) = {SE:.6f}")
print(f"q(0.95, k=3, nu=27) = {q:.6f}")
print(f"HSD = q*sqrt(MSE/n) = {HSD:.6f}   (= q/sqrt(2) * SE = {q / np.sqrt(2) * SE:.6f})")

print(f"\n{'쌍':>12}{'차이':>9}{'|차이|>HSD':>11}{'t':>9}{'p-adj':>9}{'하한':>9}{'상한':>9}")
labels = ['ctrl', 'trt1', 'trt2']
for i in range(k):
    for j in range(i + 1, k):
        d = m[j] - m[i]
        t = d / SE
        p = stats.studentized_range.sf(abs(t) * np.sqrt(2), k, nu)
        print(f"{labels[i] + '-' + labels[j]:>12}{d:>9.4f}{str(abs(d) > HSD):>11}"
              f"{t:>9.4f}{p:>9.4f}{d - HSD:>9.4f}{d + HSD:>9.4f}")

MSE = 0.388596,  SE = sqrt(MSE*(1/n+1/n)) = 0.278782
q(0.95, k=3, nu=27) = 3.506426
HSD = q*sqrt(MSE/n) = 0.691216   (= q/sqrt(2) * SE = 0.691216)

           쌍       차이   |차이|>HSD        t    p-adj       하한       상한
   ctrl-trt1  -0.3710      False  -1.3308   0.3909  -1.0622   0.3202
   ctrl-trt2   0.4940      False   1.7720   0.1980  -0.1972   1.1852
   trt1-trt2   0.8650       True   3.1028   0.0120   0.1738   1.5562

from statsmodels.stats.multitest import multipletests
from itertools import combinations
from scipy.stats import ttest_ind

# 같은 일을 본페로니로 해 본다. 쌍마다 t-검정을 하고 p-값에 비교 횟수를 곱한다.
groups = df['group'].unique()

# 쌍 세 개를 모두 돌며 보정 전 p-값을 모은다.
p_values = []
comparisons = []

for group1, group2 in combinations(groups, 2):
    data1 = df[df['group'] == group1]['weight']
    data2 = df[df['group'] == group2]['weight']
    stat, p_val = ttest_ind(data1, data2)
    p_values.append(p_val)
    comparisons.append(f"{group1} vs {group2}")

# 본페로니는 Tukey 보다 보수적이다. 분산분석의 구조를 쓰지 않고
# 검정 수만으로 문턱을 낮추기 때문이다.
_, p_values_corrected, _, _ = multipletests(p_values, alpha=0.05, method='bonferroni')

# 보정 전과 뒤를 나란히 찍어 무엇이 달라지는지 본다.
print("Bonferroni-Corrected Pairwise Comparisons:")
for comparison, p_val, p_val_corr in zip(comparisons, p_values, p_values_corrected):
    print(f"{comparison}: p-value = {p_val:.4f}, Bonferroni-corrected p-value = {p_val_corr:.4f}")

Bonferroni-Corrected Pairwise Comparisons:
ctrl vs trt1: p-value = 0.2490, Bonferroni-corrected p-value = 0.7471
ctrl vs trt2: p-value = 0.0469, Bonferroni-corrected p-value = 0.1406
trt1 vs trt2: p-value = 0.0075, Bonferroni-corrected p-value = 0.0226

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])
n, k, nu = 10, 3, 27

print(f"{'k':>4}{'쌍 수 m':>9}{'1-(1-a)^m':>12}{'합집합 한계 ma':>14}")
for kk in range(2, 7):
    mm = kk * (kk - 1) // 2
    print(f"{kk:>4}{mm:>9}{1 - 0.95 ** mm:>12.4f}{min(1, 0.05 * mm):>14.4f}")

print(f"\n세 문턱 (|t| 척도, nu = {nu})")
print(f"  보정 없음  t(0.975, 27)        = {stats.t.ppf(0.975, nu):.4f}")
print(f"  투키       q(0.95,3,27)/sqrt2  = {stats.studentized_range.ppf(0.95, k, nu) / np.sqrt(2):.4f}")
print(f"  본페로니   t(1-0.05/6, 27)     = {stats.t.ppf(1 - 0.05 / (2 * 3), nu):.4f}")

g = df.groupby('group').weight
m, s = g.mean().values, g.std().values
MSE = (s ** 2).mean()
SE = np.sqrt(2 * MSE / n)
labels = ['ctrl', 'trt1', 'trt2']
print(f"\n같은 합동 MSE(자유도 27)로 세 방법의 p-값을 나란히")
print(f"{'쌍':>12}{'|t|':>8}{'보정없음':>10}{'투키':>9}{'본페로니':>10}")
for i in range(k):
    for j in range(i + 1, k):
        t = abs(m[j] - m[i]) / SE
        p_raw = 2 * stats.t.sf(t, nu)
        p_tuk = stats.studentized_range.sf(t * np.sqrt(2), k, nu)
        print(f"{labels[i] + '-' + labels[j]:>12}{t:>8.4f}{p_raw:>10.4f}"
              f"{p_tuk:>9.4f}{min(1, 3 * p_raw):>10.4f}")

print(f"\n쪽의 코드가 쓴 두 집단만의 t-검정 (자유도 18)")
for i in range(k):
    for j in range(i + 1, k):
        a = df[df.group == labels[i]].weight
        b = df[df.group == labels[j]].weight
        _, p = stats.ttest_ind(a, b)
        print(f"{labels[i] + '-' + labels[j]:>12}  p = {p:.4f},  x3 = {min(1, 3 * p):.4f}")

   k    쌍 수 m   1-(1-a)^m     합집합 한계 ma
   2        1      0.0500        0.0500
   3        3      0.1426        0.1500
   4        6      0.2649        0.3000
   5       10      0.4013        0.5000
   6       15      0.5367        0.7500

세 문턱 (|t| 척도, nu = 27)
  보정 없음  t(0.975, 27)        = 2.0518
  투키       q(0.95,3,27)/sqrt2  = 2.4794
  본페로니   t(1-0.05/6, 27)     = 2.5525

같은 합동 MSE(자유도 27)로 세 방법의 p-값을 나란히
           쌍     |t|      보정없음       투키      본페로니
   ctrl-trt1  1.3308    0.1944   0.3909    0.5832
   ctrl-trt2  1.7720    0.0877   0.1980    0.2630
   trt1-trt2  3.1028    0.0045   0.0120    0.0134

쪽의 코드가 쓴 두 집단만의 t-검정 (자유도 18)
   ctrl-trt1  p = 0.2490,  x3 = 0.7471
   ctrl-trt2  p = 0.0469,  x3 = 0.1406
   trt1-trt2  p = 0.0075,  x3 = 0.0226

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats
import pandas as pd

def load_data():
    """PlantGrowth 자료를 읽어 집단별로 나누고 자유도까지 함께 돌려준다."""
    url = 'https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/1dcc2bf5f955cc1224a3e1307256e1fe86b68dae/csv/datasets/PlantGrowth.csv'
    df = pd.read_csv(url, usecols=[1, 2])

    grouped_data = df.groupby('group')
    data_ctrl = grouped_data.get_group('ctrl').weight
    data_trt1 = grouped_data.get_group('trt1').weight
    data_trt2 = grouped_data.get_group('trt2').weight
    data = (data_ctrl, data_trt1, data_trt2)

    total_samples = data_ctrl.shape[0] + data_trt1.shape[0] + data_trt2.shape[0]
    num_groups = len(data)
    df1 = num_groups - 1
    df2 = total_samples - num_groups

    return df, data, df1, df2

def perform_anova(data_ctrl, data_trt1, data_trt2):
    """세 집단에 일원배치 분산분석을 수행한다."""
    statistic, p_value = stats.f_oneway(data_ctrl, data_trt1, data_trt2)
    print("\nOne-way ANOVA Results:")
    print(f"F-statistic = {statistic:.4f}")
    print(f"P-value = {p_value:.4f}\n")
    return statistic, p_value

def perform_tukey_hsd(data_ctrl, data_trt1, data_trt2, confidence_level=0.95):
    """Tukey HSD 사후비교를 수행하고 쌍별 신뢰구간을 보여 준다.

    구간이 0 을 품으면 그 쌍은 유의하지 않다고 읽는다. 신뢰수준을 높이면
    구간이 넓어지므로 유의하다고 판정되는 쌍이 줄어든다.
    """
    result = stats.tukey_hsd(data_ctrl, data_trt1, data_trt2)
    print(result)

    print(f"\nTukey's HSD Pairwise Group Comparisons ({confidence_level:.0%} Confidence Interval)")
    print("Comparison    Lower CI   Upper CI")
    confidence_interval = result.confidence_interval(confidence_level=confidence_level)
    for ((i, j), low) in np.ndenumerate(confidence_interval.low):
        if i < j:
            high = confidence_interval.high[i, j]
            print(f" ({i} - {j})   {low:>10.3f}   {high:>9.3f}")
    print()

# 자료 읽기 → 분산분석 → 사후비교 순으로 돌린다.
df, (data_ctrl, data_trt1, data_trt2), df1, df2 = load_data()

# 먼저 분산분석.
statistic, p_value = perform_anova(data_ctrl, data_trt1, data_trt2)

# 기본 95% 신뢰수준으로 사후비교.
perform_tukey_hsd(data_ctrl, data_trt1, data_trt2)

# 같은 자료를 99% 로 다시 본다. 구간이 넓어지는 만큼 결론이 보수적이 된다.
perform_tukey_hsd(data_ctrl, data_trt1, data_trt2, confidence_level=0.99)

One-way ANOVA Results:
F-statistic = 4.8461
P-value = 0.0159

Tukey's HSD Pairwise Group Comparisons (95.0% Confidence Interval)
Comparison  Statistic  p-value  Lower CI  Upper CI
 (0 - 1)      0.371     0.391    -0.320     1.062
 (0 - 2)     -0.494     0.198    -1.185     0.197
 (1 - 0)     -0.371     0.391    -1.062     0.320
 (1 - 2)     -0.865     0.012    -1.556    -0.174
 (2 - 0)      0.494     0.198    -0.197     1.185
 (2 - 1)      0.865     0.012     0.174     1.556


Tukey's HSD Pairwise Group Comparisons (95% Confidence Interval)
Comparison    Lower CI   Upper CI
 (0 - 1)       -0.320       1.062
 (0 - 2)       -1.185       0.197
 (1 - 2)       -1.556      -0.174

Tukey's HSD Pairwise Group Comparisons (95.0% Confidence Interval)
Comparison  Statistic  p-value  Lower CI  Upper CI
 (0 - 1)      0.371     0.391    -0.320     1.062
 (0 - 2)     -0.494     0.198    -1.185     0.197
 (1 - 0)     -0.371     0.391    -1.062     0.320
 (1 - 2)     -0.865     0.012    -1.556    -0.174
 (2 - 0)      0.494     0.198    -0.197     1.185
 (2 - 1)      0.865     0.012     0.174     1.556


Tukey's HSD Pairwise Group Comparisons (99% Confidence Interval)
Comparison    Lower CI   Upper CI
 (0 - 1)       -0.515       1.257
 (0 - 2)       -1.380       0.392
 (1 - 2)       -1.751       0.021

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])
n, k, nu = 10, 3, 27
g = df.groupby('group').weight
m, s = g.mean().values, g.std().values
MSE = (s ** 2).mean()
labels = ['ctrl', 'trt1', 'trt2']

print(f"{'신뢰수준':>8}{'q':>10}{'반폭':>10}")
for lev in [0.95, 0.99]:
    q = stats.studentized_range.ppf(lev, k, nu)
    print(f"{lev:>8.2f}{q:>10.6f}{q * np.sqrt(MSE / n):>10.6f}")

print(f"\n{'쌍':>12}{'차이':>9}{'95% 구간':>22}{'99% 구간':>22}{'p-adj':>9}")
for i in range(k):
    for j in range(i + 1, k):
        d = m[j] - m[i]
        t = abs(d) / np.sqrt(2 * MSE / n)
        p = stats.studentized_range.sf(t * np.sqrt(2), k, nu)
        h95 = stats.studentized_range.ppf(0.95, k, nu) * np.sqrt(MSE / n)
        h99 = stats.studentized_range.ppf(0.99, k, nu) * np.sqrt(MSE / n)
        print(f"{labels[i] + '-' + labels[j]:>12}{d:>9.3f}"
              f"{f'({d - h95:>7.3f}, {d + h95:>7.3f})':>22}"
              f"{f'({d - h99:>7.3f}, {d + h99:>7.3f})':>22}{p:>9.4f}")

print("\n쌍대성 확인: 구간이 0 을 품지 않는다  <=>  p-adj < alpha")
for lev, a in [(0.95, 0.05), (0.99, 0.01)]:
    h = stats.studentized_range.ppf(lev, k, nu) * np.sqrt(MSE / n)
    for i in range(k):
        for j in range(i + 1, k):
            d = m[j] - m[i]
            t = abs(d) / np.sqrt(2 * MSE / n)
            p = stats.studentized_range.sf(t * np.sqrt(2), k, nu)
            print(f"  alpha={a}  {labels[i] + '-' + labels[j]:>10}: "
                  f"0 제외={abs(d) > h},  p-adj<alpha={p < a}")

    신뢰수준         q        반폭
    0.95  3.506426  0.691216
    0.99  4.494842  0.886061

           쌍       차이                95% 구간                99% 구간    p-adj
   ctrl-trt1   -0.371    ( -1.062,   0.320)    ( -1.257,   0.515)   0.3909
   ctrl-trt2    0.494    ( -0.197,   1.185)    ( -0.392,   1.380)   0.1980
   trt1-trt2    0.865    (  0.174,   1.556)    ( -0.021,   1.751)   0.0120

쌍대성 확인: 구간이 0 을 품지 않는다  <=>  p-adj < alpha
  alpha=0.05   ctrl-trt1: 0 제외=False,  p-adj<alpha=False
  alpha=0.05   ctrl-trt2: 0 제외=False,  p-adj<alpha=False
  alpha=0.05   trt1-trt2: 0 제외=True,  p-adj<alpha=True
  alpha=0.01   ctrl-trt1: 0 제외=False,  p-adj<alpha=False
  alpha=0.01   ctrl-trt2: 0 제외=False,  p-adj<alpha=False
  alpha=0.01   trt1-trt2: 0 제외=False,  p-adj<alpha=False

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

# ToothGrowth 자료. 보충제 종류(supp)와 투여량(dose) 두 요인이 있다.
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)

                   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
import pandas as pd

url = 'https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/1dcc2bf5f955cc1224a3e1307256e1fe86b68dae/csv/datasets/ToothGrowth.csv'
df = pd.read_csv(url, usecols=[1, 2, 3])
a, b, n = 2, 3, 10                      # supp 2 수준, dose 3 수준, 칸마다 10 개

gm = df.len.mean()
mi = df.groupby('supp').len.mean()      # 행 평균
mj = df.groupby('dose').len.mean()      # 열 평균
cell = df.groupby(['supp', 'dose']).len.mean()

SSA = b * n * ((mi - gm) ** 2).sum()
SSB = a * n * ((mj - gm) ** 2).sum()
SScell = n * ((cell - gm) ** 2).sum()
SSAB = SScell - SSA - SSB
SSE = ((df.len - df.groupby(['supp', 'dose']).len.transform('mean')) ** 2).sum()
SST = ((df.len - gm) ** 2).sum()

print(f"전체평균 = {gm:.4f}")
print(f"행 평균 (supp) = {np.round(mi.values, 4)}")
print(f"열 평균 (dose) = {np.round(mj.values, 4)}")
print(f"칸 평균 =\n{cell.round(3)}")
print(f"\nSS_supp     = {SSA:.6f}   (df {a - 1})")
print(f"SS_dose     = {SSB:.6f}   (df {b - 1})")
print(f"SS_교호작용 = {SSAB:.6f}   (df {(a - 1) * (b - 1)})")
print(f"SS_잔차     = {SSE:.6f}   (df {a * b * (n - 1)})")
print(f"SS_전체     = {SST:.6f}   (df {a * b * n - 1})")
print(f"네 조각의 합 = {SSA + SSB + SSAB + SSE:.6f},  차 = {SST - SSA - SSB - SSAB - SSE:.2e}")

MSE = SSE / (a * b * (n - 1))
print(f"\nMSE = {MSE:.6f}")
print(f"F_supp = {(SSA / (a - 1)) / MSE:.6f}")
print(f"F_dose = {(SSB / (b - 1)) / MSE:.6f}")
print(f"F_교호 = {(SSAB / ((a - 1) * (b - 1))) / MSE:.6f}")

전체평균 = 18.8133
행 평균 (supp) = [20.6633 16.9633]
열 평균 (dose) = [10.605 19.735 26.1  ]
칸 평균 =
supp  dose
OJ    0.5     13.23
      1.0     22.70
      2.0     26.06
VC    0.5      7.98
      1.0     16.77
      2.0     26.14
Name: len, dtype: float64

SS_supp     = 205.350000   (df 1)
SS_dose     = 2426.434333   (df 2)
SS_교호작용 = 108.319000   (df 2)
SS_잔차     = 712.106000   (df 54)
SS_전체     = 3452.209333   (df 59)
네 조각의 합 = 3452.209333,  차 = 1.14e-12

MSE = 13.187148
F_supp = 15.571979
F_dose = 91.999965
F_교호 = 4.106991

from statsmodels.stats.multicomp import pairwise_tukeyhsd

# 주효과에 대한 사후비교. 교호작용이 유의하면 주효과를 이렇게 읽는 것이
# 오해를 부를 수 있어, 아래 단순효과 분석으로 넘어가는 편이 낫다.
tukey_dose = pairwise_tukeyhsd(endog=df['len'], groups=df['dose'], alpha=0.05)
print("Post-Hoc Test for Dose:")
print(tukey_dose)

# 보충제 종류에 대한 주효과 사후비교.
tukey_supp = pairwise_tukeyhsd(endog=df['len'], groups=df['supp'], alpha=0.05)
print("Post-Hoc Test for Supplement:")
print(tukey_supp)

Post-Hoc Test for Dose:
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
---------------------------------------------------
Post-Hoc Test for Supplement:
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
import pandas as pd
from scipy import stats

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

gm = df.len.mean()
mi = df.groupby('supp').len.mean()
SST = ((df.len - gm) ** 2).sum()
SSA = 3 * 10 * ((mi - gm) ** 2).sum()
SSE2 = ((df.len - df.groupby(['supp', 'dose']).len.transform('mean')) ** 2).sum()

MSE1 = (SST - SSA) / 58                    # supp 만 넣은 일원배치
MSE2 = SSE2 / 54                           # 이원배치 + 교호작용
print(f"일원배치 MSE = (SST - SS_supp)/58 = {MSE1:.6f}   (df 58)")
print(f"이원배치 MSE = SS_E/54            = {MSE2:.6f}   (df 54)")
print(f"비 = {MSE1 / MSE2:.4f},  제곱근 = {np.sqrt(MSE1 / MSE2):.4f}")

d = mi['VC'] - mi['OJ']
print(f"\nVC - OJ = {d:.4f}  (어느 쪽에서나 같다)")
for lab, mse, nu in [("일원배치", MSE1, 58), ("이원배치", MSE2, 54)]:
    se = np.sqrt(mse * (1 / 30 + 1 / 30))
    t = d / se
    print(f"  {lab}: SE = {se:.4f},  t = {t:.4f},  p = {2 * stats.t.sf(abs(t), nu):.6f}")

print(f"\n일원배치 MSE 안에 들어 있는 것")
SSB = 2 * 10 * ((df.groupby('dose').len.mean() - gm) ** 2).sum()
SSAB = SST - SSA - SSB - SSE2
print(f"  SS_dose   = {SSB:.4f}  ({SSB / (SST - SSA):.1%})")
print(f"  SS_교호   = {SSAB:.4f}  ({SSAB / (SST - SSA):.1%})")
print(f"  SS_잔차   = {SSE2:.4f}  ({SSE2 / (SST - SSA):.1%})")
print(f"  합        = {SSB + SSAB + SSE2:.4f}  = SST - SS_supp = {SST - SSA:.4f}")

일원배치 MSE = (SST - SS_supp)/58 = 55.980333   (df 58)
이원배치 MSE = SS_E/54            = 13.187148   (df 54)
비 = 4.2451,  제곱근 = 2.0604

VC - OJ = -3.7000  (어느 쪽에서나 같다)
  일원배치: SE = 1.9318,  t = -1.9153,  p = 0.060393
  이원배치: SE = 0.9376,  t = -3.9461,  p = 0.000231

일원배치 MSE 안에 들어 있는 것
  SS_dose   = 2426.4343  (74.7%)
  SS_교호   = 108.3190  (3.3%)
  SS_잔차   = 712.1060  (21.9%)
  합        = 3246.8593  = SST - SS_supp = 3246.8593

# 두 요인을 붙여 하나의 요인으로 만든다. 이러면 여섯 칸을 서로 견줄 수 있다.
df['supp_dose'] = df['supp'].astype(str) + "_" + df['dose'].astype(str)

# 칸 여섯 개의 모든 쌍을 견주므로 비교 횟수가 15 로 늘어난다. 그만큼 보수적이 된다.
tukey_interaction = pairwise_tukeyhsd(endog=df['len'], groups=df['supp_dose'], alpha=0.05)
print("Post-Hoc Test for Interaction (Supplement x Dose):")
print(tukey_interaction)

Post-Hoc Test for Interaction (Supplement x Dose):
 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
import pandas as pd
from scipy import stats

url = 'https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/1dcc2bf5f955cc1224a3e1307256e1fe86b68dae/csv/datasets/ToothGrowth.csv'
df = pd.read_csv(url, usecols=[1, 2, 3])
n, K, nu = 10, 6, 54
cell = df.groupby(['supp', 'dose']).len.mean()
MSE = ((df.len - df.groupby(['supp', 'dose']).len.transform('mean')) ** 2).sum() / nu
m = K * (K - 1) // 2
print(f"칸 수 K = {K},  쌍 수 m = {m},  보정 없는 FWER (독립 가정) = {1 - 0.95 ** m:.4f}")

q = stats.studentized_range.ppf(0.95, K, nu)
HSD = q * np.sqrt(MSE / n)
SE = np.sqrt(2 * MSE / n)
tb = stats.t.ppf(1 - 0.05 / (2 * m), nu)
print(f"\nMSE = {MSE:.6f},  SE = {SE:.6f}")
print(f"투키    q(0.95,6,54) = {q:.4f}  ->  |t| 문턱 {q / np.sqrt(2):.4f},  임계차 {HSD:.4f}")
print(f"본페로니 t(1-0.05/30, 54) = {tb:.4f}  ->  임계차 {tb * SE:.4f}")
print(f"투키가 낮은 폭: |t| 에서 {tb - q / np.sqrt(2):.4f},  임계차에서 {tb * SE - HSD:.4f}")

print(f"\n같은 용량끼리의 세 비교")
for d in [0.5, 1.0, 2.0]:
    diff = cell[('VC', d)] - cell[('OJ', d)]
    t = abs(diff) / SE
    p = stats.studentized_range.sf(t * np.sqrt(2), K, nu)
    print(f"  dose={d}: OJ-VC 차 = {-diff:>6.2f},  |t| = {t:.4f},  p-adj = {p:.4f},  "
          f"구간 ({diff - HSD:>7.3f}, {diff + HSD:>7.3f})")

칸 수 K = 6,  쌍 수 m = 15,  보정 없는 FWER (독립 가정) = 0.5367

MSE = 13.187148,  SE = 1.624017
투키    q(0.95,6,54) = 4.1783  ->  |t| 문턱 2.9545,  임계차 4.7981
본페로니 t(1-0.05/30, 54) = 3.0714  ->  임계차 4.9880
투키가 낮은 폭: |t| 에서 0.1169,  임계차에서 0.1899

같은 용량끼리의 세 비교
  dose=0.5: OJ-VC 차 =   5.25,  |t| = 3.2327,  p-adj = 0.0243,  구간 (-10.048,  -0.452)
  dose=1.0: OJ-VC 차 =   5.93,  |t| = 3.6514,  p-adj = 0.0074,  구간 (-10.728,  -1.132)
  dose=2.0: OJ-VC 차 =  -0.08,  |t| = 0.0493,  p-adj = 1.0000,  구간 ( -4.718,   4.878)

# 단순효과 분석: 보충제를 하나로 고정해 두고 투여량 효과만 본다.
# 교호작용이 있을 때 결과를 말이 되게 읽는 방법이다.
oj_data = df[df['supp'] == 'OJ']
vc_data = df[df['supp'] == 'VC']

# 보충제별로 따로 일원배치 분산분석을 돌린다.
oj_model = ols('len ~ C(dose)', data=oj_data).fit()
vc_model = ols('len ~ C(dose)', data=vc_data).fit()

# 두 결과를 견주면 교호작용이 무엇을 뜻하는지 드러난다.
print("ANOVA for Dose within Supplement OJ:")
print(anova_lm(oj_model))

print("ANOVA for Dose within Supplement VC:")
print(anova_lm(vc_model))

# 보충제별 사후비교.
print("Tukey HSD for Dose within Supplement OJ:")
print(pairwise_tukeyhsd(endog=oj_data['len'], groups=oj_data['dose'], alpha=0.05))

print("Tukey HSD for Dose within Supplement VC:")
print(pairwise_tukeyhsd(endog=vc_data['len'], groups=vc_data['dose'], alpha=0.05))

ANOVA for Dose within Supplement OJ:
            df      sum_sq     mean_sq          F        PR(>F)
C(dose)    2.0  885.264667  442.632333  31.441504  8.887164e-08
Residual  27.0  380.105000   14.077963        NaN           NaN
ANOVA for Dose within Supplement VC:
            df       sum_sq     mean_sq          F        PR(>F)
C(dose)    2.0  1649.488667  824.744333  67.072379  3.357317e-11
Residual  27.0   332.001000   12.296333        NaN           NaN
Tukey HSD for Dose within Supplement OJ:
Multiple Comparison of Means - Tukey HSD, FWER=0.05 
====================================================
group1 group2 meandiff p-adj   lower   upper  reject
----------------------------------------------------
   0.5    1.0     9.47    0.0  5.3096 13.6304   True
   0.5    2.0    12.83    0.0  8.6696 16.9904   True
   1.0    2.0     3.36 0.1309 -0.8004  7.5204  False
----------------------------------------------------
Tukey HSD for Dose within Supplement VC:
Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj  lower   upper  reject
---------------------------------------------------
   0.5    1.0     8.79   0.0  4.9018 12.6782   True
   0.5    2.0    18.16   0.0 14.2718 22.0482   True
   1.0    2.0     9.37   0.0  5.4818 13.2582   True
---------------------------------------------------

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

url = 'https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/1dcc2bf5f955cc1224a3e1307256e1fe86b68dae/csv/datasets/ToothGrowth.csv'
df = pd.read_csv(url, usecols=[1, 2, 3])
n = 10
gm = df.len.mean()
SSB = 2 * n * ((df.groupby('dose').len.mean() - gm) ** 2).sum()
SSE = ((df.len - df.groupby(['supp', 'dose']).len.transform('mean')) ** 2).sum()
SSA = 3 * n * ((df.groupby('supp').len.mean() - gm) ** 2).sum()
SST = ((df.len - gm) ** 2).sum()
SSAB = SST - SSA - SSB - SSE

tot_ss, tot_e = 0.0, 0.0
print(f"{'보충제':>6}{'SS_dose':>12}{'SSE':>12}{'MSE':>11}{'df':>5}{'Tukey 반폭':>12}")
for s in ['OJ', 'VC']:
    sub = df[df.supp == s]
    mu = sub.groupby('dose').len.mean()
    ss = n * ((mu - sub.len.mean()) ** 2).sum()
    sse = ((sub.len - sub.groupby('dose').len.transform('mean')) ** 2).sum()
    mse = sse / 27
    half = stats.studentized_range.ppf(0.95, 3, 27) * np.sqrt(mse / n)
    tot_ss += ss
    tot_e += sse
    print(f"{s:>6}{ss:>12.4f}{sse:>12.4f}{mse:>11.6f}{27:>5}{half:>12.4f}")

print(f"\n단순효과 SS 의 합   = {tot_ss:.4f}")
print(f"SS_dose + SS_교호   = {SSB:.4f} + {SSAB:.4f} = {SSB + SSAB:.4f}")
print(f"단순효과 SSE 의 합  = {tot_e:.4f}   (= 이원배치 SSE = {SSE:.4f})")

MSE2 = SSE / 54
half6 = stats.studentized_range.ppf(0.95, 6, 54) * np.sqrt(MSE2 / n)
print(f"\n여섯 칸 Tukey 의 반폭 = {half6:.4f}  (비교 15 개, df 54)")
print(f"단순효과 Tukey 의 반폭 = 위 표 (비교 3 개씩, df 27)")

print(f"\nOJ 안에서 dose 1.0 vs 2.0")
sub = df[df.supp == 'OJ']
mu = sub.groupby('dose').len.mean()
mse = ((sub.len - sub.groupby('dose').len.transform('mean')) ** 2).sum() / 27
d = mu[2.0] - mu[1.0]
t = d / np.sqrt(2 * mse / n)
print(f"  차 = {d:.2f},  |t| = {abs(t):.4f},  p-adj = {stats.studentized_range.sf(abs(t) * np.sqrt(2), 3, 27):.4f}")
sub = df[df.supp == 'VC']
mu = sub.groupby('dose').len.mean()
mse = ((sub.len - sub.groupby('dose').len.transform('mean')) ** 2).sum() / 27
d = mu[2.0] - mu[1.0]
t = d / np.sqrt(2 * mse / n)
print(f"VC 안에서 dose 1.0 vs 2.0")
print(f"  차 = {d:.2f},  |t| = {abs(t):.4f},  p-adj = {stats.studentized_range.sf(abs(t) * np.sqrt(2), 3, 27):.6f}")

   보충제     SS_dose         SSE        MSE   df    Tukey 반폭
    OJ    885.2647    380.1050  14.077963   27      4.1604
    VC   1649.4887    332.0010  12.296333   27      3.8882

단순효과 SS 의 합   = 2534.7533
SS_dose + SS_교호   = 2426.4343 + 108.3190 = 2534.7533
단순효과 SSE 의 합  = 712.1060   (= 이원배치 SSE = 712.1060)

여섯 칸 Tukey 의 반폭 = 4.7981  (비교 15 개, df 54)
단순효과 Tukey 의 반폭 = 위 표 (비교 3 개씩, df 27)

OJ 안에서 dose 1.0 vs 2.0
  차 = 3.36,  |t| = 2.0024,  p-adj = 0.1309
VC 안에서 dose 1.0 vs 2.0
  차 = 9.37,  |t| = 5.9750,  p-adj = 0.000007

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats
from statsmodels.stats.libqsturng import qsturng

rng = np.random.default_rng(13001)
M = 200_000
print("모의실험으로 만든 0.95 분위수와 qsturng 비교")
print(f"{'k':>3s} {'df':>4s} {'모의 0.95 분위수':>14s} {'qsturng':>9s} "
      f"{'√2·본페로니 z':>14s}")
for k, df in [(3, 12), (4, 20), (5, 30), (8, 40), (10, 60)]:
    z = rng.standard_normal((M, k))
    s = np.sqrt(rng.chisquare(df, M) / df)
    q = (z.max(1) - z.min(1)) / s
    zb = stats.norm.ppf(1 - 0.05 / (2 * (k * (k - 1) // 2)))
    print(f"{k:3d} {df:4d} {np.quantile(q, 0.95):14.4f} "
          f"{qsturng(0.95, k, df):9.4f} {np.sqrt(2) * zb:14.4f}")

print("\nq 임계값이 k 에 따라 얼마나 자라는가 (df=60 고정)")
print(f"{'k':>4s} {'비교 수':>7s} {'q(0.95)':>9s} {'q/√2':>8s} "
      f"{'본페로니 t':>10s} {'차이':>7s}")
for k in [2, 3, 4, 5, 6, 8, 10, 15, 20]:
    M2 = k * (k - 1) // 2
    q = qsturng(0.95, k, 60)
    bc = stats.t.ppf(1 - 0.05 / (2 * M2), 60)
    print(f"{k:4d} {M2:7d} {q:9.4f} {q / np.sqrt(2):8.4f} {bc:10.4f} "
          f"{100 * (bc - q / np.sqrt(2)) / (q / np.sqrt(2)):6.2f}%")

모의실험으로 만든 0.95 분위수와 qsturng 비교
  k   df    모의 0.95 분위수   qsturng      √2·본페로니 z
  3   12         3.7771    3.7711         3.3856
  4   20         3.9532    3.9585         3.7311
  5   30         4.1114    4.1020         3.9697
  8   40         4.5204    4.5206         4.4176
 10   60         4.6474    4.6461         4.6114

q 임계값이 k 에 따라 얼마나 자라는가 (df=60 고정)
   k    비교 수   q(0.95)     q/√2     본페로니 t      차이
   2       1    2.8288   2.0003     2.0003   0.00%
   3       3    3.3986   2.4032     2.4630   2.49%
   4       6    3.7372   2.6426     2.7286   3.25%
   5      10    3.9774   2.8124     2.9146   3.63%
   6      15    4.1630   2.9437     3.0573   3.86%
   8      28    4.4408   3.1401     3.2697   4.13%
  10      45    4.6461   3.2853     3.4260   4.28%
  15     105    5.0010   3.5362     3.6960   4.52%
  20     190    5.2412   3.7061     3.8789   4.66%

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from statsmodels.stats.libqsturng import qsturng

def tk_any(gs, alpha=0.05):
    k = len(gs)
    n = np.array([len(g) for g in gs], float)
    m = np.array([g.mean() for g in gs])
    v = np.array([g.var(ddof=1) for g in gs])
    MSE = ((n - 1) * v).sum() / (n.sum() - k)
    qc = qsturng(1 - alpha, k, n.sum() - k)
    for i in range(k):
        for j in range(i + 1, k):
            if abs(m[i] - m[j]) / np.sqrt(MSE / 2 * (1 / n[i] + 1 / n[j])) > qc:
                return True
    return False

rng = np.random.default_rng(13002)
B = 5_000
print("등분산·완전 귀무, 명목 0.05")
print(f"{'집단 크기':>26s} {'k':>3s} {'FWER':>8s}")
for ns in [[15, 15, 15], [5, 15, 25], [3, 10, 30],
           [20, 20, 20, 20], [5, 10, 20, 40], [4, 4, 4, 40]]:
    a = sum(tk_any([rng.normal(0, 1, n) for n in ns]) for _ in range(B))
    print(f"{str(ns):>26s} {len(ns):3d} {a / B:8.4f}")

등분산·완전 귀무, 명목 0.05
                     집단 크기   k     FWER
              [15, 15, 15]   3   0.0502
               [5, 15, 25]   3   0.0506
               [3, 10, 30]   3   0.0444
          [20, 20, 20, 20]   4   0.0514
           [5, 10, 20, 40]   4   0.0408
             [4, 4, 4, 40]   4   0.0420

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats
from statsmodels.stats.libqsturng import qsturng

def tk_any(gs, alpha=0.05):
    k = len(gs)
    n = np.array([len(g) for g in gs], float)
    m = np.array([g.mean() for g in gs])
    v = np.array([g.var(ddof=1) for g in gs])
    MSE = ((n - 1) * v).sum() / (n.sum() - k)
    qc = qsturng(1 - alpha, k, n.sum() - k)
    return any(abs(m[i] - m[j]) / np.sqrt(MSE / 2 * (1 / n[i] + 1 / n[j])) > qc
               for i in range(k) for j in range(i + 1, k))

def snk_any(gs, alpha=0.05):
    """뉴먼-쿨스: 순위 간격 p 에 따라 q(α, p, df) 를 쓴다."""
    k = len(gs)
    n = np.array([len(g) for g in gs], float)
    m = np.array([g.mean() for g in gs])
    v = np.array([g.var(ddof=1) for g in gs])
    MSE = ((n - 1) * v).sum() / (n.sum() - k)
    dfe = n.sum() - k
    o = np.argsort(m)
    for a in range(k):
        for b in range(a + 1, k):
            i, j = o[a], o[b]
            p = b - a + 1
            if (abs(m[i] - m[j]) / np.sqrt(MSE / 2 * (1 / n[i] + 1 / n[j]))
                    > qsturng(1 - alpha, p, dfe)):
                return True
    return False

def lsd_any(gs, alpha=0.05):
    k = len(gs)
    n = np.array([len(g) for g in gs], float)
    m = np.array([g.mean() for g in gs])
    v = np.array([g.var(ddof=1) for g in gs])
    MSE = ((n - 1) * v).sum() / (n.sum() - k)
    tc = stats.t.ppf(1 - alpha / 2, int(n.sum() - k))
    return any(abs(m[i] - m[j]) / np.sqrt(MSE * (1 / n[i] + 1 / n[j])) > tc
               for i in range(k) for j in range(i + 1, k))

rng = np.random.default_rng(13002)
B = 5_000
print("등분산·균형·완전 귀무에서의 FWER, 명목 0.05")
print(f"{'k':>3s} {'n':>4s} {'튜키':>8s} {'뉴먼-쿨스':>10s} {'LSD(무보정)':>12s}")
for k, n in [(3, 12), (4, 12), (6, 12), (8, 12)]:
    a = b = c = 0
    for _ in range(B):
        gs = [rng.normal(0, 1, n) for _ in range(k)]
        a += tk_any(gs)
        b += snk_any(gs)
        c += lsd_any(gs)
    print(f"{k:3d} {n:4d} {a / B:8.4f} {b / B:10.4f} {c / B:12.4f}")