Calculating by-group estimates and CIs for a random intercept grouping variable in `lmer()`

Hi there,

I am trying to extract by-group marginal means and their CIs of an lmer() model with group as a random effect. For example:

library(lme4)
library(tidyverse)
library(parameters)

data <- read_csv("https://vincentarelbundock.github.io/Rdatasets/csv/palmerpenguins/penguins.csv") %>%
  drop_na()

model <- lmer(flipper_length_mm ~ bill_depth_mm + (bill_depth_mm | species),
              data = data)

I can access the bill_depth_mm estimate for each species by using coef(model), but I was wondering if it is possible to get the CIs with parameters. Reading the tutorials, I realized that I can get SEs using:

standard_error(model, effects= "random")

I was wondering if these are SEs for the by-group estimates and if calculating CIs manually makes sense here. Is there any function in parameters that I am missing for calculating CIs?

Start with the reproducible lmer() example in the issue and compare coef(model) with standard_error(model, effects = "random"). Determine whether those results represent by-group estimates and whether parameters has an existing confidence-interval entry point for them. Done means the supported behavior and CI approach are established and documented or implemented, with the example used for verification.