R – How to Plot LMER Confidence Intervals per Faceted Group?

Question

I am using lme4 package to run a Mixed-Effects Model followed by the predict function ot obtain fitting lines per invidual level and group level. Yet, I am struggling to get the confidence interval of the fitting line per group level to represent them in a ggplot. My code is as follows:

The data: myData is a data frame containing:
totalVol: numeric vector // Years: numeric vector // TIV: numeric vector // Group: ordered factor with 5 levels:0<1<2<3<4 // subject: factor with 46 levels

#The model
ModelLME <- lmer(totalVol ~ Years*Group + TIV + (1|subject),data=myData)

#New data frame to keep one of the variables fixed for the group prediction lines
newdat = data.frame(totalvol=myData$totalVol, Years=myData$Years, Group=myData$Group, TIV=median(myData$TIV))

#Prediction lines
pred1a <-predict(ModelLME, newdata=newdat, re.form=NA) #predict at group level
pred1b <-predict(ModelLME) #predict at individual level

#Plot
plot <- ggplot(myData, aes(x=Years, y=totalVol)) + facet_grid(~Group) +
    geom_line(aes(y=pred1b, group=subject)) + 
    geom_smooth(aes(y=pred1a), method = "glm", size=1.5, color="orange") + theme_light() + geom_point(size=1, shape=22) +
    xlab("- Years") + ylab("Total Volume (mm3)") + theme(strip.text.x = element_text(size = 12, color ="black", face = "bold"), strip.background = element_rect(
    color="black", fill="light grey",linetype="solid")) 

This gives the following:

To add the confidence interval of the group prediction line, I have tried:

myData <- cbind(myData,predictInterval(ModelLME,which="fixed")

#add the interval to the plot
plotCI <- plot + geom_ribbon(aes(Years,ymin=lwr, ymax=upr,group=Group),alpha =.2)

But in this way, I get the confidence interval adjusted to each subject, not to each group:
.
I expected that restricting predictInterval to fixed effects only would give the interval per group, as occurs in predict. But my understanding of how these functions work is quite limited, so I am probably missing something. I have also tested the confintfunction but this gives a confidence interval per estimate of the model, so not sure if that is useful for the plot.

I will greatly appreciate any advice on how to represent the CI properly or other suggestions on better ways to approach this prediction! Thank you so much in advance!
Diana.

Answer 1

You may want to double-check this is correct... You could extract the means and CI of each group at a reference grid and plot the results. The package emmeans should do most of the job.

I use here the Orthodont dataset with a model similar to yours:

library(lme4)
library(ggplot2)
library(emmeans)

data(Orthodont,package="nlme")

fit <- lmer(distance ~ age * Sex + (1|Subject), data= Orthodont)

gr <- ref_grid(fit, cov.keep= c('age', 'Sex'))
emm <- emmeans(gr, spec= c('age', 'Sex'), level= 0.95)
emm
 age Sex    emmean    SE   df lower.CL upper.CL
   8 Male     22.6 0.539 37.1     21.5     23.7
  10 Male     24.2 0.492 26.3     23.2     25.2
  12 Male     25.8 0.492 26.3     24.7     26.8
  14 Male     27.3 0.539 37.1     26.2     28.4
   8 Female   21.2 0.650 37.1     19.9     22.5
  10 Female   22.2 0.594 26.3     20.9     23.4
  12 Female   23.1 0.594 26.3     21.9     24.3
  14 Female   24.1 0.650 37.1     22.8     25.4

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

gg <- ggplot(data= Orthodont, aes(x= age, y= distance)) +
    geom_ribbon(data= data.frame(emm), aes(ymin= lower.CL, ymax= upper.CL, y= NULL), fill= 'grey80') +
    geom_line(data= data.frame(emm), aes(y= emmean)) +
    geom_point() +
    facet_wrap(~Sex)

---

I haven't done this before but it seems to me that to get the means at points other than the original ones, you can use the argument at in ref_grid and proceed in the same way as above. For example, we want age from 0 to 24:

gr <- ref_grid(fit, at= list(age= 0:24), cov.keep= c('age', 'Sex'))
emm <- ...