I have 3 level nested data that is causing convergence problems (glmer function in lme4, R for multilevel logistic regression).

The Data

Country(4 groups) -> School(100 groups) -> respondent
A respondent is nested within school, and the school (id number) is nested within country. Each school id.number is only within one country. So school-1 is only in Country-A.
Data is c.a. 8000 rows. Final model has 10 variables (three 2-level factors, 6 covariates + Country)

Cntry-A    1       m       0
Cntry-A    1       f       0
Cntry-A    2       m       1
Cntry-B    10      f       0
Cntry-B    10      f       1
Cntry-B    11      m       1
Cntry-C    100     f       0
Cntry-C    100     f       1
...        ...     ...     ...

The problem

  • Country has too few levels, four is not enough to have as a random intercept. But I wonder if it is okay to have so few groups for level-3 if level-2 has so many groups, because School is nested within Country?
  • Because I am confident that four levels is not enough, and because I am actually interested in the country variable, I have included it as a fixed effect rather than random.
  • I hypothesize that the effect of covariate1 on response.var is different within each country.

The model

     glmer(response ~ 
             covar1*country +
             gender +
             ... +

The question

  • Why am I getting convergence warning for the final model? Can I trust the results?
    • I have read a lot about this. Using more iterations and other optimizers sometimes fixes the problem. And using optimx fixes the problem as well. Running the model with many differnt optimizers (as suggested by lme4 auth) shows nearly identical results, so according to them I can ignore the convergence warnings.
    • But I would like to use multiple imputation for the data. With imputed data I cant use optimx.
    • I have read through this article https://rstudio-pubs-static.s3.amazonaws.com/33653_57fc7b8e5d484c909b615d8633c01d51.html and checked off everything, singularity etc.
    • I use grand mean centering for covariates.
  • Am I defining the model correctly? Adding Country (interaction or just as main effect) creates convergence issues.
    • Is it because the model has 100+ intercepts, one for each school. And a school is only within one country?

1 Answer 1


The default estimation procedure of glmer() is a Laplace approximation that is known not to be optimal for binary data. You could try using adaptive Gaussian quadrature by adjusting the nAGQ argument.

In addition, you could also have a look at the GLMMadaptive package that implements the adaptive Gausian quadrature rule in more general settings.

  • $\begingroup$ Thank you for the suggestion of nAGQ argument, I was unaware of this. I'll see if it fixes the problem. I think the GLMMadaptive actually solved my problem. Super fast and no warnings while I got convergence warnings for same model in lme4 with nAGQ=20. Thanks a lot! Are there any major differences between Laplace approx and adaptive Gaussian q that affect the results, or something I should be aware of? Or can I treat the results the same way from lme4 and GLMMadaptive? $\endgroup$ Sep 4, 2018 at 19:28
  • 1
    $\begingroup$ The difference is in the approximation of the integrals over the random effects in the definition of the marginal log-likelihood. Laplace does a normal approximation to the integrand, whereas the adaptive Gaussian quadrature (AGQ) approximates the integrals by a finite sum. The advantage of the AGQ is that you can control the error of the approximation by increasing the number of quadrature points (though at the expense of computational time). If you're interested, you can find more info in my course notes available at drizopoulos.com/courses/EMC/CE08.pdf (Section 5.3). $\endgroup$ Sep 4, 2018 at 19:40
  • $\begingroup$ Thanks a lot for your help. The notes will definitely come in handy. $\endgroup$ Sep 4, 2018 at 20:30
  • $\begingroup$ Hey, sorry to bother you again. I've been using your package (excellent package btw, solved most of my model problems so far). But is there a way for me to define a Null-model and get Std.Error and OddsRatio for the Null-model? I defined it as: mixed_model(fixed = response ~ 1, random = ~ 1|country, data = Dat, family = binomial()) $\endgroup$ Sep 23, 2018 at 21:18
  • $\begingroup$ That would be a null model for the mean part, but it still assumes that outcome measurements within country are correlated. $\endgroup$ Sep 24, 2018 at 6:42

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.