基于TukeyHSD输出为ggplot分组条形图自动添加显著性字母

依赖包安装与加载 #

# 安装缺失包
if (!requireNamespace("tidyverse", quietly = TRUE)) install.packages("tidyverse")
if (!requireNamespace("emmeans", quietly = TRUE)) install.packages("emmeans")
if (!requireNamespace("multcompView", quietly = TRUE)) install.packages("multcompView")

# 加载包
library(tidyverse)
library(emmeans)
library(multcompView)

单物种核心实现代码 #

# 1. 按Zone分组独立做检验、提取显著性字母与绘图统计量
stat_result <- hogs.sample %>%
  # 按Zone拆分数据，保证跨Zone无比较
  nest_by(Zone) %>%
  mutate(
    # 替换your_response_col为你实际用于方差分析的原始观测值列名
    aov_mod = list(aov(your_response_col ~ Levelname, data = data)),
    # Tukey事后检验
    emm = list(emmeans(aov_mod, ~ Levelname, adjust = "tukey")),
    # 生成显著性字母，按组均值从高到低排序，alpha=0.05
    cld = list(cld(emm, alpha = 0.05, Letters = letters, decreasing = TRUE))
  ) %>%
  # 合并检验结果
  unnest(cld) %>%
  # 计算绘图用的Y轴均值与误差（标准误，如需标准差可自行修改计算逻辑）
  group_by(Zone, Levelname) %>%
  summarise(
    y_mean = exp(mean(hogs.fit)),
    y_se = exp(sd(hogs.fit)/sqrt(n())),
    .groups = "drop"
  ) %>%
  # 合并显著性字母，去除默认生成的空格
  left_join(stat_result %>% select(Zone, Levelname, .group), by = c("Zone", "Levelname")) %>%
  mutate(.group = str_squish(.group))

# 2. 绘制带自动显著性标记的分组条形图
ggplot(stat_result, aes(x = Zone, y = y_mean, fill = Levelname)) +
  geom_col(position = position_dodge(width = 0.9), width = 0.8) +
  # 误差棒与条形对齐
  geom_errorbar(
    aes(ymin = y_mean - y_se, ymax = y_mean + y_se),
    position = position_dodge(width = 0.9),
    width = 0.2
  ) +
  # 显著性字母自动定位在误差棒上方，与对应条形对齐
  geom_text(
    aes(label = .group, y = y_mean + y_se),
    position = position_dodge(width = 0.9),
    vjust = -0.3, # 可根据出图效果微调上下偏移
    size = 3.5
  ) +
  labs(x = "Zone", y = "exp(hogs.fit)", fill = "Levelname") +
  theme_bw() +
  theme(panel.grid = element_blank())

多物种批量适配版本 #

如果数据集包含species列区分不同物种，可直接用以下代码一次性完成所有物种的检验与绘图：

# 批量统计检验
stat_result_multi <- hogs.sample %>%
  # 按物种+Zone双重分组，完全独立处理每个物种的每个Zone
  nest_by(species, Zone) %>%
  mutate(
    aov_mod = list(aov(your_response_col ~ Levelname, data = data)),
    emm = list(emmeans(aov_mod, ~ Levelname, adjust = "tukey")),
    cld = list(cld(emm, alpha = 0.05, Letters = letters, decreasing = TRUE))
  ) %>%
  unnest(cld) %>%
  group_by(species, Zone, Levelname) %>%
  summarise(
    y_mean = exp(mean(hogs.fit)),
    y_se = exp(sd(hogs.fit)/sqrt(n())),
    .groups = "drop"
  ) %>%
  left_join(stat_result_multi %>% select(species, Zone, Levelname, .group), by = c("species", "Zone", "Levelname")) %>%
  mutate(.group = str_squish(.group))

# 批量分面绘图
ggplot(stat_result_multi, aes(x = Zone, y = y_mean, fill = Levelname)) +
  geom_col(position = position_dodge(width = 0.9), width = 0.8) +
  geom_errorbar(
    aes(ymin = y_mean - y_se, ymax = y_mean + y_se),
    position = position_dodge(width = 0.9),
    width = 0.2
  ) +
  geom_text(
    aes(label = .group, y = y_mean + y_se),
    position = position_dodge(width = 0.9),
    vjust = -0.3,
    size = 3
  ) +
  facet_wrap(~species, scales = "free_y") + # 每个物种Y轴独立适配
  labs(x = "Zone", y = "exp(hogs.fit)", fill = "Levelname") +
  theme_bw()

注意事项：

1. 请将代码中your_response_col替换为你实际用于方差分析的观测值列名
2. 如误差棒需用标准差，只需将y_se计算逻辑改为sd(hogs.fit)即可
3. 字母大小、偏移量可根据出图需求调整geom_text的size和vjust参数

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