CODEWITHSUNDEEP
rcorrelationlme4
Jan Teunis
Jan Teunis

Reputation: 77

lmer - Model failed to converge with 1 negative eigenvalue

I'm trying to get some results in R, but after hours of searching I can't get past this error.

full code

Montenegro_Model2C2 <- lmer(ThuisvoelenZ ~ 1 + Gepest_voelenRZ + Instructietaal_Andere + (1 + Gepest_voelenRZ|IDSCHOOL) + (1 + Instructietaal_Andere|IDSCHOOL) + (1|IDCLASS), data = Montenegro, REML = F) 
summary(Montenegro_Model2C2)

Gives me the following warning:

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
3: Model failed to converge with 1 negative eigenvalue: -1.0e-01

Followed by an output with wrong correlation values (1.00)

Random effects:
 Groups     Name                  Variance Std.Dev. Corr
 IDCLASS    (Intercept)           0.103856 0.32227      
 IDSCHOOL   (Intercept)           0.002433 0.04933      
            Instructietaal_Andere 0.025042 0.15825  1.00
 IDSCHOOL.1 (Intercept)           0.020872 0.14447      
            Gepest_voelenRZ       0.010772 0.10379  0.46
 Residual                         0.771272 0.87822      
Number of obs: 4814, groups:  IDCLASS, 359; IDSCHOOL, 140

Later in the script, at Model2C4, the correlation is even worse. What am I doing wrong here?

The full script can be found here. OLP Functions right here, and the dataset here.

Upvotes: 1

Views: 9957

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226162

I don't think you're doing anything wrong. There are a few things you can do to confirm:

  • use allFit() to try a variety of optimizers and make sure the results are all sufficiently similar in whatever aspects you're interested in. The development version of broom.mixed() has a tidy() method for these results (you can use remotes::install_github("bbolker/broom.mixed")).
aa <- allFit(Montenegro_Model2C2)

Compare fits:

glance(aa) |> select(optimizer, AIC, NLL_rel) |> arrange(NLL_rel)
# A tibble: 7 × 3
  optimizer                        AIC      NLL_rel
  <chr>                          <dbl>        <dbl>
1 nloptwrap.NLOPT_LN_NELDERMEAD 12849. 0           
2 nlminbwrap                    12849. 0.0000000443
3 nloptwrap.NLOPT_LN_BOBYQA     12849. 0.000000231 
4 nmkbw                         12849. 0.00000143  
5 bobyqa                        12850. 0.194       
6 Nelder_Mead                   12851. 0.812       
7 optimx.L-BFGS-B               12851. 0.812

NLL_rel gives the negative log-likelihood relative to the best fit (NLL_rel = 0). bobyqa, Nelder-Mead, and L-BFGS-B are pretty bad, all the rest are close (nloptwrap.NLOPT_LN_BOBYQA is the default optimizer, which is what was used in the original fit).

Comparing parameters:

tidy(aa, conf.int = TRUE) |> 
     arrange(effect, term, estimate) |> 
     select(-c(std.error, statistic))
# A tibble: 77 × 7
   optimizer                     effect group term   estimate conf.low conf.high
   <chr>                         <chr>  <chr> <chr>     <dbl>    <dbl>     <dbl>
 1 optimx.L-BFGS-B               fixed  NA    (Inte…    0.102   0.0470     0.157
 2 Nelder_Mead                   fixed  NA    (Inte…    0.102   0.0470     0.157
 3 nloptwrap.NLOPT_LN_NELDERMEAD fixed  NA    (Inte…    0.105   0.0497     0.161
 4 nlminbwrap                    fixed  NA    (Inte…    0.105   0.0497     0.161
 5 nloptwrap.NLOPT_LN_BOBYQA     fixed  NA    (Inte…    0.105   0.0497     0.161
 6 nmkbw                         fixed  NA    (Inte…    0.105   0.0497     0.161
 7 bobyqa                        fixed  NA    (Inte…    0.105   0.0494     0.162
 8 optimx.L-BFGS-B               fixed  NA    Gepes…   -0.214  -0.249     -0.178
 9 Nelder_Mead                   fixed  NA    Gepes…   -0.214  -0.249     -0.178
10 nloptwrap.NLOPT_LN_BOBYQA     fixed  NA    Gepes…   -0.208  -0.243     -0.173

We can see (at least looking at the first few rows) that the 'bad' optimizers give slightly different parameter values, the others all look close.

If you don't want to work any harder, that should be good enough.

In any case, you have a singular fit (correlation in the IDSCHOOL term is estimated as 1.00), so depending on your modeling approach you might want to simplify the model anyway.

We can also look at some diagnostics to see if there's anything else concerning about the model:

aa <- augment(Montenegro_Model2C2)
ggplot(aa,
       aes(Gepest_voelenRZ, .resid, colour = factor(Instructietaal_Andere))) + geom_point() + geom_smooth()

enter image description here

This is a little worrying: there is a cluster of values at the maximum level of Gepest_voelenRZ for Instructietaal Andere==0 (Google translate says this means "feeling bullied", "instruction language 'other'" in Dutch??) that are badly fitted. Try checking these data points to see if there's anything suspicious? (Poor model fits don't necessarily lead to numerical instability in model fitting, but they don't help ...)

Upvotes: 6

Related Questions