R:使用 lmer 和预测(或其他)模拟数据

R: simulating data with lmer and predict (or else)

我正在拟合以下模型

fit<- lmer(y ~ a + b + (1|c) + (1|a:d) , data=inputdata)

到 "inputdata" 中收集的真实观察结果。 现在我想根据模型参数和确定的误差为模拟生成各种(1000)建模数据集。我可以使用

pred <- predict(fit, newdata=list(a=val_a1, b=val_b1, c=val_c1, d = val_d1),
                allow.new.levels = TRUE)

但这总是提供相同的(最有可能的平均值)。有没有办法获得值的分布,即从预测分布中提取?

正如@Adam Quek 所问,一个可重现的例子:

    #creating dataset
    a <- as.factor(sort(rep(1:4,5 )))
    b <- rep(1:2,10)+0.5
    c <- as.factor(c( sort(rep(1:2,5)),sort(rep(1:2,5)) ))
    d <- as.factor(rep(1:5,4 ))
    a <- c(a,a,a)
    b <- c(b,b,b)
    c <- c(c,c,c)
    d <- c(d,d,d)
    y <- rnorm(60)
    inputdata = data.frame(y,a,b,c,d)

    # fitting the model   
    fit<- lmer(y ~ a + b + (1|c) + (1|a:d) , data=inputdata)

    # making specific predictions for a parameter set
    val_a1 = 1
    val_b1 = 2
    val_c1 = 1
    val_d1 = 4
    pred <- predict(fit, newdata=list(a=val_a1, b=val_b1, c=val_c1, d = val_d1),
                    allow.new.levels = TRUE)
    pred

我得到的是: 0.2394255

如果我再来一次

pred <- predict(fit, newdata=list(a=val_a1, b=val_b1, c=val_c1, d = val_d1),
                allow.new.levels = TRUE)
pred

我当然知道: 0.2394255

但我正在寻找的是一个 R 函数或例程,它可以轻松提供一组遵循我的输入值分布的预测。像

for (i in 1:1000){
    pred[i] <- predict(fit, newdata=list(a=val_a1, b=val_b1, c=val_c1, d = 
               val_d1),allow.new.levels = TRUE)
}

mean(pred) = 0.2394255但是sd(pred) != 0

感谢@Alex W! bootMer 完成这项工作。下面为那些对示例的解决方案感兴趣的人:

m1 <- function(.) {
   predict(., newdata=inputdata, re.form=NULL)
}
boot1 <- lme4::bootMer(fit, m1, nsim=1000, use.u=FALSE, type="parametric")
boot1$t[,1]

其中 boot1$t[,1] 现在包含使用 inputdata[1,] 中定义的参数值时的 1000 个预测。 https://cran.r-project.org/web/packages/merTools/vignettes/Using_predictInterval.html 很有帮助 link。