如何用R中brms建模时间自相关？解决组内时间点不唯一错误

问题背景 #

现有两个位点7个物种的eDNA qPCR相对密度数据集，样本在两天内采集，怀疑采样设备残留DNA导致样本间存在自相关。使用brms包的brm()函数构建gamma hurdle模型（处理零膨胀），以物种ID为随机效应，目前遇到两个核心问题：

1. 如何通过残差图判断自相关存在的证据？
2. 如何在模型中纳入自相关结构，同时解决Error: Time points within groups must be unique的报错？

可复现数据集：

data <- data.frame(site = c(rep("siteA", times=35), rep("siteB", times=35)),
                  seq = c(rep(1:5, each=7), rep(6:10, each=7)),
                  species = rep(c("a", "b", "c", "d","e", "f", "g")),
                  dnaquant = rgamma(70, 5))

其中seq为采样顺序，原始无自相关的模型可正常运行：

m4 <- brm(bf(dnaquant ~ site + (1 + site|species), hu ~ site), 
          data = data, family = hurdle_gamma(), chain = 2, cores = 2)

尝试添加自相关组件时触发报错：

m4 <- brm(bf(dnaquant ~ site + (1 + site|species), hu ~ site) + acformula(~arma(time = seq, gr = species, cov=TRUE)), data = data, family = hurdle_gamma(), chain = 2, cores = 2)

错误信息：Error: Time points within groups must be unique

疑问1：解读残差图寻找自相关证据 #

通过以下步骤分析残差判断自相关：

• 提取条件残差：使用residuals(m4, type = "conditional")获取扣除随机效应后的残差，更能反映模型未解释的变异。
• 残差时序图：绘制残差随seq（采样顺序）的折线图。若存在自相关，残差会呈现连续的正/负集群（如连续多个正残差后接连续负残差），而非完全随机的散点分布。
• 自相关函数(ACF)图：对残差运行acf(residuals(m4, type = "conditional"))，若滞后1阶或多阶的自相关系数超出置信区间（显著不为0），则说明存在自相关。

疑问2：纳入自相关结构并解决报错 #

报错原因 #

acformula中指定gr = species时，要求每个物种分组内的time（即seq）必须唯一，但你的数据中同一个seq对应所有7个物种，导致每个物种的seq值重复，触发报错。

解决方案 #

根据设备残留的自相关逻辑（样本间交叉污染），推荐以下两种合理建模方式：

方案1：全局采样顺序自相关（推荐）

自相关源于采样设备残留，是样本间的关联（相邻采样的样本无论物种都存在污染），因此直接基于全局seq建模自相关，无需按物种分组：

m4_ac <- brm(
  bf(dnaquant ~ site + (1 + site|species), hu ~ site) + 
    acformula(~arma(time = seq, cov=TRUE)),
  data = data, 
  family = hurdle_gamma(), 
  chain = 2, 
  cores = 2
)

方案2：按样本分组建模自相关

每个seq对应一个独立样本，所有物种共享该样本的自相关结构，可通过创建样本ID分组实现：

# 创建唯一样本ID
data$sample_id <- as.factor(data$seq)

# 基于样本ID分组建模自相关
m4_ac_sample <- brm(
  bf(dnaquant ~ site + (1 + site|species), hu ~ site) + 
    acformula(~arma(time = seq, gr = sample_id, cov=TRUE)),
  data = data, 
  family = hurdle_gamma(), 
  chain = 2, 
  cores = 2
)

备选：物种内时间自相关（仅适用于物种自身时序关联）

若需建模同一物种在不同采样时间的自相关（非设备残留场景），可创建物种-采样顺序的唯一标识作为时间变量：

# 创建物种-采样顺序唯一ID并转为数值型
data$sp_seq <- as.numeric(factor(paste(data$species, data$seq, sep = "_")))

# 按物种分组，用唯一ID作为时间变量
m4_ac_sp <- brm(
  bf(dnaquant ~ site + (1 + site|species), hu ~ site) + 
    acformula(~arma(time = sp_seq, gr = species, cov=TRUE)),
  data = data, 
  family = hurdle_gamma(), 
  chain = 2, 
  cores = 2
)

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