Confidence Interval in mixed effect models

Question

library(lme4)

fm1 <- lmer(Reaction ~ Days + (Days|Subject), data = sleepstudy)

To generate a 95% CI, I can use the predictInterval() function from the package merTools.

library(merTools)
head(predictInterval(fm1, level = 0.95, seed = 123, n.sims = 100))
# fit      upr      lwr
# 1 255.4179 313.8781 184.1400
# 2 273.2944 333.2005 231.3584
# 3 291.8451 342.8701 240.8226
# 4 311.3562 359.2908 250.4980
# 5 330.3671 384.2520 270.7094
# 6 353.4378 409.9307 289.4760

In the documentation, it says about the predictInterval() function

This function provides a way to capture model uncertainty in predictions from multi-level models fit with lme4. By drawing a sampling distribution for the random and the fixed effects and then estimating the fitted value across that distribution, it is possible to generate a prediction interval for fitted values that includes all variation in the model except for variation in the covariance parameters, theta. This is a much faster alternative than bootstrapping for models fit to medium to large datasets.

My goal is to get all the fitted values instead of the the upper and lower CI i.e. for each row, I need the original n simulations from which these 95% CI are calculated. I checked the argument in the documentation and followed this:

head(predictInterval(fm1, n.sims = 100, returnSims = TRUE, seed = 123, level = 0.95))
# fit      upr      lwr
# 1 255.4179 313.8781 184.1400
# 2 273.2944 333.2005 231.3584
# 3 291.8451 342.8701 240.8226
# 4 311.3562 359.2908 250.4980
# 5 330.3671 384.2520 270.7094
# 6 353.4378 409.9307 289.4760

Instead of getting the 100 simulations, it still gives me the same output. What is it I am doing wrong here?

A second question though I believe this is more of a StatsExchange one.

"By drawing a sampling distribution for the random and the fixed effects and then."`

How does it draws the sampling distribution if some could explain me?

Lstat Lstat · Accepted Answer · 2018-07-10T13:05:32

You can get simulated values if you specify newdata in the predictInterval() function.

predInt <- predictInterval(fm1, newdata = sleepstudy, n.sims = 100,
           returnSims = TRUE, seed = 123, level = 0.95)

simValues <- attr(predInt, "sim.results")

Details on how to create sampling distributions of parameters are given in the Detail section of the help page.You can get the estimates of fit, lower and upper boundaries as:

fit <- apply(simValues, 1, function(x){quantile(x, probs=0.500) } )
lwr <- apply(simValues, 1, function(x){quantile(x, probs=0.025) } )
upr <- apply(simValues, 1, function(x){quantile(x, probs=0.975) } )

Confidence Interval in mixed effect models

1 Answers