12 Chapter 12: 多樣本推論 (multisample inference)

三組以上的比較：ANOVA、事後比較與不要亂開比較清單

12.1 本章學習目標

前面我們學過兩組平均值比較。可是醫學研究常常不只兩組：標準照護、飲食介入、藥物介入、合併介入；或低劑量、中劑量、高劑量；或三種病房、四種治療策略、五個年齡層。當組別變多，直接做很多次兩兩 t 檢定會讓第一型錯誤率一路升高，像是一直抽獎券，遲早抽到一張看起來很顯著的。

多樣本推論 (multisample inference) 的核心問題是：三組以上的母群體平均值是否有差異？最常見的方法是一因子變異數分析 (one-way analysis of variance, one-way ANOVA)。ANOVA 先回答「是否至少有一組不同」，再視需要進行事後比較 (post hoc comparison) 找出差異在哪裡。

讀完本章後，你應該能夠：

說明為什麼三組以上比較不宜只做多次 t 檢定。
解釋 one-way ANOVA 的研究問題、虛無假設與 F 統計量。
檢查 ANOVA 的基本假設，包括獨立性、常態性與變異數同質性。
使用 Levene test 評估變異數同質性。
使用事後比較與 Bonferroni correction 控制多重比較。
報告效果大小，例如 eta-squared。
認識二因子設計與交互作用 (interaction) 的基本直覺。

12.2 為什麼不要一直做 t 檢定？

假設有四組介入，每次兩兩比較共有 6 次 t 檢定。如果每次檢定都使用 \(\alpha = 0.05\)，整體出現至少一個偽陽性的機率會超過 0.05。這稱為多重比較問題 (multiple comparisons problem)。

ANOVA 的做法是先用一個整體檢定檢查所有組平均值是否完全相同。若整體檢定沒有顯著證據，通常不急著進行大量兩兩比較。若整體檢定顯著，再根據研究問題做事後比較，並使用校正方法控制錯誤率。

12.3 One-way ANOVA 的基本想法

One-way ANOVA 比較三組以上獨立樣本的平均值。虛無假設是：

\[ H_0: \mu_1 = \mu_2 = \cdots = \mu_k \]

對立假設則是至少有一組平均值不同。ANOVA 使用 F 統計量 (F statistic)，概念上比較「組間變異」與「組內變異」。如果組間差異相對於組內雜訊很大，F 值會變大，p 值會變小。

ANOVA 的常見假設包括：

各觀察值彼此獨立。
各組資料大致來自常態分布。
各組變異數大致相同，也稱變異數同質性 (homogeneity of variance)。

常態性在每組樣本數適中時通常有一定穩健性，但獨立性很重要，不能靠樣本數救回來。若同一位病人被重複測量，應考慮 repeated measures ANOVA 或混合模型，而不是把資料當成獨立樣本。

12.4 範例資料：四種介入後的收縮壓

假設高血壓門診比較四種 12 週介入後的收縮壓：標準照護、飲食介入、藥物介入、合併介入。每組 24 位病人。

