VA underestimates effects in binary models

[BertvanderVeen]

After implementing glmmVA() and a VA-logit formulation, I have repeated the simulations by Ormerod and Wand (2012) to compare performance of our implementation against their GVA with adaptive-GH quadrature.
There are 3 sets of simulations, each with 2 trials that determine the cluster size of the random intercept. The following figure warranted investigation of our implementation:

As is clear in set 2 and set 3 (the logistic model), our implementation severely underestimates the standard deviation of the random effect. Note,the VA implementation of Ormerod and Wand, which combines logistic regression with AGH-quadrature, does not have same issue.

This is something I have encountered before in probit models in gllvm (diagnosed as estimated LVs taking very small values), but attributed in the usual "VA underestimates variance" way. However, since Ormerod & Wand's GVA implementation does not have the same issue, something else entirely seems to be going on.

The following application to simulated data from set 3 trial 1 (with true parameters beta0 = 0, beta1 = 5, sigma = sqrt(1.5), 15 levels in the random effect and 8 data points per level) demonstrates the issue further:

data.csv

glmmVA.fit <- glmmVA(formula=vy~vx + (1|id), data = data, family = binomial)
glmmEVA.fit <- glmmVA(formula=vy~vx + (1|id), data = data, family = binomial, method = "EVA")
glmmLA.fit <- glmmVA(formula=vy~vx + (1|id), data = data, family = binomial, method = "LA")

where glmmVA.fit produces a very different result to EVA and LA. The model produces the following likelihood profile (which looks slightly different depending on the numerical details of the implementation):

plot(TMB::tmbprofile(glmmVA.fit$TMBfn, "log_sigma"))

Clearly, negative values for log(sigma) receive a better negative log-likelihood. In contrast to EVA and LA on the same data, the VA objective does not seem to have a unique minimum for sigma. This seems odd, given that EVA is supposed to be a less accurate version of VA.

par(mfrow=c(1,2))
plot(TMB::tmbprofile(glmmEVA.fit$TMBfn,"log_sigma", parm.range = c(-10,10)), main = "EVA")
plot(TMB::tmbprofile(glmmLA.fit$TMBfn,"log_sigma", parm.range = c(-10,10)), main = "LA")

Investigating this fit with tmbstan in the following manner, produces the following posterior for the log standard deviation:

library(tmbstan)
m <- tmbstan(glmmVA.fit$TMBfn, seed = 101, iter = 10e3)
hist(m@sim$samples[[1]]$log_sigma, col="transparent")
hist(m@sim$samples[[2]]$log_sigma, col="transparent", border = "blue", add = TRUE, lty = "dotted")
hist(m@sim$samples[[3]]$log_sigma, col="transparent", border = "red", add = TRUE, lty = "dotdash")
hist(m@sim$samples[[4]]$log_sigma, col="transparent", border = "orange", add = TRUE, lty = "longdash")
traceplot(m, "log_sigma")

So that it is clear the model formulation is not the issue: MCMC shows the posterior to be centered far away from zero, and it does not even have a mode at zero as one might expect. At this point, it is not clear to me what is going on, or how to fix it.

[BertvanderVeen]

Bumping up the number of observations per level of the random effect, from 8 to 20, does better. So, data deficiency definitely triggers this. That is not to say this is not an issue, more that it is such a (relatively) minor issue that more information in the likelihood improves it somewhat.

[BertvanderVeen]

Exploring also the LA fit by MCMC shows that, in the case of LA, the MCMC explores the region close to the maximum likelihood estimate well, so that the posterior mean for log(sigma) is -0.16, which corresponds much better to the MLE (0.06) than in the VA case.

In the VA case, the MCMC does not explore the negative log-likelihood at all in the region of negative log(sigma); i.e., the posterior mean does not correspond to the MLE at all, despite the fact that the sampler explores the exact same objective as the numerical optimisation. Weird.

[fhui28]

