How to get 95% CI for each level of a random effect in a lmer mixed model

I have run several mixed models in lme4 with some main effects, a few interaction terms, and two random effects (animalID and year). My model looks like this:

model <- lmer(response ~ ns(area_scale, df = 3) + 
   area_scale:spat_scale + area_scale:temp_scale + 
   NDVI_scale:spat_scale + NDVI_scale:temp_scale + 
   spat_scale*temp_scale + NDVI_scale + (1|animal) + (1|year), 
   data = data, REML=TRUE)

I am trying to find a way to determine confidence intervals for each level of my animal and year random effects for the purposes of creating an effects plot to depict how these effects change over years and differ between animals.

When I run ranef(model), I am able to find the mean effect value for each level of the random effects, and I have been able to find the overall standard deviation and variance for the entire random effect using VarCorr(model).

I have seen ways to extract CIs for each level of a RE in an lme model, but am not finding any way to extract them for lmer models.

Any help?

With the caveat that the intervals you get are not technically confidence intervals (i.e., random effects values are not 'parameters', you can't use the CIs to test hypotheses about the group-level values), there are several built-in tools:

library(lme4)
example(lmer)
rr <- ranef(fm1)
class(rr) ## ranef.mer
methods(class = "ranef.mer")
## [1] as.data.frame plot          print

as.data.frame gives you the conditional standard deviations, which you can use however you want (e.g. compute condval ± 1.96*condsd):

head(as.data.frame(rr), 3)
##    grpvar        term grp    condval   condsd
## 1 Subject (Intercept) 308   2.258551 12.07086
## 2 Subject (Intercept) 309 -40.398738 12.07086
## 3 Subject (Intercept) 310 -38.960409 12.07086

This gives you a built-in effects plot:

lattice::dotplot(rr)

(lattice is a built-in R package, you don't need to install it separately)

You can also get random effects values and confidence intervals with the broom.mixed package:

broom.mixed::tidy(fm1, effects = "ran_vals", conf.int = TRUE)

You can also get effects plot for the random effects terms, and information about the conditional standard deviations, from the ggeffects package and/or the sjPlot package ...

These plots are also often called 'caterpillar plots', so searching SO with this keyword will get you more information ...