# 寫死隨機生成的數值以確保範例數值不變
standard <- c(143.62770068055144, 144.35485428641238, 158.29722837368692, 125.93375277589783, 147.92646948365834, 145.91442545692033, 158.31017799196644, 135.6855988160743, 132.5608290574504, 155.7940012422921, 136.75734105068648, 152.425558864462, 145.5844440816137, 133.12419847502582, 150.69045124434243, 147.0876403521739, 151.6840161920338, 154.16901353442105, 159.3568664548177, 151.56918325702338, 129.83889117555742, 147.29108298898976, 131.13327881735768, 154.84031718221294)
diet <- c(145.33471710719624, 133.94370683756154, 155.8418981980534, 160.1112071226486, 138.9835817015546, 131.19894940513947, 150.28037166180346, 142.31581493306763, 149.3482227482019, 148.1347651463523, 129.5243402415302, 151.15097014011067, 139.30923865923018, 141.91543847730713, 132.26166514992303, 133.61936931930168, 136.3122065881174, 129.47943470938003, 131.76419591261308, 141.84872065240393, 148.4669567919367, 129.40305066316668, 151.08682907831908, 151.36485550649147)
drug <- c(144.2691583754587, 134.10008889928585, 141.12483564542896, 136.55070983350387, 143.7274164239949, 130.11867365194308, 124.98963283010454, 125.19077537188383, 147.11663446939514, 138.91715447563033, 116.08776675386767, 131.80315574119336, 132.28535731652957, 119.72013842683478, 136.71301554033374, 119.1819497326891, 126.86181847974814, 126.42518852959648, 120.82864887853506, 142.1871867515295, 112.22782235032372, 140.83276031505136, 132.19959657054372, 119.99829432535205)
combo <- c(133.43457901674338, 148.5771920830947, 114.15182115283538, 116.8997237560121, 120.3836580143441, 128.93691179249942, 119.41605955867131, 157.9409412362761, 127.14674004297783, 117.39415691263339, 124.7790860626084, 137.61325264307024, 123.659543428558, 148.97995579602878, 125.88769327271092, 113.85160264198318, 123.63753408322685, 123.62778955652638, 137.93790308030165, 129.29136140390045, 120.11168224135015, 143.20638962166552, 118.94113360622664, 129.8736783324609)

sbp_df <- data.frame(
  sbp = c(standard, diet, drug, combo),
  group = c(rep("標準照護", length(standard)), rep("飲食介入", length(diet)), rep("藥物介入", length(drug)), rep("合併介入", length(combo))),
  stringsAsFactors = FALSE
)
sbp_df$group <- factor(sbp_df$group, levels = c("標準照護", "飲食介入", "藥物介入", "合併介入"))

sbp_df %>%
  group_by(group) %>%
  summarise(
    count = n(),
    mean = mean(sbp),
    std = sd(sbp),
    .groups = "drop"
  ) %>%
  mutate(across(c(mean, std), ~ round(., 2)))

# A tibble: 4 × 4
  group    count  mean   std
  <fct>    <int> <dbl> <dbl>
1 標準照護    24  146.  9.96
2 飲食介入    24  142.  9.22
3 藥物介入    24  131.  9.86
4 合併介入    24  129. 11.8

p1 <- ggplot(sbp_df, aes(x = group, y = sbp, fill = group)) +
  geom_boxplot(alpha = 0.7, show.legend = FALSE) +
  geom_jitter(color = "#333333", alpha = 0.55, width = 0.15, height = 0) +
  scale_fill_brewer(palette = "Set2") +
  labs(
    title = "四組介入後收縮壓分布",
    x = "介入組別",
    y = "12 週收縮壓 (mmHg)"
  )

p1

Figure 12.1: 四組介入後收縮壓分布。箱型圖搭配原始資料點，可同時看位置、變異與離群值。

圖形比 p 值更早說話。這張圖讓我們先看見合併介入組平均收縮壓較低，但各組之間仍有重疊。ANOVA 會協助我們評估這樣的組間差異是否大到超過組內變異能合理解釋的程度。

12.5 執行 one-way ANOVA

model_anova <- lm(sbp ~ group, data = sbp_df)
anova_res <- anova(model_anova)

anova_test_df <- data.frame(
  quantity = c("F 統計量", "p 值"),
  value = c(anova_res[["F value"]][1], anova_res[["Pr(>F)"]][1]),
  stringsAsFactors = FALSE
)
anova_test_df$value <- round(anova_test_df$value, 4)

anova_test_df

  quantity   value
1 F 統計量 15.5121
2     p 值  0.0000

若 ANOVA p 值很小，表示至少有一組平均值不同。但它不會告訴你是哪一組與哪一組不同。這很像值班醫師說「病房有狀況」，但還沒告訴你是哪一床。

12.6 各組平均值與信賴區間

除了整體檢定，我們也應報告各組平均值與信賴區間。這比只報 p 值更能呈現臨床差異大小。

