添加与 GLM 的显着差异？

Question

在 R 中，我正在查看一个数据集，其中计算了多年来一个月内的总丰度（总计）。每年从多个地点收集数据；根据平均温度，这些年份被分为两个不同的组（即温暖组或寒冷组）。 数据不正常，因此使用 GLM(negbin) 来确定组之间和年份温度范围内的显着性。 temp 的 p 值显着。如果绘制了显示每年平均丰度的图，您将如何将这些年份的显着差异添加到图中？

数据基本上看起来像这样：

library(tidyverse)
library(glmmTMB)

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))

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)
fake_bmir<-fake_bmir %>% 
  mutate(Group=if_else(Year=="2000"|Year=="2004"|Year=="2010"|Year=="2017",1,0))
#GLM
NB1<- glmmTMB(Totals~ Group+Temp,
              family=nbinom2, data=fake_data) 
#Basic plot
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)) +
  labs(title = "Mean Totals by Year with Standard Error", x = "Year", y = "Mean Totals", fill = "Group") +
  theme_minimal()

我不太明白如何将差异添加到绘图上的多个条形图上。我想要的输出是显示对它们有一定程度的显着差异的年份。

Answer 1

假设您想每年分别拟合相同的模型（

Totals ~ Group + Temp

），您可以这样开始：

p_val_info <- (fake_data
  ## split data by year
  |> split(fake_data$Year)
  ## fit model to each chunk
  |> map(\(d) glmmTMB(Totals~ Group+Temp, family=nbinom2, data = d))
  ## extract coefficients/p-vals/etc. for each model, aggregate to a single df
  |> map_dfr(\(m) broom.mixed::tidy(m, effect = "fixed"), .id = "Year")
  ## keep only relevant stuff
  |> select(Year, term, estimate, p.value)
  ## convert Year back to numeric in case we want to merge
  |> mutate(across(Year, as.numeric))
  ## we're not usually interested in the significance of the intercept term ...
  |> filter(term != "(Intercept)")
)

（使用

nest

/

unnest

并在列表列上进行操作会更惯用...）

这给我们留下了一个数据框，其中包含每年

Group

和

Temp

的 p 值（以及估计值）。由于您的样本数据集很小，所有值都是

NA

- 要么您每年都有多个

Group

和

Temp

值的观察，或（唉）我不明白这个问题（也就是说，如果年内测试不可能的话，我不知道你可以“显示出对他们有一定程度的显着差异的年份”......）。

一旦你做到了这一步，你就可以

• 通过

  group_by()

  和

  summarise()

  的某种组合进一步减少数据（例如，您想要每年效果的最小 p 值吗？）
• 将此信息与

  fake_bmir

  合并（通过

  full_join(..., by = "Year")

• 使用

  geom_label()

  或

  geom_text()

  ，其中 x 和 y 由年份和

  fake_bmir

中误差条的顶部确定

ggpubr

包有一些自动将p值添加到绘图中的机制（例如基于组的成对对比），但我怀疑它可以自动处理这样的情况......