R 中嵌套随机效应模型的模型预测？

我正在尝试拟合具有随机效应的混合模型：lmer基于数据集 df 的模型这个例子在这里 https://stackoverflow.com/a/76510072/7258020。但是，我遇到了一个错误，指出“times”参数无效。有关如何修复它的任何建议。

set.seed(123)
df <- data.frame(x.axis = rep(c("A","B","C","D"), each = 50),
                 y.axis = rnorm(1000, 1.2, 2),
                 rf1 = rep(c("W","X","Y","Z"), times = 50),
                 rf2 = rep(c("1","2","1","2"), times = 50))

(preds_lm <- (lmer(y.axis ~ x.axis + (1|rf1/rf2), data = df) |>
            predict(se = TRUE, newdata = list(x.axis = c("A","B","C","D"), 
                                              rf1 = c("W","X","Y","Z"),
                                              rf2 = c("1","2"))) |>
            as.data.frame() |>
            cbind(x.axis = c("A","B","C","D")) |>
            within(y.min <- fit - qnorm(0.975) * se.fit) |>
            within(y.max <- fit + qnorm(0.975) * se.fit) |>
            within(y <- fit))[c(5, 8, 6, 7)])

rep(0, nobs) 中的错误：无效的“times”参数另外：警告 消息：在predict.merMod(lmer(y.axis ~ x.axis + (1 | rf1/rf2))中，数据 = df), : 忽略未使用的参数

调试时，“取消管道”代码很有用，这样您就可以看到哪里出了问题（我认为有一些帮助程序包用于此目的，但不记得它们是什么......）

p1 <- lmer(y.axis ~ x.axis + (1|rf1/rf2), data = df)
## fine
p2 <- predict(p1, se = TRUE, newdata = list(x.axis = c("A","B","C","D"), 
                                  rf1 = c("W","X","Y","Z"),
                                  rf2 = c("1","2")))

所以我们可以看到问题发生在predict。第一个明显的问题是predict方法用于merMod对象（即由lmer: see ?predict.merMod）没有se参数（如果你想要预测的标准错误，你可能需要看看像这样的包ggeffects）；这会导致“忽略未使用的参数”警告。

然而，采取se = TRUEout 仍然给我们同样的错误......但是newdata数据框而不是列表似乎可以工作（这也在?predict.merMod).

newdat <- data.frame(x.axis = c("A","B","C","D"), 
                                              rf1 = c("W","X","Y","Z"),
                                              rf2 = c("1","2"))
p2 <- predict(p1, newdata = newdat)

但是，您的其余代码将无法工作，因为（与predict线性模型的方法），predict()不返回包含预测和标准错误的列表。

仔细一看，我不认为ggpredict会做你想做的事。您可以使用GLMM 常见问题解答中的食谱 https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#lme用于计算预测的 SE（请注意，这仅包含固定效应中的不确定性——剩下的就是一堆蠕虫）。

mm <- model.matrix(delete.response(terms(p1)), newdat)
se.fit <- sqrt(diag(mm %*% tcrossprod(vcov(p1),mm)))

或者，您可以使用参数引导（slow但准确）：

bb <- bootMer(p1, function(m) predict(m, newdat), nsim = 1000, seed = 101)
confint(bb)
     2.5 %   97.5 %
1 0.7713446 1.451530
2 1.1581111 1.825785
3 0.9152968 1.570501
4 0.7163327 1.446359