ANOVA 후 "어느 집단이 다른가"

앞 장에서 PlantGrowth 데이터로 ANOVA를 수행해 p = .016을 얻었습니다. "세 집단 중 적어도 한 쌍이 다르다"는 것을 알았지만, ctrl vs trt1이 다른지, ctrl vs trt2가 다른지, trt1 vs trt2가 다른지는 모릅니다.

이것을 알아내는 것이 사후검정(post-hoc test)입니다.

왜 다중비교 보정이 필요한가

세 집단이면 비교 쌍이 3개(A-B, A-C, B-C)입니다. k개 집단이면 k(k-1)/2개의 쌍을 비교해야 합니다.

각 검정을 α = 0.05로 하면, m번 검정했을 때 가족별 오류율(familywise error rate, FWER)이 높아집니다.

$FWER = 1 - (1 - \alpha)^m$

20번 비교하면 실제 효과가 없는데도 1종 오류가 발생할 확률이 64%까지 올라갑니다. 다중비교 보정은 이 누적 오류를 제어합니다.

TukeyHSD() — 가장 많이 쓰는 사후검정

Tukey의 HSD(Honestly Significant Difference) 검정은 모든 쌍 비교를 수행하면서 FWER을 α로 통제합니다.

data(PlantGrowth)
aov_result <- aov(weight ~ group, data = PlantGrowth)

tukey_result <- TukeyHSD(aov_result)
print(tukey_result)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = weight ~ group, data = PlantGrowth)

$group
            diff        lwr        upr     p adj
trt1-ctrl -0.371 -1.0622161  0.3202161 0.3921943
trt2-ctrl  0.494 -0.1972161  1.1852161 0.1981290
trt2-trt1  0.865  0.1737839  1.5562161 0.0120064

결과를 해석합니다.

• trt1-ctrl: trt1이 ctrl보다 0.371 작지만 유의하지 않습니다 (p = .392).
• trt2-ctrl: trt2가 ctrl보다 0.494 크지만 유의하지 않습니다 (p = .198).
• trt2-trt1: trt2가 trt1보다 0.865 크고 유의합니다 (p = .012, 95% CI [0.17, 1.56]).

ANOVA는 유의했지만, 사후검정에서 ctrl vs trt 비교는 유의하지 않고 trt2 vs trt1만 유의합니다.

TukeyHSD 시각화

plot(tukey_result, las = 1)

또는 ggplot2로 더 보기 좋게 만들 수 있습니다.

library(ggplot2)

tukey_df <- as.data.frame(tukey_result$group)
tukey_df$comparison <- rownames(tukey_df)

ggplot(tukey_df, aes(x = comparison, y = diff,
                      ymin = lwr, ymax = upr,
                      color = `p adj` < 0.05)) +
  geom_pointrange(size = 0.8) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  scale_color_manual(values = c("FALSE" = "gray60", "TRUE" = "tomato"),
                     labels = c("유의하지 않음", "유의함 (p < .05)")) +
  labs(
    title = "TukeyHSD 사후검정 결과",
    x = "비교 쌍",
    y = "평균 차이 (95% 신뢰구간)",
    color = "유의성"
  ) +
  coord_flip() +
  theme_minimal()

0을 포함하지 않는 신뢰구간(빨간 점)이 유의한 차이를 나타냅니다.

Bonferroni 보정

Bonferroni 보정은 유의수준을 비교 횟수로 나누어 각 검정에 적용합니다.

$\alpha_{adjusted} = \frac{\alpha}{m}$

3번 비교하면 각 검정의 유의수준이 0.05/3 = 0.0167이 됩니다. 매우 보수적(conservative)인 방법입니다. 모든 집단 쌍 비교가 아닌 특정 비교를 계획하고 있을 때(planned comparison) 사용합니다.

# p.adjust()로 Bonferroni 보정
p_values <- c(0.03, 0.04, 0.012)  # 세 쌍의 t-검정 p-값
p.adjust(p_values, method = "bonferroni")

[1] 0.090 0.120 0.036

원래 모두 0.05 미만이었지만, 보정 후 세 번째만 유의합니다.

FDR 보정 — Benjamini-Hochberg (BH)

Bonferroni는 너무 보수적이어서 실제 효과를 놓칠 수 있습니다. 탐색적 연구에서는 False Discovery Rate(FDR) 보정을 사용합니다.

FDR은 FWER 대신 "유의하다고 선언된 결과 중 거짓양성의 비율"을 제어합니다. Bonferroni보다 통계적 검정력이 높습니다.

# pairwise.t.test()로 다중비교 + BH 보정
pairwise.t.test(PlantGrowth$weight, PlantGrowth$group,
                p.adjust.method = "BH")

	Pairwise comparisons using t tests with pooled SD 

data:  PlantGrowth$weight and PlantGrowth$group 

     ctrl  trt1 
trt1 0.194 -    
trt2 0.194 0.012

BH 보정 후에도 trt2 vs trt1이 유의합니다(p = .012).

방법별 p.adjust() 비교

# 같은 원시 p-값에 대한 다양한 보정 비교
raw_p <- c(0.001, 0.008, 0.039, 0.047, 0.082, 0.29)

data.frame(
  원시_p값     = raw_p,
  Bonferroni   = round(p.adjust(raw_p, method = "bonferroni"), 4),
  BH           = round(p.adjust(raw_p, method = "BH"), 4),
  Holm         = round(p.adjust(raw_p, method = "holm"), 4)
)

  원시_p값 Bonferroni     BH   Holm
1    0.001     0.0060 0.0060 0.0060
2    0.008     0.0480 0.0240 0.0400
3    0.039     0.2340 0.0780 0.1560
4    0.047     0.2820 0.0705 0.1410
5    0.082     0.4920 0.0984 0.2460
6    0.290     1.0000 0.2900 0.2900

Bonferroni가 가장 엄격(보수적)이고, BH가 가장 관대합니다. 연구 목적에 따라 적절한 방법을 선택합니다.

방법 선택 가이드

일반적으로 확증적 연구(가설 검증)에는 TukeyHSD나 Bonferroni를, 탐색적 연구(발견 지향)에는 BH를 사용합니다.