Calculating Log-Likelihood of Logistic Adaptive-Quadrature GLMM for Comparison with Fixed Model

Fitting a binary logistic GLMM here, with ungrouped data (all responses either 0 or 1).

It says in this thread and in the documentation of anova.merMod that the significance of the random intercept cannot be tested with a LRT against the fixed model if the random-intercept model was fit using adaptive quadrature with more than 1 quadrature point. More specifically, the documentation says that the output of logLik() for a GLMM is 'only proportional to the log probability density of the response variable'. Here's a silly question coming from a non-statistician:

Why not just calculate the log-likelihood for the random-intercept model manually, using the response vector and the fitted values, and compare that to the corresponding log-likelihood of the fixed-effects model? Like so:

LL.mixed <- sum(dbinom(data$$response, size = 1, prob = fitted(mixed.model), log = TRUE)) LL.fixed <- sum(dbinom(data$$response, size = 1, prob = fitted(fixed.model), log = TRUE))
pchisq((-2*LL.fixed)-(-2*LL.mixed), df = 1, lower.tail = FALSE) / 2

It seems straightforward enough, since the model's log-likelihood is just a sum over the log probability densities of every data point given their respective fitted values, no? Why should the number of quadrature points used in estimating the random effect matter at all? It doesn't enter in any way into this simple calculation of the model's log likelihood, so what's the complication?

I now notice that while the above calculation using sum() and dbinom() replicates exactly the output of logLik() for any logistic GLM, these two outputs do not match for GLMMs, not even if the number of quadrature points is 1 (using glmer()). What in the world is going on? Is not the log-likelihood of any logistic model, mixed or fixed, based quite simply on the responses, the fitted values, and the binomial PMF? Is it me who is f**ing up, or is it logLik()?

The log-likelihood function of a mixed-effects logistic regression model with normal random effects does not have a closed-form. In particular, it is defined as:

$$\ell(\theta) = \sum_{i = 1}^n \log \int p(y_i \mid b_i; \theta) \, p(b_i; \theta) \; db_i,$$

where

• $$y_i$$ is the multivariate response vector for the $$i$$-th subject.

• $$b_i$$ are the random effects, which have a normal distribution with mean zero and covariance matrix $$D$$.

• $$\theta$$ is the parameters vector that includes the fixed effects regression coefficients, and the unique elements of the $$D$$ matrix.

• $$p(y_i \mid b_i; \theta)$$ is the binomial probability mass function (i.e., what you get from dbinom()) conditional on the random effects.

• $$p(b_i; \theta)$$ is the probability density function of the multivariate normal distribution for the random effects.

The integral in the definition of the log-likelihood above does not have a closed-form solution and needs to be approximated with numerical methods, such as the adaptive Gaussian quadrature. Hence, the log-likelihood of a mixed-effects logistic regression model cannot be calculated as sum(dbinom(..., log = TRUE)).

• Are you saying, sir, that if a logistic GLMM was fit using more than 1 quadrature point, its fit can no longer be checked using dbinom() because the binomial PMF no longer applies? If so, that's kind of hard to wrap one's head around. It entails that if I fit a standard logistic GLM to the same data and obtain the EXACT SAME fitted values as the GLMM, dbinom() is applicable for calculating the Log-Likelihood. Yet dbinom() doesn't work for a GLMM with the same fitted values? Same fitted values, same data, yet the binomial distribution works in one case but not in the other? – blokeman Dec 11 '18 at 16:08
• Even with 1 quadrature point the likelihoods are not the same. They are only intrinsically the same if the variance of the random effects is zero. – Dimitris Rizopoulos Dec 11 '18 at 16:13
• Right, but with 1 quadrature point the likelihoods are comparable, right? That is presumably why anova() allows comparing a Laplace-approximated GLMM with one random intercept with an fixed-effects model without the random intercept. By contrast, if more than one quadrature point was used, there is only the error message Error in anova.merMod(mix4, mod4b) : GLMMs with nAGQ>1 have log-likelihoods incommensurate with glm() objects – blokeman Dec 11 '18 at 16:18
• AFAIK there is no theoretical reason that disallows the comparison of the log-likelihoods between a logistic regression and a mixed effects logistic regression fitted with more than 1 quadrature points. I'd say this is a feature of the lme4 package. – Dimitris Rizopoulos Dec 11 '18 at 20:58
• As an alternative to the lme4 package you can try the GLMMadaptive package that implements the adaptive Gaussian quadrature and can do this test. For an example, check: drizopoulos.github.io/GLMMadaptive/articles/… – Dimitris Rizopoulos Dec 14 '18 at 16:26