mean_ci <- function(values) {
  n <- length(values)
  m_val <- mean(values)
  se <- sd(values) / sqrt(n)
  margin <- qt(0.975, df = n - 1) * se
  return(c(mean = m_val, lower = m_val - margin, upper = m_val + margin))
}

summary_rows <- list()
idx <- 1
for (g in unique(sbp_df$group)) {
  sub_vals <- sbp_df$sbp[sbp_df$group == g]
  ci_vals <- mean_ci(sub_vals)
  summary_rows[[idx]] <- data.frame(
    group = g,
    mean = ci_vals["mean"],
    lower = ci_vals["lower"],
    upper = ci_vals["upper"],
    stringsAsFactors = FALSE
  )
  idx <- idx + 1
}
summary_df <- do.call(rbind, summary_rows)
summary_df$group <- factor(summary_df$group, levels = c("標準照護", "飲食介入", "藥物介入", "合併介入"))
summary_df <- summary_df %>% arrange(group)

summary_df %>% mutate(across(c(mean, lower, upper), ~ round(., 2)))

         group   mean  lower  upper
mean  標準照護 145.58 141.37 149.79
mean1 飲食介入 141.79 137.90 145.69
mean2 藥物介入 130.98 126.81 135.14
mean3 合併介入 128.57 123.60 133.54

p2 <- ggplot(summary_df, aes(x = group, y = mean)) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1, color = "#2a9d8f", linewidth = 1) +
  geom_point(color = "#2a9d8f", size = 3) +
  labs(
    title = "各組平均值與信賴區間",
    x = "介入組別",
    y = "平均收縮壓與 95% 信賴區間"
  )

p2

Figure 12.2: 各組平均收縮壓與 95% 信賴區間。信賴區間用來描述估計不確定性，不是正式事後比較。

12.7 變異數同質性與 Levene test

ANOVA 假設各組變異數相近。Levene test 可用來檢查變異數是否有明顯差異。它的虛無假設是各組變異數相同。

levene_test_manual <- function(df, val_col, grp_col) {
  medians <- tapply(df[[val_col]], df[[grp_col]], median)
  df$z <- abs(df[[val_col]] - medians[df[[grp_col]]])
  fit <- lm(z ~ group, data = df)
  an <- anova(fit)
  return(list(statistic = an[["F value"]][1], p.value = an[["Pr(>F)"]][1]))
}

levene_result <- levene_test_manual(sbp_df, "sbp", "group")

levene_test_df <- data.frame(
  quantity = c("Levene 統計量", "p 值"),
  value = c(levene_result$statistic, levene_result$p.value),
  stringsAsFactors = FALSE
)
levene_test_df$value <- round(levene_test_df$value, 4)

levene_test_df

       quantity  value
1 Levene 統計量 0.1264
2          p 值 0.9443

若 Levene test p 值很小，代表變異數同質性可能不成立。此時可考慮 Welch ANOVA、資料轉換、穩健方法，或改用無母數方法如 Kruskal-Wallis。方法選擇要回到資料型態與研究問題，不是看到一個檢定顯著就立刻換整套人生。

12.8 效果大小：eta-squared

p 值受樣本數影響很大。多樣本比較最好同時報告效果大小。ANOVA 常用 eta-squared (\(\eta^2\))，表示總變異中有多少比例可由組別解釋。

eta_squared <- function(df, val_col, grp_col) {
  fit <- lm(as.formula(paste(val_col, "~", grp_col)), data = df)
  an <- anova(fit)
  ss_between <- an[["Sum Sq"]][1]
  ss_total <- sum(an[["Sum Sq"]])
  return(ss_between / ss_total)
}

eta2 <- eta_squared(sbp_df, "sbp", "group")

eta_df <- data.frame(
  quantity = c("eta-squared"),
  value = c(eta2),
  stringsAsFactors = FALSE
)
eta_df$value <- round(eta_df$value, 4)

eta_df

     quantity  value
1 eta-squared 0.3359

