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

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),
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")),
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).

• Why do you have both mean stress and centered stress entered in your model? Jan 14 at 16:51
• @ShawnHemelstrand Because I want to disentangle the within- and between-person associations. Jan 14 at 19:33
• Gotcha. Im just going to assume that based off the fact you have tried repeatedly to get the model to converge that this is due to the model not matching the data. In other words, I would attempt to reduce the complexity of the model and see if it has problems, as this is a common issue with GLMMs. Jan 15 at 1:42
• I don't get this warning when removing the autocorrelation. But what can I do if I still want to account for autocorrelation? Is there another approach? Initially, I used lme and it worked well, but the residuals were not normally distributed, then I can't use lme. Jan 15 at 8:45
• That would make sense given the autocorrelation in your summary is almost zero. What happens when you remove random slopes and keep in the autocorrelation? Jan 15 at 8:46

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.