如何将带相关结构的nlme混合模型转换为glmmTMB模型？

背景 #

原本使用nlme构建含时间自相关的混合效应模型，希望通过MuMIn::model.avg得到模型平均结果后调用predict绘制置信区间，但nlme的模型平均结果无法直接实现该功能，计划转用glmmTMB，但不清楚如何对应设置模型结构；若无法转换，希望了解如何通过bootstrap获取模型平均结果的置信区间。

原始nlme模型及数据 #

数据为不完整时间序列，随机结构为嵌套在id下的testforid，同时设置指数自相关结构（含块金效应）：

library(nlme)
library(glmmTMB)

mydata <- structure(list(id = c("F530", "F530", "F530", "F530", "F530", "M391", "M391", "M391", "M391", "M391", "M391", "M391"),testforid = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), levels = c("1", "2"), class = c("ordered", "factor")), time = c(12.043, 60.308, 156.439, 900.427, 1844.542, 42.095, 61.028, 130.627, 194.893, 238.893, 905.282, 1859.534), a = c(35.5786398928957, 35.4973671656257, 36.7414694383557, 37.4316029157078, 36.0805603474457, 38.892219234833, 37.081136308003, 37.339272893363, 36.744902161663, 36.741897283613, 38.158072893363, 38.946697283613), b = c(0.0079975108148372, 0.0151689857479705, 0.0275942757878888, 0.0125676102827941, 0.0352227834243443, 0.0195902976534779, 0.0118588484445401, 0.0069799148425349, 0.00723445099500534, 0.00787758751826021, 0.0162518412492866, 0.0127526068249484), c = c(1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0)), row.names = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L), class = "data.frame")

model.lme <- lme(a ~ b + c,
                 random = list(id = ~1, testforid = ~1),
                 correlation = corExp(metric = "maximum", nugget = TRUE),
                 method = "ML",
                 data = mydata)

错误的glmmTMB尝试及问题 #

尝试将time转为毫秒间隔因子并设置分组变量，构建模型时出现错误：

# 错误的转换代码
mydata$times <- factor(mydata$time,
                       levels = seq(from = min(mydata$time),
                                    to = max(mydata$time),
                                    by = 0.001))
mydata$group <- 1

# 错误的模型代码
model.glmmTMB <- glmmTMB(a ~ b + c + exp(times + 0 | group) + (1|id/testforid), data = mydata)

报错信息：

Error in parseNumLevels(reTrms$cnms[[i]]) : 
  Failed to parse numeric levels: times12.043times42.095times60.308times61.028times130.627times156.439times194.893times238.893times900.427times905.282times1844.542times1859.534
In addition: There were 12 warnings (use warnings() to see them)
> warnings()
Warning messages:
1: In lapply(strsplit(tmp, ","), as.numeric) : NAs introduced by coercion

解决方案 #

一、正确转换为glmmTMB模型 #

glmmTMB中处理时间自相关（指数结构+块金效应），不需要将time转为因子，而是通过corr参数直接指定自相关结构，同时正确设置嵌套随机效应：

1. 模型结构说明：

   • 固定效应：a ~ b + c
   • 嵌套随机效应：(1|id/testforid)（等价于(1|id) + (1|id:testforid)）
   • 指数自相关结构（含块金效应）：通过corExp(form = ~time, nugget = TRUE)基于原始时间变量计算自相关

2. 正确代码：

# 加载必要包
library(glmmTMB)
library(nlme) # 用于corExp结构

# 构建glmmTMB模型
model.glmmTMB <- glmmTMB(
  a ~ b + c,
  random = ~ (1|id/testforid),
  corr = corExp(form = ~time, nugget = TRUE), # 指定时间自相关+块金效应
  method = "ML",
  data = mydata
)

# 查看模型结果
summary(model.glmmTMB)

二、模型平均后获取置信区间（若保留nlme模型） #

如果无法转用glmmTMB，可通过bootstrap对MuMIn::model.avg的结果进行抽样，获取预测值的置信区间：

1. 步骤说明：

   • 构建候选模型集并完成模型平均
   • 自定义bootstrap函数，每次抽样后重新拟合模型、计算模型平均预测值
   • 基于bootstrap结果计算95%置信区间

2. 示例代码：

library(MuMIn)
library(nlme)
library(boot)

# 1. 构建候选模型集（可根据需求扩展）
model1 <- lme(a ~ b + c, random = list(id = ~1, testforid = ~1), correlation = corExp(metric = "maximum", nugget = TRUE), method = "ML", data = mydata)
model2 <- lme(a ~ b, random = list(id = ~1, testforid = ~1), correlation = corExp(metric = "maximum", nugget = TRUE), method = "ML", data = mydata)
model3 <- lme(a ~ c, random = list(id = ~1, testforid = ~1), correlation = corExp(metric = "maximum", nugget = TRUE), method = "ML", data = mydata)

# 2. 模型平均（使用AIC权重）
model_avg <- model.avg(model1, model2, model3)

# 3. 自定义bootstrap函数
boot_pred <- function(data, indices) {
  # 抽样数据
  d <- data[indices, ]
  # 重新拟合模型（跳过拟合失败的情况）
  m1 <- try(lme(a ~ b + c, random = list(id = ~1, testforid = ~1), correlation = corExp(metric = "maximum", nugget = TRUE), method = "ML", data = d), silent = TRUE)
  m2 <- try(lme(a ~ b, random = list(id = ~1, testforid = ~1), correlation = corExp(metric = "maximum", nugget = TRUE), method = "ML", data = d), silent = TRUE)
  m3 <- try(lme(a ~ c, random = list(id = ~1, testforid = ~1), correlation = corExp(metric = "maximum", nugget = TRUE), method = "ML", data = d), silent = TRUE)
  
  if(inherits(m1, "try-error") || inherits(m2, "try-error") || inherits(m3, "try-error")) return(NA)
  
  # 重新计算模型平均并生成预测值
  avg <- model.avg(m1, m2, m3)
  predict(avg, newdata = data)
}

# 4. 运行bootstrap（设置抽样次数，示例为1000次）
set.seed(123)
boot_result <- boot(data = mydata, statistic = boot_pred, R = 1000)

# 5. 计算95%置信区间
ci <- apply(boot_result$t, 2, function(x) quantile(x, c(0.025, 0.975), na.rm = TRUE))

# 查看置信区间
print(ci)

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