若 \(\eta^2 = 0.25\)，可粗略解讀為組別解釋了結果總變異的 25%。這不是因果證明，但能幫助判斷組別差異是否有實質意義。

12.9 事後比較與 Bonferroni correction

ANOVA 顯著後，常見下一步是比較哪些組別不同。若有 4 組，兩兩比較共有 6 次。Bonferroni correction 的想法很直接：把每個 p 值乘以比較次數，或把顯著水準除以比較次數。

Bonferroni 方法保守，但容易理解。保守的意思是比較不容易宣稱顯著，代價是可能降低檢定力。

pairwise_ttests <- function(df, value_col, group_col) {
  groups <- unique(df[[group_col]])
  n_groups <- length(groups)
  n_tests <- n_groups * (n_groups - 1) / 2
  
  rows <- list()
  idx <- 1
  for (i in 1:(n_groups-1)) {
    for (j in (i+1):n_groups) {
      g_a <- groups[i]
      g_b <- groups[j]
      vals_a <- df[[value_col]][df[[group_col]] == g_a]
      vals_b <- df[[value_col]][df[[group_col]] == g_b]
      t_res <- t.test(vals_a, vals_b, var.equal = TRUE)
      mean_diff <- mean(vals_a) - mean(vals_b)
      
      rows[[idx]] <- data.frame(
        comparison = paste(g_a, "-", g_b),
        mean_difference = mean_diff,
        raw_p = t_res$p.value,
        bonferroni_p = min(t_res$p.value * n_tests, 1.0),
        stringsAsFactors = FALSE
      )
      idx <- idx + 1
    }
  }
  return(do.call(rbind, rows))
}

posthoc_df <- pairwise_ttests(sbp_df, "sbp", "group")
posthoc_df %>% mutate(across(c(mean_difference, raw_p, bonferroni_p), ~ round(., 4)))

           comparison mean_difference  raw_p bonferroni_p
1 標準照護 - 飲食介入          3.7899 0.1781       1.0000
2 標準照護 - 藥物介入         14.6041 0.0000       0.0000
3 標準照護 - 合併介入         17.0115 0.0000       0.0000
4 飲食介入 - 藥物介入         10.8143 0.0003       0.0017
5 飲食介入 - 合併介入         13.2217 0.0001       0.0005
6 藥物介入 - 合併介入          2.4074 0.4464       1.0000

posthoc_plot <- posthoc_df %>% arrange(bonferroni_p)
posthoc_plot$comparison <- factor(posthoc_plot$comparison, levels = rev(posthoc_plot$comparison))

p3 <- ggplot(posthoc_plot, aes(x = bonferroni_p, y = comparison)) +
  geom_col(fill = "#457b9d", width = 0.6) +
  geom_vline(xintercept = 0.05, color = "#d62828", linetype = "dashed", linewidth = 1) +
  labs(
    title = "事後比較：多重比較校正後 p 值",
    x = "Bonferroni 校正後 p 值",
    y = "組間比較"
  )

p3

Figure 12.3: 事後兩兩比較的 Bonferroni 校正後 p 值。紅線為 0.05。

事後比較要與研究問題相配。若研究一開始只關心「合併介入是否優於標準照護」，就不一定需要報告所有可能兩兩比較。統計表格不是越長越有學問，有時只是讓讀者在 p 值森林裡迷路。

12.10 二因子設計與交互作用

有些研究不只一個因子。例如同時考慮治療組別與性別，或劑量與年齡層。二因子 ANOVA (two-way ANOVA) 可同時評估兩個因子的主效應 (main effect)，以及兩者是否有交互作用。

交互作用 (interaction) 的意思是：某個因子的效果會因另一個因子的水準而改變。例如合併介入對男性的降壓效果可能比女性更明顯。這不是說哪一性別「比較重要」，而是治療效果可能不完全一致。

