merTools predictInterval() 用于具有嵌套随机效应的模型

merTools predictInterval() for model with nested random effect

merTools 包中的 predictInterval() 是否不像嵌套随机效应?例如,使用 ggplot2 包中的 msleep 数据集:

library("lme4")
library("merTools")
library("ggplot2")
mod <- lmer(sleep_total ~ bodywt + (1|vore/order), data=msleep)
predInt <- predictInterval(merMod=mod, newdata=msleep) 

Returns 一个错误:

Error in '[.data.frame'(newdata, , j) : undefined columns selected

这运行良好没问题:

mod <- lmer(sleep_total ~ bodywt + (1|vore) + (1|order), data=msleep)
predInt <- predictInterval(merMod=mod, newdata=msleep)

(实际上它给出了关于随机效应变量中 NA 水平的警告,但我并不担心)

更新

正如下面 Ben Bolker 的回答的评论中所讨论的,merTools 的新版本说明了嵌套随机效应。但是,当我尝试预测包含新级别嵌套随机效应的数据时,出现错误。

这个有效:

mod <- lmer(sleep_total ~ bodywt + (1|vore/order), data=msleep)
predInt <- predictInterval(merMod=mod, newdata=msleep) 

虽然有几个警告,但它有效(请参阅下面有关警告的其他问题*):

mod <- lmer(sleep_total ~ bodywt + (1|vore) + (1|order), data=msleep)
msleep2 <- msleep %>% mutate(vore = "omni")
predInt <- predictInterval(merMod=mod, newdata=msleep2) 

但这不起作用:

mod <- lmer(sleep_total ~ bodywt + (1|vore/order), data=msleep)
msleep2 <- msleep %>% mutate(vore = "omni")
predInt <- predictInterval(merMod=mod, newdata=msleep2) 

出现以下错误:

Error in `[.data.frame`(tmp, alllvl) : undefined columns selected
In addition: Warning message:
In predictInterval(merMod = mod, newdata = msleep3) :
  newdata is tbl_df or tbl object from dplyr package and has been
              coerced to a data.frame

在这里,"omni" 实际上并不是 vore 的新级别,但当与 order 结合时,它会创建新的变量嵌套组合。

如果我使用 "new" 或任何其他未观察到的 vore 水平的东西,我会得到类似的结果:它适用于模型的非嵌套版本,但不适用于嵌套版本。


*此外,我是否应该关注上面第二个模型块给出的警告:

> mod <- lmer(sleep_total ~ bodywt + (1|vore) + (1|order), data=msleep)
> msleep2 <- msleep %>% mutate(vore = "omni")
> predInt <- predictInterval(merMod=mod, newdata=msleep2)
Warning messages:
  1: In predictInterval(merMod = mod, newdata = msleep2) :
     newdata is tbl_df or tbl object from dplyr package and has been
       coerced to a data.frame
  2: In chol.default(sigma, pivot = TRUE) :
     the matrix is either rank-deficient or indefinite

我猜第二个是 vore 每次观察都采用相同值的结果,但这不应该成为预测的问题,不是吗?如果变量在我拟合模型时取相同的值,我可以认为这是一个问题,但我认为在预测新观察时这不应该是一个问题吗?

可以(显然)通过显式写出交互项来解决这个问题。 警告:我还没有实际检查以确保结果预测是正确的,只是看到没有产生错误并且结果对象大致合理...

msleep <- transform(msleep,voreOrder=interaction(vore,order,drop=TRUE))
mod2 <- lmer(sleep_total ~ bodywt + (1|vore)+(1|voreOrder), data=msleep)
predInt <- predictInterval(merMod=mod2, newdata=msleep) 

这确实会生成警告消息,但显然它们是由于 vore 变量中的 <NA> 值(我不知道这个数据集...)