Search code examples
lme4mixed-modelsrandom-effects

"iteration limit reached" in lme4 GLMM - what does it mean?


I constructed several glmer.nb models with different combinations of random intercepts, and for one of the models (nested random intercepts, with the lowest AICc), I consistently get: "iteration limit reached", without the usual "Warning message: In theta.ml(Y, mu, weights = object@resp$weights, limit = limit, :..."

Here's what I know:

  1. it is a warning (from the color) but not labeled as such
  2. you can also have that warning with GLMs and LMERs

Here's what I don't know:

  1. does it mean the model is invalid?
  2. what causes that issue?
  3. what could I do to resolve that issue?

Here's what I searched:

  1. https://stats.stackexchange.com/questions/67287/very-large-theta-values-using-glm-nb-in-r-alternative-approaches (no explanation as to the why and how)
  2. GLMM FAQ: no mention
  3. I am not the only regularly running into that or similar problems: Using glmer.nb(), the error message:(maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate is returned https://stats.stackexchange.com/questions/40647/lme-error-iteration-limit-reached/40664

Here's what would be highly appreciated:

A more informative warning message: did the model converge? what caused this? What can one do to fix it? Can we read more about this (link to GLMM FAQ - brms-style)?

This is a general question. I did not provide reproducible code because an answer that is generalisable would be most useful.


Solution

  • library(lme4)
    dd <- data.frame(f = factor(rep(1:20, each = 20)))
    dd$y <- simulate(~ 1 + (1|f), family = "poisson",
                     newdata = dd,
                     newparam = list(beta = 1, theta = 1),
                     seed = 101)[[1]]
    
    m1 <- glmer.nb(y ~ 1 + (1|f), data = dd)
    

    Warning message: In theta.ml(Y, mu, weights = object@resp$weights, limit = limit, : iteration limit reached

    It's a bit hard to tell, but this warning occurs in MASS::theta.ml(), which is called to get an initial estimate of the dispersion parameter. (If you set options(error = recover, warn = 2), warnings will be converted to errors and errors will dump you into a debugger, where you can see the sequence of calls that were active when the warning/error occurred).

    This generally occurs when the data (specifically, the conditional distribution of the data) is actually equidispersed (variance == mean) or underdispersed (i.e. variance < mean), which can't be achieved by a negative binomial distribution. If you run getME(m1, "glmer.nb.theta") you'll generally get a very large value (in this case it's 62376), which indicates where the optimizer gave up while it was trying to send the dispersion parameter to infinity.

    You can:

    • ignore the warning (the negative binomial isn't a good choice, but the model is effectively converging to a Poisson solution anyway).
    • revert to a Poisson model (the CV question you link to does say "a Poisson model might be a better choice")
    • People often worry less about underdispersion than overdispersion (because underdispersion makes results of a Poisson model conservative), but if you want to take underdispersion into account you can fit your model with a conditional distribution that allows underdispersion as well as overdispersion (not directly possible within lme4, but see here)

    PS the "iteration limit reached without convergence" warning in one of your linked answers, from nlminb within lme, is a completely different issue (except that both situations involve some form of iterative solution scheme with a set maximum number of iterations ...)