# 寫死隨機生成的數值以確保範例數值不變
sex_vec <- c(rep("女性", 54), rep("男性", 54))
treatment_vec <- c(rep("標準照護", 18), rep("藥物介入", 18), rep("合併介入", 18), rep("標準照護", 18), rep("藥物介入", 18), rep("合併介入", 18))
sbp_vec <- c(145.62289188027577, 135.25639563946706, 145.21215059044005, 146.3849198545578, 144.27569953128852, 154.385072345434, 131.17300107621958, 134.24180253476123, 141.7277577158328, 131.39708774512732, 148.67449487182196, 138.2634747167062, 137.6580083691114, 142.01101341377122, 132.89463679954503, 128.1892546421823, 131.16011626172372, 141.33326345713908, 125.43426299936148, 138.46395303773355, 136.54016117968754, 144.64754516901468, 133.13187005321933, 116.58270352700106, 122.10618876947494, 145.1899716432932, 144.44216782127432, 129.7775676235261, 123.39023009836049, 129.40153164094357, 127.1698743184818, 135.65241443135602, 126.7156572525416, 126.67224962344373, 130.3785186817352, 143.15404431447814, 146.93934515950912, 141.2567163074632, 116.82603803630637, 123.70243902217464, 143.17072920733455, 133.80071341979283, 141.64330898583384, 114.21374140897233, 140.3320271781515, 123.91490797120757, 136.780088170077, 121.51917386023342, 121.85505881205869, 130.1080134826732, 131.17620392089418, 133.84556354461807, 126.7009826436455, 116.78175367698181, 144.61561198668448, 153.92940519546204, 147.12782246219237, 141.495647497892, 132.64129543960766, 161.5702923686874, 143.8424888899133, 160.48188560145095, 145.05541753846913, 144.55702529366081, 148.83609998807333, 145.53082270096522, 153.7423290474123, 151.0840913756423, 133.44667795080332, 158.2560904088788, 134.7778245198447, 148.07914187883827, 129.33714501814748, 137.06407298831772, 127.8987221749782, 138.7249194472063, 141.99691684336938, 127.04985120103173, 138.83390300472774, 146.00372257375025, 146.24483693282735, 122.62284361649904, 142.6069669100332, 134.60067562515448, 143.56387594719263, 138.1954621104294, 135.40071477997134, 135.39896378007305, 133.12510104607946, 135.22188616713674, 133.801977574994, 122.4868943430265, 130.7806701271456, 131.0880869019809, 116.43539877514914, 133.89433959233392, 131.592914568407, 132.0413786011036, 124.78792227175232, 126.46728223168563, 112.1662614155949, 124.90702219254028, 126.02291109088692, 127.05828949868916, 116.22556087422495, 128.13977633181682, 141.75461949458614, 117.97023736560335)

interaction_df <- data.frame(
  sex = sex_vec,
  treatment = treatment_vec,
  sbp = sbp_vec,
  stringsAsFactors = FALSE
)

head(interaction_df)

   sex treatment      sbp
1 女性  標準照護 145.6229
2 女性  標準照護 135.2564
3 女性  標準照護 145.2122
4 女性  標準照護 146.3849
5 女性  標準照護 144.2757
6 女性  標準照護 154.3851

interaction_summary <- interaction_df %>%
  group_by(treatment, sex) %>%
  summarise(
    mean_sbp = mean(sbp),
    se_sbp = sd(sbp) / sqrt(n()),
    lower = mean_sbp - qt(0.975, df = n() - 1) * se_sbp,
    upper = mean_sbp + qt(0.975, df = n() - 1) * se_sbp,
    .groups = "drop"
  )

interaction_summary$treatment <- factor(interaction_summary$treatment, levels = c("標準照護", "藥物介入", "合併介入"))
interaction_summary$sex <- factor(interaction_summary$sex, levels = c("女性", "男性"))

p4 <- ggplot(interaction_summary, aes(x = treatment, y = mean_sbp, color = sex, group = sex)) +
  geom_line(position = position_dodge(width = 0.15), linewidth = 1.2) +
  geom_errorbar(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.15), width = 0.1, linewidth = 1) +
  geom_point(position = position_dodge(width = 0.15), size = 3) +
  labs(
    title = "二因子情境：治療效果可能因性別而不同",
    x = "介入組別",
    y = "平均收縮壓 (mmHg)",
    color = "性別"
  )

