0

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

Dimitris Rizopoulos
  • 17,519
  • 2
  • 16
  • 37
blokeman
  • 136
  • 1
  • 8

1 Answers1

3

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

Dimitris Rizopoulos
  • 17,519
  • 2
  • 16
  • 37
  • 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: https://drizopoulos.github.io/GLMMadaptive/articles/GLMMadaptive_basics.html#mixed-effects-logistic-regression – Dimitris Rizopoulos Dec 14 '18 at 16:26
  • Sir, in your original answer you explained that the logistic (non-Laplace) GLMM's likelihood function is different from that of a logistic GLM (which simply uses the binomial PMF). Now it seems to me that you contradict your original answer when you say there's "no theoretical reason" disallowing the comparison of a binomial (fixed-effects) log likelihood to the likelihood of a logistic non-Laplace GLMM. My brain melts. I mean OK, any two things can be compared. The question is: is it helpful or informative to compare those two? And your first post seems to suggest that it it's not. – blokeman Dec 17 '18 at 06:46
  • To compare two models with a likelihood ratio test you need the two models to be nested, i.e., the one being a special case of the other. The simple logistic regression model is a special case of the mixed-effects logistic regression model if the variance of the random effects is set to 0. Hence, the two are comparable with a LRT. The LRT does **not** require that both likelihoods are written in the same manner (i.e., indeed the mixed-effects lositic regression likelihood cannot be written in terms of the binomial PMF). – Dimitris Rizopoulos Dec 17 '18 at 07:17
  • Yes sir. Two nested logistic models, the only difference being the random intercept. Therefore the LRT compares -2(LL(fixed.mod) - LL(mixed.mod)) to a chi-squared distribution with df = 1, no? But what does it use for LL(mixed.mod)? Does it use `sum(dbinom(data$y, size = 1, prob = fitted(mixed.mod), log = TRUE))` or does it use the LL that was used to arrive at the parameter estimates (whose calculation requires an integral)? The two are different. If it uses the latter and that's fine, why does `anova()` refuse to do so whenever >1 quadrature points were used? I assume there's a good reason. – blokeman Dec 17 '18 at 09:04
  • You can do the LRT with more than one quadrature points using the GLMMadaptive package - check the link provided above. – Dimitris Rizopoulos Dec 17 '18 at 09:23
  • You've been very helpful, sir. I just gave GLMMadaptive a shot, and when trying to fit an equivalent of my lme4 logistic GLMM with nAGQ = 20 (1 random intercept), I get an error message stating that 'a large coefficient value has been detected during the optimization'. But the largest fixed-effect coefficient in the lme4 model is -4.1, i.e. well within the bounds of the default limit of 10, so I'm not sure how to proceed. On a second note, can't we also perform the LRT manually by calling `pchisq(fixed.model$deviance - -2*logLik(mixed.model), df = 1, lower.tail = F) / 2` ? – blokeman Dec 18 '18 at 16:38
  • Try adding the option iter_EM = 0 – Dimitris Rizopoulos Dec 19 '18 at 10:54