将具有相关结构的混合模型从 nlme 转换为 glmmTMB

Question

我想在

predict

中构建的一些平均模型上使用

nlme

来绘制建模关系的置信区间。但是，我发现使用

nlme

和

MuMIn::model.avg

是不可能的。相反，我打算按照

here

的建议使用glmmTMB。但是，我正在努力研究如何在

glmmTMB

中设置相关结构。

以下是我的一小部分数据，以及

nlme

中的模型说明。数据是一个不完整的时间序列，随机结构是给定ID在序列中的测试位置，嵌套在ID中。

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)

我尝试按照此vignette中的说明进行操作，将时间转换为以单位间隔时间点为级别（在本例中为毫秒）的因子，并设置单个分组因子：

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

我的猜测是问题是时间序列不完整，但我不确定。

关于我是否/如何正确地将模型从

nlme

转换为

glmmTMB

的任何想法/建议，或者关于我如何从平均

nlme

模型（使用

MuMIn::model.avg

平均）引导置信区间失败的任何想法/建议都会非常有用赞赏。谢谢！

Answer 1

有两点很重要：

• 您需要使用

  numFactor()

  而不是

  factor

  ：对于一维结构（例如时间），这基本上只是使您的变量成为一个水平与唯一值相对应的因子（与您使用

   相比factor

  ，它给你一个超过百万级别的变量...）
• 正如 Kasper Kristensen（TMB/glmmTMB 的作者）所指出的，您应该使用

  ou()

  (Ornstein-Uhlenbeck) 来计算 time 中相关性的指数衰减；

  exp()

  用于 space 中相关性的指数衰减（并且要慢得多......）

所以这有效：

mydata$times <- numFactor(mydata$time)
mydata$group <- 1
model.glmmTMB <- glmmTMB(a ~ b + c + ou(times + 0 | group) + (1|id/testforid), 
                data = mydata)

但它与

lme

模型拟合度不太对应（甚至抛开使用

metric = "maximum"

的问题，我认为这在当前版本的

glmmTMB

中是不可能的）。

lme

符合相关结构 within 由随机效应定义的组，因此：

model.glmmTMB <- glmmTMB(a ~ b + c + ou(times + 0 | id/testforid),
    data = mydata)

更近了。 （您不需要

nugget = TRUE

，因为

glmmTMB

默认包含一个残差项，除非您使用

dispformula = ~0

将其关闭[对应于

nugget = FALSE

]。）

这会向您提供有关非正定 Hessian 矩阵的警告消息。然而，这实际上也与

lme

结果相匹配：如果您运行

intervals(models.lme)

，您会发现除了固定效应之外的大多数参数的置信区间覆盖了很大的范围（例如 2e-17 到 8e+15

id

水平的随机效应 SD），对应于无法识别的参数。 （希望这是因为您只给了我们一小部分数据，并且不会发生在您真正的问题上。）

---

（希望尽快更新下面的模拟人生以使用

ou()

而不是

exp()

...）

更新：看起来这个模型的计算成本（带有

ou()

）大约为（独特时间点的数量）^2.5。在我的机器上，如果不打开并行化（这可能有帮助也可能没有帮助——我怀疑代码的相关部分没有并行化），运行 1500 次观察（和 1500 次唯一时间）需要 45 秒。

您还可以尝试四舍五入您的时间值，以便唯一时间值的数量更少......

library(glmmTMB)
form <- a ~ b + c + ou(times + 0 | id)

## n should be a factor of 5
simfun <- function(n, round_times = FALSE, seed = 101) {
    if (!is.null(seed)) set.seed(seed)
    bigdata <- data.frame(b = runif(n, 0.001, 0.1),
                          c = sample(0:1, n, replace = TRUE),
                          time = c(10, 60, 150, 900, 1850)*runif(n, 0.9, 1.1),
                          id = factor(rep(seq(n/5), each = 5)))
    if (round_times) bigdata$time <- round(bigdata$time)
    bigdata$times <- numFactor(bigdata$time)
    bigdata$a <- simulate_new(RHSForm(form, as.form = TRUE),
                              ## show_pars = TRUE,
                              newdata = bigdata,
                              newparams = list(beta = c(35, 100, 1),
                                               betad = 1,
                                               theta = c(1,1)))[[1]]
    bigdata
}

nvec <- seq(50, 1500, by = 50)
pb <- txtProgressBar(max = length(nvec), style = 3)
elapsed <- rep(NA, length(nvec))
for (i in seq_along(nvec)) {
    setTxtProgressBar(pb, i)
    elapsed[i] <- system.time(simfun(nvec[i]))[["elapsed"]]
}
close(pb)

plot(nvec, elapsed, log = "xy")
lm(log(elapsed) ~ log(nvec))

elapsed_rnd <- n_unique <- rep(NA, length(nvec))

for (i in seq_along(nvec)) {
    setTxtProgressBar(pb, i)
    elapsed_rnd[i] <- system.time(res <- simfun(nvec[i], round_times = TRUE))[["elapsed"]]
    n_unique[i] <- length(unique(res$time))
}
close(pb)
lm(log(elapsed_rnd) ~ log(n_unique))
plot(n_unique, elapsed_rnd, log = "xy")