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

Question

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

3

我目前正试图帮助同事，但我无法找到解决方案。所以我希望有人能够帮助我们。

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

模型如下，并且运行良好：

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

我试图为不同的物种做出预测，但只进行了一项通用研究：我创建了一个名为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)

- I Ho

你可能需要查看 merTools 包，通过模拟生成混合效应模型的预测结果。文档链接在这里。 - ulfelder

“m”中“species”的方差项大小是多少？ - user20650

var(ranef(m)$species:study)的值为8.15——如果这是问题的正确答案。 - I Ho

不，我的意思是输出或VarCorr(m)返回的方差 -- 只是一个想法，如果方差非常小/接近零，那么预测结果可能没有什么差异。只是一个想法。 - user20650

这似乎与 re.form=NA 或 re.form=~0 相同，即没有随机效应。我想知道这是否相关：“对于以前未观察到的水平（或NAs）的数据，预测将使用无条件（人口水平）值。” - user20650

2个回答

0

您可以使用ggeffects包，它允许您获取固定效应的预测（包括置信区间）或在随机效应组水平上进行条件预测（这里不返回置信区间）。

以下是一个使用您的数据的示例，更多示例可以在此vignette中找到。

library(data.table)
library(lme4)
#> Loading required package: Matrix

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)

library(ggeffects)
ggpredict(m, "design")
#> 
#> # Predicted values of weight
#> # x = design
#> 
#> x | Predicted |   SE |         95% CI
#> -------------------------------------
#> a |     72.64 | 6.57 | [59.77, 85.52]
#> b |     82.66 | 6.57 | [69.78, 95.54]
#> 
#> Adjusted for:
#> * species = 0 (population-level)
#> *  report = 0 (population-level)

ggpredict(m, c("design", "report"), type = "re")
#> 
#> # Predicted values of weight
#> # x = design
#> 
#> # report = 1
#> 
#> x | Predicted
#> -------------
#> a |     80.78
#> b |     90.80
#> 
#> # report = 2
#> 
#> x | Predicted
#> -------------
#> a |     64.91
#> b |     74.92
#> 
#> # report = 3
#> 
#> x | Predicted
#> -------------
#> a |     72.24
#> b |     82.26
#> 
#> Adjusted for:
#> * species = 0 (population-level)

plot(ggpredict(m, c("design", "report"), type = "re"))
#> Loading required namespace: ggplot2

^{2020年02月07日使用reprex package (v0.3.0)创建}

- Daniel

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Oliver · Accepted Answer

“相同性”来自于你设置了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。由于species是带有权重b的Z的一部分，因此您无法看到对预测的物种特定效应。

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

附注： data.table包具有与expand.grid等价的函数，称为CJ，对于较大的数据集来说，这个函数更快且更节省内存。