How to get confidence intervals for lmer object?

Question

I am trying to get confidence intervals for predictions on the mixed model. The predict function does not output any confidence intervals. Few StackOverflow answers suggested using predictInterval function from the merTools package to obtain the intervals but there is a discrepancy between the prediction estimates from these two function which I am trying to compare in the plot below. Can someone kindly let me know what I am doing wrong here? Also, the actual model I am trying to build is similar to the one showed in the snippet below where I only have no fixed effect components except for the intercept.

library(merTools)
library(lme4)

dat <- iris
mod <- lmer(Sepal.Length ~ 1 + (1 + Sepal.Width + Petal.Length + 
                                  Petal.Width|Species), data=dat)

c1 <- predict(mod, dat)
c2 <- predictInterval(mod, dat)

plot_data <- cbind(c1, c2)
plot_data$order <- c(1:nrow(plot_data))

library(ggplot2)
ggplot(plot_data) + geom_line(aes(x=order, y=c1), color='red') + 
  geom_ribbon(aes(x=order, ymin=lwr, ymax=upr), color='blue', alpha=0.2) +  
  geom_line(aes(x=order, y=fit), color='blue')

The red line indicates the predictions 'c1' and the blue line indicates the predictions 'c2'

Weihuang Wong Weihuang Wong · Accepted Answer · 2018-02-08T14:11:04

I haven't quite been able to isolate the part of predictInterval that's causing the problem, but a solution to your specific problem is to note that if what you want is group-varying intercepts and slopes, then you can fit the following equivalent model

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

Now if we apply predictInterval to this fitted model

c2 <- predictInterval(mod2, dat)

keeping the rest of your example, we get the following plot:

which is what we want. (To emphasize, the red line represents predictions from your original model specification i.e. with only the intercept in the "fixed" component.)

How to get confidence intervals for lmer object?

1 Answers