r - 具有嵌套随机效应的 lmer 模型的预测

Question

我目前正在尝试帮助一位同事，但我根本找不到解决方案。所以我希望其他人可以帮助我们。

我有一个数据集，其中包含使用不同研究设计评估的权重数据，针对不同研究中的不同物种(一项研究包括多种设计和多种物种)。我想调查权重与研究设计之间的关系，使用研究和物种作为嵌套随机效应。

模型看起来像这样并且运行良好:

m <- lmer(weight ~ design +(1|study/species), data=dataset)

我试图对不同的物种做出预测，但进行了一般性研究:我创建了一个新的 data.table new.dt，其中包含原始数据集的独特设计-物种组合，并为报告添加了一列。

new.dt <- unique(dataset[,.(design, species))
new.dt$study <- "xyz"

然后我使用预测函数并允许新级别。

new.dt$p <- predict(m, newdata=new.dt, re.form= NULL, allow.new.levels=TRUE) 

我没有得到错误，但我对设计中的每个物种都得到了相同的预测。

有没有办法保持嵌套随机效应的一部分的原始水平，而另一部分成为新的水平？

提前致谢!

更新 - 工作示例:此问题不依赖于数据集。

library(data.table)
library(lme4)

dt <- data.table(expand.grid(design=c("a", "b"), species=c("x", "y", "z"), report=c("1", "2", "3"), count=seq(1, 10, 1)))
dt$weight <- 0

dt[species=="x"]$weight <- rnorm(60, 70, 10)
dt[species=="y"]$weight <- rnorm(60, 80, 15)
dt[species=="z"]$weight <- rnorm(60, 90, 20)

dt[design=="a"]$weight <- dt[design=="a"]$weight- 0.1*dt[design=="a"]$weight

dt[report=="1"]$weight <- dt[report=="1"]$weight+0.15*dt[report=="1"]$weight
dt[report=="2"]$weight <- dt[report=="2"]$weight-0.15*dt[report=="1"]$weight

m <-lmer(weight~design+(1|report/species), data=dt)

dt.pred <- unique(dt[,c(1:2)])
dt.pred$report<- "xyz"
dt.pred$pred<-predict(m, newdata=dt.pred, re.form= NULL, allow.new.levels=TRUE) 

Answer 1

“相同性”来自于您正在设置 re.form = NULL 或等效的 re.form = ~ 0。

线性混合效应模型对 Y|beta,b ~ intercept + X %*% beta + Z %*% b + e 进行建模，并通过设置 re.form = NULL 您正在预测期间设置 Z %*% b = 0 的定义。由于这是模型的随机部分(即 (1|report/species) )，您要删除 species 和 report 的随机效应>。

在混合模型中，您会称这种预测为“无条件预测”(或边际预测)[虽然它在实践中更像是伪无条件的]。它通常用于随机效应包含 individual 的模型中。在这种情况下，当您观察一个新个体时，您会产生未知的随机效应，但根据您的研究，您可能只对“系统”或“固定”效应感兴趣(即，该个体在被汽车？他骑自行车吗？)。在这里只看无条件/边际效应是有意义的。

换句话说，通过设置 re.form = NULL 就是说 Z %*% b = 0。由于物种是 Z 的一部分，权重为 b，因此您无法看到特定物种对预测的影响。

只有当您了解物种并且可以在预测中使用随机效应时，您才能在具有相同固定效应的物种之间得到不同的预测。

Ps.data.table 包具有与 expand.grid 等效的功能，称为 CJ，对于更大的集合，它会更快，内存更大高效的。