我想从 lmer() 模型中获得预测的预测区间。我发现了一些关于此的讨论：
http://rstudio-pubs-static.s3.amazonaws.com/24365_2803ab8299934e888a60e7b16113f619.html
但他们似乎没有考虑随机效应的不确定性。
这是一个具体的例子。我在赛金鱼。我有过去 100 场比赛的数据。我想预测第 101 个，考虑到我的 RE 估计和 FE 估计的不确定性。我包括鱼的随机截距（有 10 种不同的鱼），以及重量的固定效果（较轻的鱼更快）。
library("lme4")
fish <- as.factor(rep(letters[1:10], each=100))
race <- as.factor(rep(900:999, 10))
oz <- round(1 + rnorm(1000)/10, 3)
sec <- 9 + rep(1:10, rep(100,10))/10 + oz + rnorm(1000)/10
fishDat <- data.frame(fishID = fish,
raceID = race, fishWt = oz, time = sec)
head(fishDat)
plot(fishDat$fishID, fishDat$time)
lme1 <- lmer(time ~ fishWt + (1 | fishID), data=fishDat)
summary(lme1)
现在，预测第 101 场比赛。鱼已经称重，准备出发：
newDat <- data.frame(fishID = letters[1:10],
raceID = rep(1000, 10),
fishWt = 1 + round(rnorm(10)/10, 3))
newDat$pred <- predict(lme1, newDat)
newDat
fishID raceID fishWt pred
1 a 1000 1.073 10.15348
2 b 1000 1.001 10.20107
3 c 1000 0.945 10.25978
4 d 1000 1.110 10.51753
5 e 1000 0.910 10.41511
6 f 1000 0.848 10.44547
7 g 1000 0.991 10.68678
8 h 1000 0.737 10.56929
9 i 1000 0.993 10.89564
10 j 1000 0.649 10.65480
鱼 D 真的放手了（1.11 盎司），实际上预计会输给鱼 E 和鱼 F，这两个鱼他都比过去好。但是，现在我想说，“鱼 E（重 0.91 盎司）将以概率 p 击败鱼 D（重 1.11 盎司）。” 有没有办法使用 lme4 做出这样的声明？我希望我的概率 p 考虑到我在固定效应和随机效应中的不确定性。
谢谢！
PS查看predict.merMod文档，它建议“没有计算预测标准误差的选项，因为很难定义一种在方差参数中包含不确定性的有效方法；我们建议bootMer执行此任务，”但是天哪，我看不到如何使用bootMer来做到这一点。它似乎bootMer可用于获取参数估计的自举置信区间，但我可能是错的。
更新问：
好吧，我想我问错了问题。我想能够说，“鱼 A，重量为 w oz，将有 90% 的时间是 (lcl, ucl) 的比赛时间。”
在我列出的示例中，重 1.0 盎司的鱼 A 的9 + 0.1 + 1 = 10.1 sec平均比赛时间为 0.1，标准偏差为 0.1。因此，他观察到的比赛时间将介于
x <- rnorm(mean = 10.1, sd = 0.1, n=10000)
quantile(x, c(0.05,0.50,0.95))
5% 50% 95%
9.938541 10.100032 10.261243
90% 的时间。我想要一个试图给我答案的预测函数。设置 all fishWt = 1.0in newDat，重新运行 sim 并使用（如下 Ben Bolker 所建议）
predFun <- function(fit) {
predict(fit,newDat)
}
bb <- bootMer(lme1,nsim=1000,FUN=predFun, use.u = FALSE)
predMat <- bb$t
给
> quantile(predMat[,1], c(0.05,0.50,0.95))
5% 50% 95%
10.01362 10.55646 11.05462
这似乎实际上以人口平均值为中心？好像没有考虑到 FishID 效应？我想这可能是样本量的问题，但是当我将观察到的比赛数量从 100 增加到 10000 时，我仍然得到类似的结果。
我会注意默认bootMer使用use.u=FALSE。另一方面，使用
bb <- bootMer(lme1,nsim=1000,FUN=predFun, use.u = TRUE)
给
> quantile(predMat[,1], c(0.05,0.50,0.95))
5% 50% 95%
10.09970 10.10128 10.10270
该区间太窄，似乎是 Fish A 平均时间的置信区间。我想要鱼 A 观察到的比赛时间的置信区间，而不是他的平均比赛时间。我怎么能得到那个？
更新 2，几乎：
我以为我在Gelman and Hill (2007)的第 273 页找到了我想要的东西。需要使用该arm软件包。
library("arm")
对于鱼 A：
x.tilde <- 1 #observed fishWt for new race
sigma.y.hat <- sigma.hat(lme1)$sigma$data #get uncertainty estimate of our model
coef.hat <- as.matrix(coef(lme1)$fishID)[1,] #get intercept (random) and fishWt (fixed) parameter estimates
y.tilde <- rnorm(1000, coef.hat %*% c(1, x.tilde), sigma.y.hat) #simulate
quantile (y.tilde, c(.05, .5, .95))
5% 50% 95%
9.930695 10.100209 10.263551
对于所有的鱼：
x.tilde <- rep(1,10) #assume all fish weight 1 oz
#x.tilde <- 1 + rnorm(10)/10 #alternatively, draw random weights as in original example
sigma.y.hat <- sigma.hat(lme1)$sigma$data
coef.hat <- as.matrix(coef(lme1)$fishID)
y.tilde <- matrix(rnorm(1000, coef.hat %*% matrix(c(rep(1,10), x.tilde), nrow = 2 , byrow = TRUE), sigma.y.hat), ncol = 10, byrow = TRUE)
quantile (y.tilde[,1], c(.05, .5, .95))
5% 50% 95%
9.937138 10.102627 10.234616
实际上，这可能不是我想要的。我只考虑了整体模型的不确定性。例如，在我观察到 Fish K 的 5 场比赛和 Fish L 的 1000 场比赛的情况下，我认为与我对 Fish K 的预测相关的不确定性应该比我对 Fish L 的预测相关的不确定性大得多。
将进一步研究 Gelman 和 Hill 2007。我觉得我可能最终不得不切换到 BUGS（或 Stan）。
更新第三个：
也许我对事物的概念化很糟糕。使用predictInterval()Jared Knowles 在下面的答案中给出的函数给出的间隔并不是我所期望的......
library("lattice")
library("lme4")
library("ggplot2")
fish <- c(rep(letters[1:10], each = 100), rep("k", 995), rep("l", 5))
oz <- round(1 + rnorm(2000)/10, 3)
sec <- 9 + c(rep(1:10, each = 100)/10,rep(1.1, 995), rep(1.2, 5)) + oz + rnorm(2000)
fishDat <- data.frame(fishID = fish, fishWt = oz, time = sec)
dim(fishDat)
head(fishDat)
plot(fishDat$fishID, fishDat$time)
lme1 <- lmer(time ~ fishWt + (1 | fishID), data=fishDat)
summary(lme1)
dotplot(ranef(lme1, condVar = TRUE))
我添加了两条新鱼。Fish K，我们观察了 995 场比赛，Fish L，我们观察了 5 场比赛。我们已经观察了 100 场 Fish AJ 的比赛。我和以前一样适合lmer()。从包装上看dotplot()：lattice
默认情况下，dotplot()按点估计对随机效应重新排序。Fish L 的估计值位于顶部，并且具有非常宽的置信区间。Fish K 在第三条线上，置信区间非常窄。这对我来说很有意义。我们有很多关于 Fish K 的数据，但没有很多关于 Fish L 的数据，所以我们对 Fish K 真实游泳速度的猜测更有信心。现在，我认为这会导致 Fish K 的预测区间变窄，而 Fish L 的预测区间变宽predictInterval()。豪瓦：
newDat <- data.frame(fishID = letters[1:12],
fishWt = 1)
preds <- predictInterval(lme1, newdata = newDat, n.sims = 999)
preds
ggplot(aes(x=letters[1:12], y=fit, ymin=lwr, ymax=upr), data=preds) +
geom_point() +
geom_linerange() +
labs(x="Index", y="Prediction w/ 95% PI") + theme_bw()
所有这些预测区间的宽度似乎都相同。为什么我们对 Fish K 的预测没有缩小其他预测？为什么我们对 Fish L 的预测不比其他人更广泛？