Thanks Bert,

You mentioned this in an email separately, but to clarify, this problem is:

• Not limited to logit link and also occurs when you set binomial(link = "probit") and use the latent normal parametrization?
• Not limited to gllmmVA and you also get this when row effects are included in binary gllvm fitted using VA. A deeper worry would be that if it occurs for variance components row effects, then perhaps it also occurs with latent variables in either the loadings (since this is reduced-rank random effects covariance matrix) and/or the scale of the latent variables when this is also estimated? But perhaps let's solve the GLMM problem first!

I also agree that it is not surprising the model formulation is not the problem, as the model after all is pretty standard. It's the VA polya-gamma set up of this model which is the problem (at least in the case of logit link), and suggests somehow introducing the latent variational distribution with these extra parameters (although there are no extra parameters though cause everything is "conjugate"?) is the problem.

Hmmm...might need to get the gllvm Jedi Council together on a Zoom meeting to stare the maths and brainstorm?

[BertvanderVeen]

Yes, I suspect it also occurs for the loadings and for correlation parameters of random effects more generally.

I will try to finish up the alternative logistic derivation that we discussed during your visit, and if that has the same problem I will plan a meeting.

Ultimately, in our VA formulation we decouple the variational distributions. Effectively, in models that include an auxiliary variable, we are overpaying the marginal posterior distribution of the random effects with a normal. In most cases we know that this is inaccurate and underestimates the posterior variance as it ignores heavy tails.

In the polya gamma case, the conditional posterior of the random effect given the polya gma random variable is normal. The marginal is (IIRC) a mixture of gaussians, which we are overlaying with a Gaussian.

So, perhaps this comes down to searching for binomial model formulations where the marginal posterior is as close to Gaussian as possible?

[BertvanderVeen]

@fhui28 I've shared the alternative approach with you on overleaf. The document is messy, but it will get the idea across.

The main point now is, I don't think that an alternative derivation can solve this. I could be wrong, but I think the issue is that the normal posterior is too simple; it doesn't capture the true shape of the posterior.

In any of the models where we introduce auxiliary latents, including the probit, logit, and negative-binomial models, we factorize the posteriors and eventually assume a marginal posterior for the random effect. The schemes we have used (including the cloglog, PolG and PoisG approaches) allow us to find a closed form approximation, but make no promises when it comes to the true shape of the posterior distribution. Especially in finite samples, that can be problematic, and probably worse when a hierarchical structure is present in the auxiliary variables (such as for the alternative logit derivation).

Assuming normality of the marginal posterior is probably fine when large sample asymptotics kick in, which clearly does not happen under the simulation shown here: a binary model with $n = 8$ observations per level of the random effect.

As I see it, there are a few ways to continue: 1) a finite sample correction (in the form of e.g., the Hessian of the joint likelihood similar to Laplace) for the variational covariance matrix to keep it from collapsing, or 2) a more complex shape for the approximate posterior (i.e., a finite mixture with a low number of components, but that very quickly blows up in the number of VA parameters), 3) use a different divergence, 4) ad-hoc scale the contribution of the likelihood and KL terms to the objective function.

Fortunately, there is a large body of literature to draw from on this from, but doing it well and in a manner that keeps the models fast needs some thinking.

It also explains why EVA might do better, as it includes local curvature information in the form of a Hessian.

[BertvanderVeen]

On that note, even adding a constant of (say) 1 to the diagonals of the variational covariance matrix already does better; the random effects take their normal sizes, same for the variance estimate of the random effects. It doesn't unbias the fixed effects, but it's something..

That confirms my suspicion that collapse of the variational distribution is the primary culprit.

[fhui28]

