2

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?

goldisfine
  • 4,742
  • 11
  • 59
  • 83
  • 3
    have you looked at http://glmm.wikidot.com/faq ? – Ben Bolker Mar 26 '15 at 01:26
  • have you though about grabbing the data you need from m1 to create your own simulations of this model. I'd imagine you could create your own sim function from there. good luck. – miles2know Mar 26 '15 at 01:40

1 Answers1

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.

jknowles
  • 479
  • 3
  • 13