p4

Figure 12.4: 治療組別與性別的交互作用圖。若線條不平行，提示治療效果可能因性別不同而改變。

交互作用圖不是正式檢定，但非常有用。若線條接近平行，代表各組差異在不同性別中大致相似；若線條明顯不平行，表示可能存在交互作用，需要進一步用適當模型檢定。

12.11 何時改用無母數方法？

若資料明顯偏態、離群值很多、樣本數小，或結果是有序尺度，Kruskal-Wallis test 可能比 one-way ANOVA 更適合。這在 Chapter 9 已介紹過。

不過，Kruskal-Wallis 也不是單純「中位數 ANOVA」。它檢定的是各組排序分布是否相同。若分布形狀差很多，解釋仍要小心。選方法時要看研究問題、資料型態與分布，而不是把「常態不常態」當成唯一開關。

12.12 常見錯誤

第一個錯誤是三組以上直接做一堆 t 檢定，卻沒有多重比較校正。這會增加偽陽性風險。

第二個錯誤是 ANOVA 顯著後就說「所有組都不同」。ANOVA 只告訴我們至少有一組不同，不能直接指出是哪幾組。

第三個錯誤是只報告 p 值，不報告各組平均值、信賴區間與效果大小。沒有這些資訊，讀者很難判斷臨床意義。

第四個錯誤是忽略研究設計。獨立樣本 ANOVA 不適用於同一批人重複測量的資料。

第五個錯誤是看到交互作用圖線條交叉就立刻過度解讀。交互作用需要模型與檢定支持，也需要臨床合理性。

12.13 本章重點整理

多樣本推論處理三組以上的比較。One-way ANOVA 先檢定所有組平均值是否相同，若整體檢定顯著，再依研究問題進行事後比較。多重比較會增加第一型錯誤率，因此需要 Bonferroni、Tukey 或其他校正方法。報告 ANOVA 時，應同時呈現各組摘要、整體檢定、效果大小與事後比較。

當研究設計有兩個因子時，除了各因子的主效應，也要留意交互作用。臨床資料常常不只問「哪個治療比較好」，還會問「對誰比較好」。這時候模型就要跟著問題長大一點。

12.14 小練習

某研究比較三種止痛藥的疼痛分數下降量。請說明 one-way ANOVA 的虛無假設。
若 ANOVA p 值為 0.01，是否表示三組兩兩都不同？為什麼？
四組比較共有幾個兩兩比較？若使用 Bonferroni correction，單次比較的顯著水準是多少？
請解釋 eta-squared 的意義。
如果 Levene test 顯著，你會考慮哪些替代方法？
請畫出一個可能有交互作用的臨床情境，並用文字描述線條為何不平行。

12.15 Glossary

中文術語	English term	說明
多樣本推論	multisample inference	處理三組以上資料比較的推論方法。
變異數分析	analysis of variance, ANOVA	比較多組平均值是否相同的統計方法。
單因子變異數分析	one-way ANOVA	只有一個分組因子的 ANOVA。
F 統計量	F statistic	ANOVA 中用於比較組間變異與組內變異的統計量。
組間變異	between-group variation	不同組平均值之間的變異。
組內變異	within-group variation	同一組內個體之間的變異。
變異數同質性	homogeneity of variance	各組變異數大致相同的假設。
Levene 檢定	Levene test	檢定各組變異數是否相同的方法。
事後比較	post hoc comparison	整體檢定後進行 of 組間比較。
多重比較問題	multiple comparisons problem	多次檢定增加偽陽性機率的問題。
Bonferroni 校正	Bonferroni correction	以校正顯著水準或 p 值控制多重比較的保守方法。
Eta-squared	eta-squared	ANOVA 中表示組別解釋總變異比例的效果大小。
主效應	main effect	某一因子平均而言對結果的影響。
交互作用	interaction	某一因子的效果會因另一因子的水準而改變。