如何将emmeans Tukey检验的P值添加至ggplot箱线图？

我正在研究蓝莓不同杂交组合与校正果实质量（实际产量的替代指标）的关系，构建混合线性模型sm <- lmer(adj_fruitmass ~ pair + (1|plantid), data = pilot)后发现显著效应，希望通过emmeans执行Tukey事后检验找出最优杂交组合。但尝试将检验的P值（或显著性标记星号）添加到ggplot箱线图时，当前代码生成的对比括号无标记且仍显示大写"NS"，需要更优解决方法。

原尝试代码 #

sm <- lmer(adj_fruitmass ~ pair + (1|plantid), data = pilot)

# Post-hoc test 
cmpr <- list(c("Bluecrop", "Duke"),
             c("Bluecrop", "Reka"), 
             c("Duke", "Bluecrop"),
             c("Duke", "Reka"))

emmeans(sm, pairwise ~ pair, adjust = "tukey")

# Change one NA to zero to fix order of boxplots
pilot$adj_fruitmass[is.na(pilot$adj_fruitmass)] <- 0

# Add post-hoc p-values to boxplot
pilot %>%
ggplot(aes(x = reorder(Pollen_donor, adj_fruitmass), y = adj_fruitmass, na.rm = T)) +
  geom_boxplot() +
  geom_point() +
  labs(x = "Pollen donor", y = "Adjusted fruit mass") +
  theme_classic() +
  facet_grid(~Variety) +
  stat_compare_means(comparisons = cmpr, 
                     label = "p.signif", 
                     hide.ns = TRUE,
                     tip.length = 0.01,
                     symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), 
         symbols = c("****", "***", "**", "*", "ns")))

问题原因 #

• stat_compare_means默认基于常规t检验/方差分析，而非你用lmer拟合的混合线性模型的emmeans结果，导致检验结果不匹配，出现无标记或错误标记的情况。
• 对比列表cmpr存在重复项（如c("Bluecrop", "Duke")和c("Duke", "Bluecrop")），干扰可视化逻辑。
• 将NA替换为0会扭曲数据分布，破坏箱线图的准确性。

解决方案 #

直接基于emmeans的Tukey检验结果生成显著性标记，借助ggsignif包实现精准匹配：

步骤1：获取emmeans的Tukey检验结果 #

library(lme4)
library(emmeans)
library(ggplot2)
library(ggsignif)
library(dplyr)
library(tidytext) # 用于分面排序

# 拟合混合线性模型
sm <- lmer(adj_fruitmass ~ pair + (1|plantid), data = pilot)

# 执行Tukey事后检验并提取结构化结果
emm <- emmeans(sm, ~ pair)
tukey_results <- pairs(emm, adjust = "tukey") %>% 
  as.data.frame() %>% 
  # 生成标准化对比组名称
  rowwise() %>% 
  mutate(comparison = list(c(group1, group2))) %>% 
  ungroup()

步骤2：生成自定义显著性标记 #

# 定义显著性规则并生成标记
tukey_results <- tukey_results %>% 
  mutate(signif_label = case_when(
    p.value < 0.0001 ~ "****",
    p.value < 0.001 ~ "***",
    p.value < 0.01 ~ "**",
    p.value < 0.05 ~ "*",
    TRUE ~ "ns"
  ))

# 提取去重后的对比列表
clean_comparisons <- tukey_results$comparison

步骤3：绘制箱线图并添加匹配的显著性标记 #

# 移除NA而非替换，保证数据真实性
pilot %>%
  drop_na(adj_fruitmass) %>%
  ggplot(aes(x = reorder_within(Pollen_donor, adj_fruitmass, Variety), 
             y = adj_fruitmass)) +
  geom_boxplot() +
  geom_point(alpha = 0.6) +
  labs(x = "花粉供体", y = "校正果实质量") +
  theme_classic() +
  facet_grid(~Variety, scales = "free_x") +
  scale_x_reordered() + # 配合分面排序使用
  # 添加基于emmeans结果的显著性标记
  geom_signif(
    comparisons = clean_comparisons,
    # 根据数据范围调整标记垂直位置
    y_position = max(pilot$adj_fruitmass, na.rm = T) + seq(0.3, 1.2, by = 0.3),
    annotations = tukey_results$signif_label,
    tip_length = 0.01,
    hide.ns = TRUE # 隐藏无显著性的标记
  )

关键改进点 #

• 直接复用emmeans的检验结果，确保显著性标记与混合模型分析完全匹配。
• 移除NA而非替换为0，保留数据原始分布。
• 使用reorder_within实现分面下的x轴合理排序，避免混乱。
• 用geom_signif替代stat_compare_means，实现自定义化的显著性标记控制。

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