convergence warnings in lme4

• Likelihood: the probability of the data given the model
• Given a model specification, find the values of the model parameters that maximizes the likelihood
• Find the version of specified model that best matches observed data

You need to remember that “likelihood” is a technical term. The likelihood of \( H \), \( Pr(O\|H) \), and the posterior probability of \( H \), \( Pr(H\|O) \), are different quantities and they can have different values. The likelihood of \( H \) is the probability that \( H \) confers on \( O \), not the probability that \( O \) confers on \( H \). Suppose you hear a noise coming from the attic of your house. You consider the hypothesis that there are gremlins up there bowling. The likelihood of this hypothesis is very high, since if there are gremlins bowling in the attic, there probably will be noise. But surely you don’t think that the noise makes it very probable that there are gremlins up there bowling. In this example, \( Pr(O\|H) \) is high and \( Pr(H\|O) \) is low. The gremlin hypothesis has a high likelihood (in the technical sense) but a low probability.

Sober, E. (2008). Evidence and Evolution: the Logic Behind the Science. Cambridge University Press.

• REML: Residualized Maximum Likelihood
  • Not actually the likelihood
  • Related to the \( n \) vs. \( n-1 \) in e.g. calculating SDs
    • If your sample sizes are small enough that this matters, then your sample sizes are too small.
• Deviance:
  • “\( -2 \log \text{Likelihood} \)”
  • (up to an additive constant, due to the difficultiy in defining the fully saturated model required for the deviance)
  • so maximizing the likelihood is the same as minimizing the deviance
  • the difference in the deviance is the \( \chi^2 \) value in the likelihood ratio test computed by anova()
    • differences on the log scale are the same as ratios on the original scale
    • the additive constant cancels out

• Once the REML/ML/deviance criterion is defined, we just need to find the value maximizes (REML,ML) / minimizes it (deviance).
• It turns out that this problem only depends on the random effects and not on the fixed effects!
• This is called function optimization.
• There are many techniques for this and many different implementations of each technique.
• There are no free lunches and no optimizer which will always work the best.
• This is why lme4 allows you to pick your optimizer and change its settings.
• Changes in the default optimization settings are often the cause of new warnings when updating your lme4 version.
• For more information, see:
  • ?convergence
  • ?lmerControl
• If there are concrete ways to improve the documentation, let us know!

m2 <- lmer(Reaction ~ 1 + Days + (1 + Days | Subject), data=sleepstudy, REML=FALSE)

anova(m,m2)

Data: sleepstudy
Models:
m: Reaction ~ 1 + Days + (1 | Subject)
m2: Reaction ~ 1 + Days + (1 + Days | Subject)
   Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
m   4 1802.1 1814.8 -897.04   1794.1                             
m2  6 1763.9 1783.1 -875.97   1751.9 42.139      2  7.072e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(VarCorr(m2))

 Groups   Name        Std.Dev. Corr 
 Subject  (Intercept) 23.7798       
          Days         5.7168  0.081
 Residual             25.5919       

Partial image for intercept with constant (optimal) slope and correlation.

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: Reaction ~ 1 + Days + ((1 | Subject) + (0 + Days | Subject))
   Data: sleepstudy

     AIC      BIC   logLik deviance df.resid 
    1762     1778     -876     1752      175 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9536 -0.4673  0.0239  0.4625  5.1884 

Random effects:
 Groups    Name        Variance Std.Dev.
 Subject   (Intercept) 584.22   24.17   
 Subject.1 Days         33.64    5.80   
 Residual              653.11   25.56   
Number of obs: 180, groups:  Subject, 18

Fixed effects:
            Estimate Std. Error t value
(Intercept)  251.405      6.708  37.481
Days          10.467      1.519   6.889

Correlation of Fixed Effects:
     (Intr)
Days -0.194

• Optimizers have ways to detect when they have reached an optimum.
• If they didn't, they wouldn't know how to stop!
• However, optimization may be stopped prematurely for various reasons :
  • limit on number of function evaluations (maxeval)
  • not enough improvement in deviance with additional steps (ftol_abs, ftol_rel)
  • not enough change in theta / RE estimates with additional steps (xtol_abs, xtol_rel)
