【发布时间】:2020-02-03 11:56:20
【问题描述】:
我目前正在尝试帮助一位同事,但我根本找不到解决方案。所以我希望其他人可以帮助我们。
我有一个数据集,其中包含针对不同研究中不同物种的不同研究设计评估的体重数据(一项研究包括多个设计和多个物种)。我想研究权重和研究设计之间的关系,使用研究和物种作为嵌套随机效应。
模型看起来像这样并且运行良好:
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)
【问题讨论】:
-
您可能需要查看
merTools包,以通过模拟从混合效应模型生成预测。文档here. -
m中species的方差项的大小是多少? -
var(ranef(m)$
species:study) 是 8.15 - 如果这是这个问题的正确答案。 -
不,我的意思是输出中显示的方差或
VarCorr(m)返回的方差——这只是一个想法,以防方差很小/几乎为零,那么预测中就没有什么不同了。只是一个想法。 -
这似乎给出了与
re.form=NA或re.form=~0相同的结果,即没有随机效应。我想知道这是否相关:“预测将使用以前未观察到的水平(或 NA)的数据的无条件(人口水平)值。”
标签: r predict lme4 mixed-models random-effects