I stumbled upon this paper yesterday, https://doi.org/10.1002/sim.8782 (Gory, J. J., Craigmile, P. F., & MacEachern, S. N. (2020). A class of generalized linear mixed models adjusted for marginal interpretability. Statistics in Medicine .) where some interesting results are derived for marginally interpretable GLMMs, or rather results on how to model parameters with some different commonly used link functions with marginally interpretable “fixed effects” while still allowing for “random effects”.
The idea is that (a bit more general than in the paper) the parameter is modeled by the restriction:
h(\mathbf{x}_i^T\mathbf{\beta})=\int h(\mathbf{x}_i^T\mathbf{\beta} + \mathbf{d}_i^T\mathbf{u}+\mathbf{d}_i^T\mathbf{a}_i)f_\mathbf{U}(\mathbf{u})d\mathbf{u}, where h() is the inverse link-function, \mathbf{u} is the random variables for the random effects, and \mathbf{a}_i is an adjustment factor for a specific \mathbf{d}, i.e. the incidence matrix for the random effects.
The paper derives some restrictions on the distribution on the random effects for some different link functions given that there exists an adjustment, e.g. that the MGF must exist for the random effect distribtuion in a log-link. The most permissive type of link functions are those that translate the real line onto a strictly bounded domain (e.g. logit, probit, cloglog, etc., typically inverses of CDFs of other distributions) where the adjustment exists for any choice of random effect distribution.
A result that is relevant for computation is the proposition 8, that we can do a change in variable in the integral of the multivariate normal distribution, making it a univariate integral. I did not know this, but I guess it is a standard result of multivariate normals.
I don’t know enough Stan to program a model putting it to use, at least not yet, but I hope I have sparked your interest in a problem I think is very interesting.