如何在facet_wrap分面geom_boxplot中添加两两比较显著性标记

问题背景 #

此前通过ggboxplot()搭配ggarrange()拼接独立子图，结合stat_pvalue_manual()已经实现带Kruskal-Wallis检验后Dunn两两比较显著性标记的箱线图，效果如下：

对应旧代码如下：

(plot1.1 <- ggboxplot(pamT1, x = "Temperature", y = "Yield") +
  stat_pvalue_manual(pwc2, hide.ns = TRUE))
(plot2.2 <- ggboxplot(pamT2, x = "Temperature", y = "Yield") +
  stat_pvalue_manual(pwc4, hide.ns = TRUE))
(plot3.3 <- ggboxplot(pamT3, x = "Temperature", y = "Yield") +
  stat_pvalue_manual(pwc6, hide.ns = TRUE))
(plot4.4 <- ggboxplot(pamT4, x = "Temperature", y = "Yield") +
  stat_pvalue_manual(pwc8, hide.ns = TRUE))

ggarrange(plot1.1, plot2.2, plot3.3, plot4.4, 
          labels = c("A", "B", "C", "D"),
          ncol = 2, nrow = 2)

现在需要改用原生geom_boxplot()搭配facet_wrap()生成分面箱线图，目标效果是在每个分面添加*、**类显著性标记，目标图效果如下：

当前使用的分面绘图代码如下，运行时无法正常加载显著性标记：

ggplot(pam2, aes(x=Temperature, y=Yield, fill=Temperature)) + 
  geom_boxplot(alpha=0.9) +
  facet_wrap(pam2$Time, scales="free_y", labeller = labeller(Time=labels)) + 
  theme(legend.position="none") +
  scale_fill_manual(values=c("dodgerblue4", "dodgerblue", "grey70", "firebrick1", "firebrick")) +
  labs(y="Photosynthetic yield", x=xlab, fill=expression("Temperature " ( degree*C))) +
  theme_classic() +
  theme(legend.position="bottom") +
  scale_y_continuous(labels = label_number(accuracy = 0.01))

失败原因 #

stat_pvalue_manual()在分面场景下，要求传入的显著性结果表必须包含和分面变量同名同值的列，才能将每组显著性比较结果匹配到对应分面。此前拆分数据集做检验得到的pwc2/pwc4/pwc6/pwc8四个结果表没有分面对应的Time列，函数无法识别结果对应的分面位置，自然无法正常渲染标记。

解决方案 #

第一步：整理符合分面要求的显著性结果表 #

你可以选择两种方式生成结果表：

• 手动拼接已有检验结果：给每个已有的Dunn检验结果表添加对应分面的Time取值，再合并为总表

  library(dplyr)
  # 注意替换等号后面的内容为每个子数据集对应的Time实际分组值
  pwc2 <- pwc2 %>% mutate(Time = "pamT1对应的Time分组值")
  pwc4 <- pwc4 %>% mutate(Time = "pamT2对应的Time分组值")
  pwc6 <- pwc6 %>% mutate(Time = "pamT3对应的Time分组值")
  pwc8 <- pwc8 %>% mutate(Time = "pamT4对应的Time分组值")
  # 合并所有结果，自动计算适配free_y坐标轴的标记位置
  pwc_all <- bind_rows(pwc2, pwc4, pwc6, pwc8) %>%
    add_xy_position(x = "Temperature", scales = "free_y")
  

• 直接对完整长表做分组检验（更推荐，避免手动匹配出错）：无需拆分数据集，直接按分面变量分组完成所有检验，输出结果自动带分面列

  library(rstatix)
  pwc_all <- pam2 %>%
    group_by(Time) %>%
    # 按Time分组做Dunn两两比较，p值校正方法可按需替换
    dunn_test(Yield ~ Temperature, p.adjust.method = "holm") %>%
    # 自动计算每个分面下显著性标记的横纵位置，适配free_y刻度
    add_xy_position(x = "Temperature", scales = "free_y")
  

第二步：在分面图中添加显著性标记 #

在原有ggplot代码中加入stat_pvalue_manual()层，传入整理好的总显著性结果表即可，修正后的完整代码如下：

ggplot(pam2, aes(x=Temperature, y=Yield, fill=Temperature)) + 
  geom_boxplot(alpha=0.9) +
  # 显著性标记层
  stat_pvalue_manual(
    pwc_all,
    hide.ns = TRUE, # 隐藏无统计学显著性的标记
    label = "p.adj.signif", # 用*、**作为显著性标记文本
    tip.length = 0.01, # 调整标记竖线的长度
    bracket.size = 0.5 # 调整标记横线的粗细
  ) +
  # 注意facet_wrap的写法不要用pam2$Time，直接写~Time即可，避免匹配异常
  facet_wrap(~Time, scales="free_y", labeller = labeller(Time=labels)) + 
  scale_fill_manual(values=c("dodgerblue4", "dodgerblue", "grey70", "firebrick1", "firebrick")) +
  labs(y="Photosynthetic yield", x=xlab, fill=expression("Temperature " ( degree*C))) +
  theme_classic() +
  theme(legend.position="bottom") +
  # 可通过expand参数适当扩大Y轴范围，避免显著性标记超出绘图边界
  scale_y_continuous(
    labels = label_number(accuracy = 0.01),
    expand = expansion(mult = c(0.05, 0.15))
  )

常见问题调整 #

• 如果标记位置错位：检查pwc_all$Time的取值是否和pam2$Time完全一致，不要存在大小写、空格、因子水平不匹配的问题
• 如果标记压到箱线图或者超出坐标轴：可以手动修改pwc_all$y.position列的数值调整标记高度，也可以调整scale_y_continuous()的expand参数扩大Y轴范围
• 如果需要显示p值而非星号：把stat_pvalue_manual()的label参数改为"p.adj"即可，也可以通过add_significance()函数自定义星号对应的p值阈值

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