multilevel modeling with lmer(): understanding failure to converge in a toy example

Question

I am trying to get a deeper understanding of failures to converge in multilevel models that I estimate with lmer(). "Failure to converge" is vague; I want to be able to specify the problem that underpins these failures, and to express it numerically in terms of the data that I'm using. I am starting with little knowledge of the math that underpins estimation of these models.

Consider the following toy example:

library(lme4)
set.seed(1234)
person <- factor(rep(c("Alice", "Bob", "Catherine", "David"), each = 3))
female <- rep(1:0, each = 3)
myDF <- data.frame(
  y      = c(1:6, 10:12, 1:3),
  person = person,
  female = female)
lmer(y ~ (1 + female | person), data = myDF) 

Running that example generates these warnings:

Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

Although I've kept the dataset small (n = 12), the convergence warnings don't seem to be due to the small sample. There are more observations than parameters to estimate, and with a small tweak, I can generate the same warnings when (say) n = 1200, while holding constant the number of parameters.

That said, I have only a superficial understanding of what the warnings mean. Can the problems be expressed in terms of the data used in this example -- in terms of the values of y, the values of female, and so on?

There are many posts here about failures to converge when using lmer(). They are useful, but they generally emphasize programming strategies –– use a different optimizer, etc. –– rather than the intuition behind those strategies. Some posts are helpful for developing intuition: for example, "Model failed to converge" warning in lmer() and this non-StackExchange post. But I fear that even those posts don't help me to understand the problem above. I haven't found any sources that work through simple numerical examples like this one.

Answer 1

Building on @Emma Jean's response, Your created data frame is not multilevel in the sense that person (the grouping variable) and female are at the same level. In order for your data to accommodate a random slope, you would need to measure a variable that varied within person. It could be anything - mood or cortisol or pocket change - as long as it is measured repeatedly. So although you measure each of your persons on multiple occasions, person Alice always has a value of 1 on female, Bob always has a 0, Catherine always has 1, and David always a 0. In a multilevel model with such data, you can model the association between female and each person's y-intercept (essentially an average of their three observed y values). Below is your code with the addition of a time-varying variable (mood):

library(lme4)
set.seed(1234)
person <- factor(rep(c("Alice", "Bob", "Catherine", "David"), each = 3))
female <- rep(1:0, each = 3)myDF <- data.frame(
  y      = c(1:6, 10:12, 1:3),
  person = person,
  female = female,
  mood = c(4.5, 3.5, 4, 2.1, 3, 3.5, 1.6, 2.8, 4.3, 4.6, 4.2, 3.9))

If you look at the data frame, you will see that mood varies within persons, so it is time-varying. In the multilevel framework, such time-varying variables can be investigated as having associations with the outcome that are allowed to vary across groups (i.e., persons). I did not create the mood variable with careful consideration of building such unique associations, and not surprisingly, none were found and I get convergence warnings when trying to estimate the random slope for mood.

Another important feature of time-varying variables, is that you can calculate summary statistics on them that then become group (person) level variables. So with time-varying mood:

myDF$pmn.mood <- with(myDF, ave(mood, person, FUN=mean))
str(myDF)
'data.frame': 12 obs. of  5 variables:
 $ y       : int  1 2 3 4 5 6 10 11 12 1 ...
 $ person  : Factor w/ 4 levels "Alice","Bob",..: 1 1 1 2 2 2 3 3 3 4 ...
 $ female  : int  1 1 1 0 0 0 1 1 1 0 ...
 $ mood    : num  4.5 3.5 4 2.1 3 3.5 1.6 2.8 4.3 4.6 ...
 $ pmn.mood: num  4 4 4 2.87 2.87 ...

These summary variables can be included, along with their time-varying analog, in your multilevel model as a person-level predictor of the random intercept:

m2 <- lmer(y ~ female + mood + pmn.mood + (1 | person), data = myDF) 
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ female + mood + pmn.mood + (1 | person)
   Data: myDF

REML criterion at convergence: 32.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.36975 -0.36411  0.08105  0.39417  1.26723 

Random effects:
 Groups   Name        Variance Std.Dev.
 person   (Intercept) 12.8519  3.5850  
 Residual              0.9209  0.9596  
Number of obs: 12, groups:  person, 4

Fixed effects:
            Estimate Std. Error t value
(Intercept)  19.6480    10.6931   1.837
female        2.5451     3.6393   0.699
mood          0.5357     0.4125   1.299
pmn.mood     -5.0845     2.9531  -1.722

Correlation of Fixed Effects:
         (Intr) female mood  
female   -0.247              
mood      0.000  0.000       
pmn.mood -0.961  0.080 -0.140

If I calculated a person mean for female, the value for the person mean variable would be the same as their value on female and would not be permissible as a predictor in a model that also had female.