Thanks Bert, * Re Overleaf: on a very surface skim that seems correct on my end, but at the risk of being a PITA, if you please double check to see if you start at the beginning by assuming phi = 0 and phi = 1 [exponential distribution], then you produce the same cloglog variational log-likelihood you derived in the your other Overleaf, and the logit variational log-likelihood at the end of your current overleaf. There are a lot of mathematical moving parts here, so worth checking the maths multiple times; * So to confirm, your previous mention to me that our truncated-Poisson cloglog approach for VA works well was optimistic, and now it does not work well. I guess the way you checked this previous was with GLLVMs, and here you are checking this with GLMMs; The main point now is, I don't think that an alternative derivation can solve this. I could be wrong, but I think the issue is that the normal posterior is too simple; it doesn't capture the true shape of the posterior. I can imagine this being a problem, as it is a known phenomenon but not necessarily a problem for VA in general. But I also need much more clarification on this first: if indeed the combination of VA auxiliary latents + assuming q(LVs/random effects) really is a problem, then we are in deep crap because: 1) this is more than half [?!] the distributions we have in the current package; 2) we have missed a big problem in numerous previous papers we have published on this before, and we risk promoting/have promoted bad science. I know I am sounding alarmist, but your findings suggest we should be alarmed :😱 * I don't understand why we did not pick this up in all the previous VA projects we worked on before. Was it because we kept hand-waving the binary response cases/make the negative binomial cases not that severe? Was it because the number of sites and/or species was large enough [but we have to introduce auxiliary latents per site-species combination!]. * Put another way, were we always in the "large sample asymptotics kicked in" empirical settings we did before? * A binary model with $n = 8$ observations per level of the random effect is not that small in many applications, and we have to emphasize that what you are saying is the combination of VA auxiliary latents + assuming q(LVs/random effects) is the problem. For example, if we use EVA or approximated the integral using quadrature then assuming q(LVs/random effects) is fine! So I am alarmed but also do not want to overblow the problem yet until we know the exact combination of things that is sourcing the problem, let alone trying a elaborate solution e.g., I tried a marginal finite mixture for the VA distribution of LVs in a different factor analysis problem I did a few months ago with some colleagues on multivariate longitudinal data. The TL;DR is that it does tend to capture the true posterior shape better, but in terms of its overall mean and covariance it was not much different from what you get just assuming marginal normality on the LVs. Plus as you say, more VA parameters, and one usually would not advocate the mixture-specific means and covariances unless they meant something. Cheerio and a bit lost...more research and playing around needed!

…