• Trying to different options or even a different optimizer are good strategies.
• If an optimizer thinks that it has converged, then that's generally the best convergence diagnostic.
• There are a few post-hoc tests we can run, but they are not guaranteed to be accurate.

Option names are all for nloptwrap, the current default lmer() optimizer. glmer still uses bobyqa, which has different names ….

0.1 + 0.1 + 0.1 - 0.3

[1] 5.551115e-17

• Computers don't store infinite digits and internally use a form of scientific notation (floating point, 1e-16 => \( 1 \times 10^{-16} \)).
• Exact equality doesn't hold in the presence of rounding.
• Flat deviance curves may not be numerically flat with coarse approximations.
• Similarly, Numerically flat curves may not be actually flat .

• Try setting stricter convergence criteria and allowing more iterations (ftol_rel, xtol_rel, maxeval).
• Check your overall model fit by plotting fitted vs. observed, etc. to make sure your model isn't misspecified.
• Consider setting control=lmerControl(calc.derivs=FALSE).
  • The derivative check takes a loooooong time
  • And tends to deliver false positives with maximal models.
• BUT PAY ATTENTION TO YOUR MODEL FIT AND YOUR OPTIMIZER'S OWN WARNINGS!!!
• Note that warnings about convergence failure, e.g. a Hessian that is not positive definite, may still indicate genuine problems – or that you should at least consider rescaling, as the warnings tell you to do.

Successfully completed optimization may not be the global optimum, but only a local one.

adapted from https://stats.stackexchange.com/questions/384528/lme-and-lmer-giving-conflicting-results/

• For a 2x2 within-subject design with “maximal” random slopes, we have 10 parameters!
  • 4 standard deviations
    • 1 intercept
    • 2 main effects
    • 1 interaction
  • 6 correlations
    • (4 SDs * 3 other SDs) / 2 because order doesn't matter
• With typical 20-30 participant designs, we only have 20-30 levels of the grouping variable
• So we only have 2-3 “samples” per parameter we're trying to estimate – expect noisy estimates and overfitting!
• For a 2x2x2 design, we have more parameters (34) than samples!
• (The same math holds for item RE.)
• It is not the total number of observations that matters for RE but rather the number of levels of the grouping variable!

• One of the innovations of lme4 compared to its predecessor nlme is its ability to fit models with optima at the boundaries, e.g. when a random effect goes to zero.
• Note that this does not mean that there is no between participant/group variation, but rather that there is no variation not captured by the residual variation.
• This lack of variation beyond the residual variation may just be due to lack of power to detect it!

• Numerically, boundary fits have non-invertible, i.e. singular matrices.
• One or more of the REs being numerically zero is one possibility; however, not the only one.
• If any linear combination of REs is zero, then this also yields a singular / boundary fit. For a geometric interpretation, think of a sloping ceiling – you may not be able to get to the edge of the floor nor the complete top of the building, but the combination of constraints leading to you hitting an edge.
• Singular fits are not problematic per se – they are mathematically well-defined and it is possible that the singular fit is actually the best description of the data.
• However, they are often indicative of overfitting, e.g. because you don't have enough data to distinguish group variation from residual variation.
• This overfitting, or equivalently, lack of power is a problem.
• Singular fits are also slow to compute. So simplify your model or gather more data.

• rePCA() performs PCA on the random-effects matrix.
• This is useful for determining the effective dimensionality of our RE – in other words, how many RE terms are supported by your data.
• For the sleepstudy, we could get 96% of the variance explained with one RE term!

summary(rePCA(m2))

$Subject
Importance of components:
                         [,1]    [,2]
Standard deviation     0.9294 0.22260
Proportion of Variance 0.9457 0.05425
Cumulative Proportion  0.9457 1.00000

• see Bates et al. (2015), “Parsimonious Mixed Models”.

• The successor to lme4 is MixedModels.jl.
• In addition to technical reasons related to the differences between R and Julia, there are a few computational innovations that are now in MixedModels.jl but not lme4.
• Not all models supported by lme4 are yet supported by MixedModels.jl, but those that are generally much faster.

• I am working on an R-package that will fit models in Julia and then turn them into lme4 fits …

• Some interesting reading on that note: https://nextjournal.com/dmbates/complexity-in-fitting-linear-mixed-models/