How to get a prediction interval for a lmer object?

Question

I've surfed the web and haven't found a satisfactory answer for this.

How would I generate a prediction interval for from a lmer object for each observation in the test dataset?

train_ind <- sample(seq(1:nrow(iris)), size = nrow(iris)/2, replace = F)
TRAIN <- iris[train_ind,]
TEST <- iris[-train_ind,]

m1 <- lmer(Sepal.Length ~ Sepal.Width + Petal.Length + (1|Species), data = TRAIN)

interval has no argument with predict. sim apparently doesn't work so I need to use an mcmc function that comes with LMER and draw from quantiles?

Answer 1

We just released a package called merTools which eases this process. It does not provide a complete prediction interval because it skips the step of simulating the theta terms in the merMod, but it does produce an interval that accounts for the variation in the fixed effect coefficients, random effect coefficients, and the residual error in the model. It's also pretty fast.

library(merTools)
preds <- predictInterval(m1, n.sims = 500, level = 0.9, stat = 'median')
head(preds)

It's that easy, you can change the values you want returned, return all simulated values of yhat, and you can also adjust the width of the prediction interval you are interested in obtaining.