R语言 具有相关结构的混合模型从nlme到glmmTMB的转换

fnvucqvd  于 2023-04-27  发布在  其他
关注(0)|答案(1)|浏览(141)

我想使用predict在一些构建在nlme中的平均模型上绘制建模关系的置信区间。但是,我发现使用nlmeMuMIn::model.avg是不可能的。相反,我计划使用glmmTMB,就像建议的here一样。然而,我正在努力解决如何在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平均)中引导置信区间,请提出任何想法/建议。谢谢!

ewm0tg9j

ewm0tg9j1#

这里有两件重要的事情:

  • 你需要使用numFactor()而不是factor:对于一维结构(例如时间),这基本上只是使变量成为一个因子,其水平与唯一值相对应(与使用factor相反,factor会给你一个超过一百万个水平的变量......)
  • 正如Kasper Kristensen(TMB/glmmTMB的作者)所指出的那样,您应该使用ou()(Ornstein-Uhlenbeck)来表示相关性在 * 时间 * 中的指数衰减;exp()表示 * 空间 * 中相关性的指数衰减(并且慢得多...)

这是可行的:

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适合由随机效应定义的组内的相关结构,所以:

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),你会发现除了固定效应之外的大多数参数的置信区间覆盖了一个很大的范围(例如,对于id水平的随机效应SD,从2 e-17到8 e +15),对应于无法识别的参数。(希望这是因为你只给了我们一小部分数据,而不会发生在你的真实的问题中。)
(Hope更新下面的西姆斯人生使用ou()而不是exp()不久...)

update:看起来这个模型(使用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")

相关问题