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()?


1 Answer 1


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,$$


  • $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)).

  • $\begingroup$ 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? $\endgroup$
    – blokeman
    Dec 11, 2018 at 16:08
  • $\begingroup$ 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. $\endgroup$ Dec 11, 2018 at 16:13
  • $\begingroup$ 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 $\endgroup$
    – blokeman
    Dec 11, 2018 at 16:18
  • $\begingroup$ 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. $\endgroup$ Dec 11, 2018 at 20:58
  • $\begingroup$ 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/… $\endgroup$ Dec 14, 2018 at 16:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.