glmmTMB - Model convergence problem; non-positive-definite Hessian matrix

Question

I am trying to conduct a model with the glmmTMB package - depression as an outcome and stress as a predictor, including age, gender, working hours, and observation number (=time) as covariates. I also want to disentangle stress into within-subject and between-subject effects, and include autocorrelation (using the number of observations) in the model. Because the data is skewed, I am using a GLMM model with a gamma distribution.

First I used this code:

dep_gamma_ar1 = glmmTMB(dep_sum ~ cmean_gender + cmean_age + 
                          cmean_workh + obs_cent 
                             + mean_stress + cent_stress
                             + (cent_stress + obs_cent |id)
                             + ar1(factor(obs) + 0 | id),
                             data=dat,family=Gamma(link = "log"))
summary(dep_gamma_ar1)

However, I got the following warning:

Warning messages:
1: In fitTMB(TMBStruc) :
  Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
2: In fitTMB(TMBStruc) :
  Model convergence problem; false convergence (8). See vignette('troubleshooting')

After searching for solutions, I found recommendations to add the following optimizer control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")) to the code:

dep_gamma_ar2 = glmmTMB(dep_sum ~ cmean_gender + cmean_age + 
                          cmean_workh + obs_cent 
                             + mean_stress + cent_stress
                             + (cent_stress + obs_cent |id)
                             + ar1(factor(obs) + 0 | id),
                           control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")),
                             data=dat,family=Gamma(link = "log"))
summary(dep_gamma_ar2)

I am still getting a non-positive-definite Hessian matrix. Here is the model summary:

     AIC      BIC   logLik deviance df.resid 
      NA       NA       NA       NA     3520 

Random effects:

Conditional model:
 Groups Name         Variance  Std.Dev. Corr                               
 id     (Intercept)  7.114e-02 0.26672                                     
        cent_stress     9.308e-04 0.03051  0.14                               
        obs_cent     4.998e-05 0.00707  0.26        0.23                   
 id.1   factor(obs)1 1.026e-01 0.32038  -0.01 (ar1) -0.01 (ar1) -0.01 (ar1)
Number of obs: 3536, groups:  id, 123

Dispersion estimate for Gamma family (sigma^2): 4.69e-06 

Conditional model:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.9568053  0.0319644  29.934  < 2e-16 ***
cmean_gender  0.0625847  0.0724055   0.864 0.387388    
cmean_age    -0.0009886  0.0018623  -0.531 0.595537    
cmean_workh     0.0050827  0.0046462   1.094 0.273970    
obs_cent     -0.0013033  0.0008871  -1.469 0.141791    
mean_stress  0.0292378  0.0100527   2.908 0.003632 ** 
cent_stress  0.0147936  0.0044367   3.334 0.000855 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 

Is the code correct? How can I solve this warning message? I read some past threads like here, and also read the 'troubleshooting' file but I still don't really understand the problem (I am pretty new to these kinds of models).

I'd kindly appreciate your help.

Answer 1

Looking at what you have so far, I see the following:

• You have modeled your syntax correctly.
• You have attempted to optimize estimation with some of the built in arguments in this package.
• You have correctly assumed there are still problems.

Looking at your regression summary, I see the following potential reasons why there are still warnings:

• The autocorrelations are so weak they are almost zero.
• The random effects are somewhat complex.
• The fixed effect point estimates are fairly low as well, which contribute to the RE problem.
• Probably of most importance: your random effect variance is super low.

All of these seem to point to the issue that your model is simply too complex and doesn't represent your data well. It may help to look at some of your scatter plots and some of your by-subject variances to see where the issue may lie. Seeing how in the comments you have mentioned the issue going away after minimizing the complexity, this seems to be the correct path, but you'll need to do some data exploration to answer that concretely.

A couple articles on model complexity have been written or co-written by Douglas Bates, the man behind the lme4 package used for mixed models. I cite them below because they illustrate how model complexity can often cause issues of convergence.

• Bates et. al., 2018: Parsimonious mixed models.
• Matuschek et al., 2017: Balancing Type I error and power in linear mixed models