如何按Day执行单因素ANOVA并在分面柱状图添加显著性括号

问题背景 #

已绘制按Day分面的DoublingTime均值柱状图，需要分别对Day 1和Day 28执行单因素ANOVA，仅关注Amoxicillin组0与60、70、80的组内比较；但当前使用aov(DoublingTime ~ Amoxicillin * Day)会输出所有组间组合的结果，不知道如何筛选目标比较结果并添加到图上。

数据结构 #

DTDataFixed <- structure(list(Amoxicillin = structure(c(1L, 1L, 1L, 2L, 2L, 
2L, 3L, 3L, 3L, 4L, 4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 
4L, 4L, 4L), levels = c("0", "60", "70", "80"), class = "factor"), 
    Day = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), levels = c("1", 
    "28"), class = "factor"), DoublingTime = c(42.27, 39.38, 
    34.83, 53.32, 42.01, 43.32, 119.51, 111.8, 121.6, 96.27, 
    101.93, 99.02, 67.96, 71.46, 106.64, 133.3, 223.6, 128.36, 
    182.41, 101.93, 187.34, 99.02, 64.78, 50.59)), class = "data.frame", row.names = c(NA, 
-24L))

解决方案 #

1. 按Day拆分，单独执行单因素ANOVA #

不需要用包含交互项的全局模型，直接按Day分组后分别做单因素ANOVA，得到每个Day下Amoxicillin组的整体差异结果：

library(tidyverse)
library(rstatix)
library(ggpubr)

# 按Day分组做ANOVA
anova_results <- DTDataFixed %>%
  group_by(Day) %>%
  anova_test(DoublingTime ~ Amoxicillin) %>%
  adjust_pvalue(method = "bonferroni") %>%
  add_significance()

# 查看结果
anova_results

2. 筛选目标两两比较（0 vs 60/70/80） #

用pairwise_t_test指定参照组为0，只保留需要的对比，同时按Day分组：

# 按Day分组，仅做0与其他组的两两t检验
pairwise_results <- DTDataFixed %>%
  group_by(Day) %>%
  pairwise_t_test(
    DoublingTime ~ Amoxicillin,
    ref.group = "0",  # 指定参照组为0
    p.adjust.method = "bonferroni"
  ) %>%
  add_significance() %>%
  # 整理标记的纵坐标位置
  mutate(
    y.position = case_when(
      Day == "1" ~ max(DTDataFixed$DoublingTime[DTDataFixed$Day == "1"]) + 10,
      Day == "28" ~ max(DTDataFixed$DoublingTime[DTDataFixed$Day == "28"]) + 15
    )
  )

# 查看筛选后的结果
pairwise_results

3. 绘制带显著性标记的分面柱状图 #

用ggpubr的stat_pvalue_manual把筛选后的比较结果添加到图中，配合分面参数实现每个Day下的标记：

# 计算均值和标准误（与原代码一致）
data_summary <- DTDataFixed %>%
  group_by(Day, Amoxicillin) %>%
  summarise(
    mean_doubling_time = mean(DoublingTime),
    se = sd(DoublingTime) / sqrt(n()),
    .groups = "drop"
  )

# 绘图并添加显著性标记
ggplot(data_summary, aes(x = Amoxicillin, y = mean_doubling_time, fill = Amoxicillin)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black") +
  geom_errorbar(
    aes(ymin = mean_doubling_time - se, ymax = mean_doubling_time + se),
    width = 0.2, position = position_dodge(width = 0.9)
  ) +
  # 添加显著性标记
  stat_pvalue_manual(
    pairwise_results,
    x = "Amoxicillin",
    y.position = "y.position",
    label = "p.signif",  # 显示*、**等显著性标记
    tip.length = 0.01,
    facet.by = "Day"  # 按Day分面匹配标记
  ) +
  facet_wrap(~ Day) +
  labs(title = "Doubling Time by Day and Treatment", x = "Amoxicillin", y = "Mean Doubling Time") +
  theme_bw() + 
  scale_fill_brewer(palette = "Greys") +
  # 扩展y轴范围，避免标记超出图外
  expand_limits(y = max(pairwise_results$y.position) + 5)

关键说明 #

• 拆分Day做独立的ANOVA，避免全局交互模型带来的多余结果；
• 通过ref.group参数精准筛选需要的对比组，不用处理无关的组间比较；
• stat_pvalue_manual支持分面匹配，能自动把每个Day的标记放到对应分面中；
• 可根据需求调整y.position的数值，让显著性标记的位置更美观。

内容的提问来源于stack exchange，提问作者RayzedByRobots