r marginal effects plot for multiple objects

My question is how to plot marginal effects or adjusted predictions using ggeffects() or ggemmeans() when I run my model on multiple subsets of data using dlply which creates summary estimates as multiple objects

First, I start with a toy dataset.

   library(lme4)
   data(Soybean)
   library(plyr)
   library(nlme)
   library(ggplot2)
   library(ggeffects)

head(Soybean)

       Plot   Variety Year Time weight  
1988.1 1988F1       F 1988   14  0.106 
1988.2 1988F1       F 1988   21  0.261 
1988.3 1988F1       F 1988   28  0.666  
1988.4 1988F1       F 1988   35  2.110  
1988.5 1988F1       F 1988   42  3.560  
1988.6 1988F1       F 1988   49  6.230  

Second, I fit model with Variety as random intercept on each subset of dataset by Year.

table(Soybean$Year)

1988 1989 1990 
 156  128  128 
 
table(Soybean$Variety)

  F   P 
204 208 

# Fit LME model for each subset(Soybean)
MxM1 <- dlply(Soybean,
                    "Year",
                    function(x)
                      lme(fixed = weight ~ Time,
                          random = ~ 1 |Variety ,
                           data=x, na.action=na.exclude))

These are the coefficients of Random intercept (Variety) from the model.

2 varieties, 3 years , 6 different intercepts,

  ldply(MxM1, coef)
  Year (Intercept)  Time
1 1988       -8.65 0.349
2 1988       -7.73 0.349
3 1989       -6.71 0.242
4 1989       -4.04 0.242
5 1990       -7.28 0.309
6 1990       -6.63 0.309

How do I plot the marginal effects or adjusted predictions using ggeffects or ggemeans for lme estimates from multiple datasets ?

Expecting a marginal effects or conditional effects plot like this.

Is something like this what you're looking for? Here's an example with the Chile data from carData.

library(ggplot2)
library(lme4)
library(dplyr)
library(ggeffects)
data("Chile", package="carData")
Chile$votey = ifelse(Chile$vote == "Y", 1, 0)

Run the models as you suggest using dlply().

mods <- plyr::dlply(Chile, "sex", function(x)glmer(votey ~ age + education + income + (1|region), 
                                           data=x,
                                           family=binomial))

Use map() from purrr to generate predictions using ggpredict() and then turn the results into raw data frames. (you could just as easily use lapply().) Then, you bind the predictions together into a single data frame.

preds <- purrr::map(mods, function(x)ggpredict(x, c("age[all]", "region"), type="re")  %>% as.data.frame())
preds <- bind_rows(preds, .id="sex")

Then, make the plot:

ggplot(preds, aes(x=x, y=predicted, ymin=conf.low, ymax=conf.high, color=group, fill=group)) + 
  #geom_ribbon(col="transparent", alpha=.2) + 
  geom_line() + 
  facet_grid(sex ~ .)

<!-- -->

^{Created on 2023-03-20 with reprex v2.0.2}