基于GLMM的鸟类样点计数数据中物种-年度随机效应建模咨询及代码指导

Great question—you’re absolutely on the right track with using random effects for shrinkage here! Let’s break this down clearly, since you’ve already done a ton of good groundwork.

First off: 用随机效应替代固定的species×year交互项完全合理且必要。你的核心问题是固定交互的参数数量爆炸（40个物种×33个年度=1320个独立参数），远远超出了数据能支撑的范围。随机效应的"信息借取"（shrinkage）会把每个物种-年度的估计向总体均值收缩，既保留物种特异性趋势，又避免模型过拟合、崩溃或得到极端不合理的结果，这正是这类多物种、稀疏结构化数据的标准解法。

一、合理的随机效应结构选择 #

你需要交叉随机效应来捕获物种和年度的独立组合效应，这里有两种适配你需求的常用结构：

1. 纯交叉截距效应（最简洁） #

如果你只需要每个物种-年度的独立丰度截距（不考虑年度对物种的斜率变化），用: 表示交叉交互的随机效应：

(1|species:year)

这个结构会为每个物种-年度组合生成一个随机截距，同时自动施加收缩，把极端值拉回合理区间。

2. 随机斜率-截距交叉效应（更灵活） #

如果想让每个物种的年度变化趋势（斜率）也随机变化（比如不同物种的种群增长/下降速率不同），用(species|year)——这是glmmTMB中"在year的每个水平下，species的截距和斜率都随机变化"的语法，等价于允许物种和年度的交互随机变化：

(species|year)

这个结构比纯交叉截距更灵活，你的26万条观测数据完全能支撑它的参数规模。

二、完整模型的代码实现 #

结合你的需求（加入site、week效应，零膨胀负二项分布），这里提供几个递进的可运行方案：

基础模型：聚焦物种-年度+site随机效应 #

# 基础零膨胀负二项模型，用纯交叉随机效应捕获物种-年度差异
m_base <- glmmTMB(
  count ~ (1|species:year) + (1|site),
  ziformula = ~ (1|species),  # 零膨胀部分考虑物种特异性的零概率
  family = nbinom2,
  data = df,
  control = glmmTMBControl(optCtrl = list(maxit = 1000))  # 增加迭代次数避免收敛警告
)

加入week效应的优化模型 #

考虑到week的非线性和物种迁徙特性，你可以选择两种week的参数化方式：

# 方案1：week作为固定效应（适合所有物种的周变化趋势类似的情况）
m_week_fixed <- glmmTMB(
  count ~ week + (1|species:year) + (1|site),
  ziformula = ~ (1|species),
  family = nbinom2,
  data = df,
  control = glmmTMBControl(optCtrl = list(maxit = 1000))
)

# 方案2：week效应随物种随机变化（适合不同物种迁徙时间差异大的情况）
m_week_random <- glmmTMB(
  count ~ (1|species:year) + (1|site) + (0 + week|species),
  ziformula = ~ (1|species),
  family = nbinom2,
  data = df,
  control = glmmTMBControl(optCtrl = list(maxit = 1000))
)

你测试过的(species|year)结构优化版 #

如果你想保留物种的年度随机斜率，这个结构完全合理，优化后代码如下：

m_slope <- glmmTMB(
  count ~ species + (species|year) + (1|site) + week,
  ziformula = ~ species,  # 零膨胀先试用固定物种效应（如果随机效应收敛困难）
  family = nbinom2,
  data = df,
  control = glmmTMBControl(optCtrl = list(maxit = 2000), nlminbwrap = list(eval.max = 2000))
)

这里把species作为固定效应，加上(species|year)随机效应，相当于每个物种有一个固定截距，同时年度对每个物种的效应是随机变化的，既保留物种特异性，又有收缩效应。

三、收敛问题的排查与优化 #

你遇到的迭代限制、Hessian非正定等警告，可以通过以下方法解决：

• 增加迭代次数：用control=glmmTMBControl(optCtrl=list(maxit=1000))甚至更高数值，默认迭代次数往往不够支撑复杂模型。
• 简化零膨胀部分：如果零膨胀的随机效应收敛困难，先换成固定效应（比如ziformula=~species），确认主模型收敛后再逐步调整。
• 检查极端值：极少数异常大的count值可能影响收敛，可以考虑截断，或者尝试nbinom1分布（nbinom1方差与均值线性相关，nbinom2是二次相关，可根据数据离散程度选择）。
• 从小模型测试：先拿10个物种跑通模型，确认结构合理后再扩展到全部40种，避免直接跑大模型浪费时间。
• 检查收敛细节：用checkConvergence(m_slope)函数（glmmTMB包内置）查看具体收敛问题，针对非正定的参数调整模型结构。

四、提取物种年度丰度估计 #

模型拟合完成后，用emmeans包可以轻松提取每个物种-年度的预测丰度：

library(emmeans)

# 提取物种-年度的原始尺度预测值（针对m_base模型）
yearly_abundance <- emmeans(m_base, ~ species:year, type = "response")
# 转换成数据框查看结果
as.data.frame(yearly_abundance)

注意指定type="response"，这样得到的是原始丰度尺度的估计，而非对数尺度。

五、经验总结 #

1. 随机效应收缩是多物种建模的标配：当因子水平组合过多时，固定交互几乎必然导致收敛问题，随机效应的收缩是唯一可行的方案，也是生态学家处理多物种动态数据的常用方法。
2. 模型复杂度循序渐进：先从最简单的结构（比如只加species:year和site）开始，确认收敛后再逐步加入week、零膨胀的复杂结构。
3. 选择合适的分布：nbinom2适合过度离散严重的数据，nbinom1适合方差与均值成比例的情况，如果nbinom2收敛困难，可以试试nbinom1或者泊松+零膨胀。
4. 零膨胀部分灵活调整：如果物种间的零概率差异很大，用随机效应；如果差异稳定，用固定效应，优先保证主模型的收敛。

内容的提问来源于stack exchange，提问作者G. Channing