This should be the last post in this lengthy thread. There was a conceptual problem in my initial approach.
Averaging over fitted values incorporating random effects does not appropriately integrate over the random effects. Specifically, you do not get an estimate of the average marginal effect (AME) in the population that was sample from. You get an estimate in the people sampled.
If the goal is an estimator of the AME in the population, the random effects need to be integrated out.
Averaging over fitted() estimates, whether or not random effects were used in these predictions, is not the same.
I have implemented this in an R package, brmsmargins. The vignette for mixed effects models
describes and outlines this more:
Marginal Effects for Mixed Effects Models • brmsmargins (joshuawiley.com)
I have done preliminary testing and some simulations, which all seems to be working.
I would really welcome anyone interested to test and provide feedback, log issues on github, etc.
@jackbailey I know you had interest in this space so tagging you here.