如何在R的年度丰度均值柱状图中添加GLM（负二项）显著性标记

针对你的需求，以下是基于R的完整实现步骤，包含数据修正、模型检验和带显著性标记的可视化：

1. 修正模拟数据与依赖包加载 #

先修正原代码中Group赋值的类型错误，同时加载所需的可视化和检验工具包：

library(tidyverse)
library(glmmTMB)
library(emmeans) # 用于事后成对检验
library(ggsignif) # 用于添加显著性标记

set.seed(123) 
n <- 72
fake_years <- sample(2000:2020, size = 12, replace = FALSE)
fake_group <- tibble(
  Year = fake_years,
  Temp = round(runif(12, min = -10, max = 40), 1), 
  Group = sample(1:2, 12, replace = TRUE) 
)

fake_data <- fake_group %>%
  slice(rep(1:n(), each = 6)) %>% 
  mutate(Site = sample(paste0("Site", 1:5), n, replace = TRUE), 
         Totals = round(runif(n, min = 0, max = 1000), 1))

# 直接关联原Group数据，避免手动赋值错误
fake_bmir <- fake_data %>% 
  group_by(Year) %>%
  summarise(mean = mean(Totals),num_obs=n(),sum_year_totals=sum(Totals),
            sd_year_totals= sd(Totals),se_mean=sd_year_totals/sqrt(num_obs),
            se_upper=mean+se_mean,se_lower=mean-se_mean) %>%
  left_join(fake_group %>% select(Year, Group), by = "Year")

2. GLM模型与显著性检验 #

基于你的负二项GLM模型，用emmeans做事后成对检验，筛选出显著差异的年份对：

# 拟合原负二项GLM模型
NB1 <- glmmTMB(Totals~ Group+Temp, family=nbinom2, data=fake_data) 

# 控制Temp变量，做年份间的成对比较
year_em <- emmeans(NB1, ~ Year, covariate = Temp)
year_pairs <- pairs(year_em)
# 提取显著差异（p<0.05）的年份对
sig_year_pairs <- as.data.frame(year_pairs) %>%
  filter(p.value < 0.05) %>%
  mutate(
    # 拆分对比的两个年份
    year1 = str_extract(contrast, "^\\d+"),
    year2 = str_extract(contrast, "\\d+$"),
    # 转换为x轴位置（对应因子年份的索引）
    x1 = match(year1, levels(factor(fake_bmir$Year))),
    x2 = match(year2, levels(factor(fake_bmir$Year))),
    # 生成显著性标记
    sig_label = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01 ~ "**",
      p.value < 0.05 ~ "*"
    )
  )

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

用ggplot2结合ggsignif绘制年度均值柱状图，并添加显著年份对的标记：

# 计算标记的y轴基准位置（误差棒上方）
y_base <- max(fake_bmir$se_upper) + 40

ggplot(fake_bmir, aes(x = factor(Year), y = mean, fill = factor(Group))) +
  geom_bar(stat = "identity", position = position_dodge()) +
  geom_errorbar(aes(ymin = se_lower, ymax = se_upper), width = 0.2, position = position_dodge(0.9)) +
  # 添加手动指定的显著性标记
  geom_signif(
    data = sig_year_pairs,
    aes(xmin = x1, xmax = x2, annotations = sig_label, y_position = y_base),
    manual = TRUE,
    tip_length = 0.02, # 标记竖线长度
    text_size = 5
  ) +
  labs(title = "年度丰度均值对比（带标准误与显著性标记）", 
       x = "年份", y = "丰度均值", fill = "温度组") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # 旋转年份标签避免重叠

关键说明 #

• ggsignif是添加显著性标记的便捷工具，需提前安装：install.packages("ggsignif")
• emmeans用于控制协变量（如你的Temp）的事后检验，确保差异结果的准确性，安装：install.packages("emmeans")
• 若需展示冷暖组的整体差异，可将geom_signif的comparisons参数设为list(c("1", "2"))，直接标记组间差异

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