________________________________ From: Bert ***@***.***> Sent: Tuesday, 19 August 2025 3:54 AM To: JenniNiku/gllvm ***@***.***> Cc: Francis KC Hui ***@***.***>; Mention ***@***.***> Subject: Re: [JenniNiku/gllvm] VA underestimates effects in binary models (Issue #237) [https://avatars.githubusercontent.com/u/19729157?s=20&v=4]BertvanderVeen left a comment (JenniNiku/gllvm#237)<#237 (comment)> On that note, even adding a constant of (say) 1 to the diagonals of the variational covariance matrix already does better; the random effects take their normal sizes, same for the variance estimate of the random effects. It doesn't unbias the fixed effects, but it's something.. — Reply to this email directly, view it on GitHub<#237 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABH6DBAZTAM6M5OS7BB7Q5L3OIHM3AVCNFSM6AAAAACD4E6KBGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTCOJXHA4DCNZVGQ>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[BertvanderVeen]

How mature is your mixture code? Could be worth finding a niche case where this happens and that your code can tackle, to see if it helps.

Wouldn't the mean of the mixture be interpreted similarly as our current VA means? I think the main issue is better capturing heavy tails in the posterior, not capturing multimodality, so there may be other tricks to try.

In short, yes I think we are often in a case where VA does well in a community ecological setting as the number of observations are given by the number of species, which is usually at least 30 or more.

I'm not completely sure how that translates to the site-species combination for the auxiliaries, to be fair. It might be that my conclusion about MF is off, and it is the information content in the likelihood that is a larger culprit, causing KL terms to dominate the objective over the likelihood, thus over-regularize the random effect estimates, and thus collapse the VA distribution.

Ultimately, I would say that collapse of the VA distribution is caused by a combination of factors, including the MF assumption, poor matching of the posterior approximation to the truth, and low information content in the likelihood. Thus there may be different solutions, each with pros and cons.

I wouldn't say we've been doing bad science. This has been (to a lesser degree) on my radar since the middle of my PhD, and the advice to people using the software when this happens has been to change to LA/EVA. However, a deeper investigation and mitigation of the issues for finite samples is worth the persuit, especially if we want the methodology to be adopted for fitting ordinary mixed-effects models.

Let's plan a meeting and discuss a bit more. I have some other things I need to do too before I am able to spend (much) more time on this.

[fhui28]

Thanks Bert, and apologies for sound alarmist. Your comments are reasonble.

But yes let's take it to the gllvm Jedi council and go from there!

[BertvanderVeen]

It might be interesting to augment the VA approach with q(u) ~MNIG(mu, Sigma, lambda), i.e., a multivariate normal-inverse gaussian (gaussian mean-variance mixture) where lambda controls heaviness of the tails, but the distribution is symmetric (a special case of MNIG for skewness parameter fixed to zero), while q(u) is then approximately normal for large lambda. The mean of the inverse gaussian needs to be a function of lambda (i.e., not fixed to 1), otherwise the first two moments of q(u) do not depend on lambda.

The entropy is, I believe, not quite available in closed form.

E(u) = mu, cov(u) = Sigma(1/1+lambda), with heavier tails than a Gaussian. MGF is available in closed form following a minor modification of the MGF of the multivariate normal.

In all fairness, this would need a lot more consideration, and the value of the log-likelihood was still about 7 points higher than EVA and LA.
It seems a bit suspicious to me at the moment, at the estimator for gamma only interacts with the variational covariance matrix.

[fhui28]

Thanks for trying Bert! Good to at least rule out "easy" solutions!

[BertvanderVeen]

The attached image shows the true posterior of the random effects for the simulation (straightforwardly calculated from the joint likelihood of the true parameters). It's safe to summarise that I'm not quite sure why VA doesn't capture this, and have not found any useful remedies.

The thing that comes to mind is the fact that there ar 15 levels in the random effect, but there are only 4 unique posteriors, so that we are effectively estimating the variance parameter with 4 levels rather than 15.
The left-most mode represents 1 level, the next 3 levels, the next 6 levels and the right-most 5 levels.

The posteriors are not particularly skewed, bimodal, or heavy-tailed, so the VA distribution doesn't seem like the issue. The only thing I can come to conclude is that KL is oversmoothing due to the likelihood not providing a particularly large amount of information. Not sure why that isn't an issue with AGH in Ormerod and Wand.

The other thing that comes to mind, in this setting with b0 = 0, b1 = 5, a large amount of the estimated p's are close to 1. That is bound to give some issues numerically.

I think I've exhausted most of my possibilities for pursuing this at present, so hopefully someone else can make a smart suggestion during Wednesday's meeting.

[fhui28]

I thought we established that the problem was not so much KL, but specifically the introduction of VA auxiliary latents, irrespective of the distribution e.g., problems in binary, negative binomial. This would explain why:

• AGH/EVA/LA does not have this problem because they do not [need to] introduce such latents in handling such problems;
• Why doing fancier q(LVs/random effects) does not necessarily ameliorate the problem.

The only persistent aspect appears to be that increasing the number of observations per random effect helps to resolve the problem. This suggests that as long as we have observations per VA auxiliary latent distribution, then we can survive this? Although, your above discussion confuses me because for things like probit binary and negbin we have as many VA auxiliary latents as observations, so then no matter what we do there we should always be screwed?

I don't have any super insightful thoughts at this point, but I am still at the point or wanting clarification as to what the exact source of the problem is!