负二项模型中emmeans事后检验报错问题求助

问题根源 #

你使用MASS::glm.nb()拟合的模型包含随机效应项(1 | replic)，但该函数仅支持固定效应负二项模型，无法识别lme4风格的随机效应语法。这导致emmeans无法正确解析模型结构，进而抛出“未定义列”的错误。

正确建模与解决步骤 #

1. 替换建模工具 #

改用glmmTMB包拟合带随机效应的负二项模型，它对混合效应、零膨胀模型的支持更完善，且与emmeans兼容性更好：

# 安装所需包（首次使用时运行）
install.packages(c("glmmTMB", "emmeans", "car"))

# 加载包
library(glmmTMB)
library(emmeans)
library(car)

2. 重新拟合模型 #

# 拟合带随机效应的负二项模型
negbin_glmm <- glmmTMB(flwrs ~ tre*wa + (1 | replic), 
                       family = nbinom2,  # 选择负二项参数化方式，nbinom1/nbinom2二选一
                       data = test)

# 查看模型结果与方差分析
summary(negbin_glmm)
Anova(negbin_glmm)

如果你的数据是零膨胀计数数据，可添加零膨胀公式，例如ziformula = ~1表示零膨胀部分仅含截距

3. 运行emmeans事后检验 #

重新拟合模型后，执行你的原代码即可正常输出结果：

# 按周数分组，计算各处理的边际均值
emm <- emmeans(negbin_glmm, ~ tre | wa)
# 两两比较并做Tukey校正
pair_results <- pairs(emm, adjust = "tukey")

# 查看结果
pair_results

绘图用显著性标记提取 #

若要获取用于绘图的显著性信息：

• 生成字母标记（组间差异可视化常用）：

  cld_results <- cld(emm, adjust = "tukey", Letters = letters)
  

• 提取两两比较的p值：

  p_value_table <- as.data.frame(pair_results)[, c("contrast", "wa", "p.value")]
  

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