R: Plotting lmer confidence intervals per faceted group 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.

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 <- ...

