Question: Is there a citation for the 'bias adjustment' used for the outcome scale predictions of a model with random effects

[statsccpr]

Hi thank you for the package and hard work, it seems flexible and helpful post model fitting.

I've been experimenting with the regrid methods. Looking at the bias adjusted incidence probabilities for the random effect logistic regression example.

https://cran.r-project.org/web/packages/emmeans/vignettes/sophisticated.html#mcmc

Will this 'marginalize/integrate out' the random effects, then return the incidence probabilities? Along the lines of https://drizopoulos.github.io/GLMMadaptive/reference/marginal_coefs.html

"Hedeker et al. (2017) to calculate marginal coefficients from mixed models with nonlinear link functions"

Is it correct to say, using emmeans + regrid + bias adjustment + sd is similar to the marginal_coefs() of GLMMadaptive where the difference is emmeans uses the taylor approximation for the adjustment while glmm adaptive is using montecarlo integration

Is this the adjustment emmeans is using even for the random effects ('CBPP' example)? The vignette below uses the 'pigs' example, (no random effects) https://cran.r-project.org/web/packages/emmeans/vignettes/transformations.html#bias-adj

I'm guessing this is the implementation talked about here

https://github.com/rvlenth/emmeans/issues/99#issuecomment-512948053 https://github.com/rvlenth/emmeans/blob/master/R/emmGrid-methods.R#L490

Thank you

[rvlenth]

The answer to the first question is yes or no, depending on whether GLMMadaptive::marginal_coefs returns comparable results to those from stan_glmer (the $stanfit slot).

It's important to understand that emmeans functions basically analyze the model -- not the data. They take whatever fixed effects are handed to it in the model, and whatever covariance matrix thereof. The bias adjustment is strictly interpreted in terms of trying to estimate E(Y) = E(h(Z)) where Y is the response, Z is the transformed response, and h is the inverse transformation.

The other thing to know is that bias adjustment should be done according to the actual random effects, as opposed to the model deviance. For example, the pigs illustration has a response transformation with random error and we use sigma(model) for adjusting bias. If we were to fit, say, a Poisson regression model using glm with family = poisson, we would not do any bias adjustment because the GLM is already framed in terms of the response mean. in a GLMM like that cbpp example, note that the bias adjustment is based on the variance of the random effects in the model, and nothing from the deviance of the GLM part of the model.

So it is important to know when to use the bias adjustment, and what variances should and should not be used. But the code for actually doing the adjustment is the same second-order approximation. The adjustment is just a second-order Taylor approximation, h(Z) ~= h(mu) + h'(mu)(Z - mu) + 0.5 h''(mu)(Z - mu)^2, so that E[h(Z)] ~= h(mu) + 0 + 0.5 h''(mu)*Var(Z)

I am aware that there are other ways of doing this. Another place I have seen this done is in the MCMCglmm package, where they use an adjustment before back-transforming whereas I back-transform then adjust with that h'' term.

I hope this helps clarify what is done.

[rvlenth]

I am closing this issue; but you can still comment if you have further questions.

[statsccpr]

Yes thank you. Your response clears things up greatly.

This is more of an aside, but would you have any insight on using the taylor approximation for a glmm, where the covariance of the MVN random effect is something more complex than sigma^2*Identity ? Say, a covariance matrix with non 0 off diagonals (so the random effects are correlated)? So far, I think it's doable and would follow along the lines of your examples.

Using the mcmc pigs example,

cbpp.rstan <- rstanarm::stan_glmer(
    cbind(incidence, size - incidence) ~ period + (1|herd) + (1|unit),
    data = cbpp, family = binomial,
    chains = 2, cores = 1, seed = 12345, iter = 500)
cbpp.rg <- ref_grid(cbpp.rstan)
cbpp.sigma = as.matrix(cbpp.rstan$stanfit)[, 78:79]
## iterations Sigma[unit:(Intercept),(Intercept)] Sigma[herd:(Intercept),(Intercept)]

totSD <- sqrt(apply(cbpp.sigma^2, 1, sum))
cbpp.rgrd <- regrid(cbpp.rg, bias.adjust = TRUE, sigma = totSD)

I understand that the key here is to get the posterior standard deviations of the unit random effects (and combine it with the std dev of the herd random effects)

If instead, the distribution of the random intercepts for the 'unit' effects was MVN(0,Sigma) Sigma = matrix(var1,cov12,cov21,var2) where lets say there's only 2 unit levels for ease of example.

My thinking is we would additionally need to account for the correlation induced by cov12 and cov21, where the variance across the unit random effects 'theta_unit' would be what's used for the taylor adjustment.

var_unit_correlated = Var(theta_unit[1], ... ,theta_unit[g])

to plug into

var_herd = cbpp.sigma[,2]^2  # just keep herd variance

total.SD = sqrt(var_herd  + var_unit_correlated )
summary(emm, type = "response", bias.adjust = TRUE, sigma = total.SD)

Or is something more complex needed that involves theta's entire Var-Cov matrix, Sigma

[rvlenth]

Yes, if two random effects are correlated, then for the total variance you need to add twice the covariance. In general, you need (in matrix form) SD_tot = 1' V 1 where 1 is a vector of all 1s and V is the covariance matrix.

[statsccpr]

Hi Russell, following up on this last thought (although I'm in no rush to use it).

After sleeping on it and digging around I think there needs to be a small tweak on the emmeans interface to properly do this (multivariate) taylor approximation for correlated random effects.

Without further modification, I think using the current interface regrid(..., bias.adjust = TRUE, sigma = totSD)

is only valid for the simple random intercept glmms, where the random intercepts have the exchangeable covariance. MVN(0, constant * Identity_nxn)

I think for any other more exotic covariances, say in (as much as i don't like sas, they have pretty good docs) https://support.sas.com/resources/papers/proceedings/proceedings/sugi30/198-30.pdf

we would need the multivariate taylor approx for the adjustment

https://stats.stackexchange.com/questions/401749/multivariate-taylor-series-for-moments-of-a-random-variable

The user (or maybe the emmeans package internally) would need to compute 1/2Tr(f′′(E[λ])Var(λ)) and pass it into the internal calculations in your helper functions below

https://github.com/rvlenth/emmeans/blob/5b02798ad96c951a0fcef1faa58b1c762ed9e718/R/summary.R#L772

I think the current framework is only setup to handle if 1/2Tr(f′′(E[λ])Var(λ)) reduces to the simple scenario where you can just plug in a constant to sigma in regrid(..., bias.adjust = TRUE, sigma = totSD)

[rvlenth]

You are correct in your assessment of the method currently used, as well as its limitations. However, expanding it to support a multivariate Taylor approximation is not what I would call a "small tweak."