Model convergence problem; non-positive-definite Hessian matrix; NA values for likelihood/AIC/etc.

This warning (Model convergence problem; non-positive-definite Hessian matrix) states that at glmmTMB's maximum-likelihood estimate, the curvature of the negative log-likelihood surface is inconsistent with glmmTMB really having found the best fit (minimum): instead, the surface is downward-curving, or flat, in some direction(s).

It will usually be accompanied by NA values for the standard errors, log-likelihood, AIC, and BIC, and deviance. When you run summary() on the resulting model, you'll get the warning In sqrt(diag(vcov)) : NaNs produced.

These problems are most likely:

• when a model is overparameterized (i.e. the data does not contain enough information to estimate the parameters reliably)
• when a random-effect variance is estimated to be zero, or random-effect terms are estimated to be perfectly correlated ("singular fit": often caused by having too few levels of the random-effect grouping variable)
• when zero-inflation is estimated to be near zero (a strongly negative zero-inflation parameter)
• when dispersion is estimated to be near zero
• when complete separation occurs in a binomial model: some categories in the model contain proportions that are either all 0 or all 1

How do we diagnose the problem?

zinbm0 = glmmTMB(count~spp + (1|site), zi=~spp, Salamanders, family=nbinom2)

fixef(zinbm0)

ff <- fixef(zinbm0)$zi
round(plogis(c(sppGP=unname(ff[1]),ff[-1]+ff[1])),3)

Salamanders <- transform(Salamanders, GP=as.numeric(spp=="GP"))
zinbm0_A = update(zinbm0, ziformula=~GP)

fixef(zinbm0_A)[["zi"]]

zinbm0_B = update(zinbm0, ziformula=~(1|spp))
fixef(zinbm0_B)[["zi"]]
VarCorr(zinbm0_B)

zinbm1 = glmmTMB(count~spp + (1|site), zi=~mined, Salamanders, family=nbinom2)
fixef(zinbm1)[["zi"]]

## at present we need to specify the parameter by number; for
##  extreme cases need to specify the parameter range
## (not sure why the upper bound needs to be so high ... ?)
cc = confint(zinbm1,method="uniroot",parm=9, parm.range=c(-20,20))
print(cc)

## data taken from gamlss.data:plasma, originally
## http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/plasma.html
gt_load("vignette_data/plasma.rda")
m4.1 <- glm(calories ~ fat*fiber, family = Gamma(link = "log"), data = plasma)
m4.2 <- glmmTMB(calories ~ fat*fiber, family = Gamma(link = "log"), data = plasma)
ps  <- transform(plasma,fat=scale(fat,center=FALSE),fiber=scale(fiber,center=FALSE))
m4.3 <- update(m4.2, data=ps)
## scaling factor for back-transforming standard deviations
ss <- c(1,
        fatsc <- 1/attr(ps$fat,"scaled:scale"),
        fibsc <- 1/attr(ps$fiber,"scaled:scale"),
        fatsc*fibsc)
## combine SEs, suppressing the warning from the unscaled model
s_vals <- cbind(glm=sqrt(diag(vcov(m4.1))),
                glmmTMB_unsc=suppressWarnings(sqrt(diag(vcov(m4.2)$cond))),
                glmmTMB_sc=sqrt(diag(vcov(m4.3)$cond))*ss)
print(s_vals,digits=3)

Example 3.

Here is another example (from Samantha Sherman):

L <- gt_load("vignette_data/sherman.rda")

The first model gives the warning: "non-integer counts in an nbinom1 model" (indicating that we probably should use a different response distribution, or round the values if that seems appropriate).

summary(mod1)

We can immediately see that the dispersion is very small and that the zero-inflation parameter is strongly negative.

Running diagnostics on the model, these are the only problems reported.

d1 <- diagnose(mod1)

Let's try dropping the zero-inflation term:

mod2 <- update(mod1, ziformula=~0)

We also get a "false convergence (8)" warning; see below.

FIXME: this anticipates/duplicates some of the discussion near the end.

The summary() and diagnose() functions reveal only the large, negative dispersion parameter:

summary(mod2)

Diagnose:

diagnose(mod2)

The "false convergence" warning comes from the nlminb() optimizer, and is difficult to interpret and resolve. The documentation says that the cause of this warning is that:

the gradient ∇f(x) may be computed incorrectly, the other stopping tolerances may be too tight, or either f or ∇f may be discontinuous near the current iterate x.

The only practical options we have for satisfying ourselves that a false convergence warning is really a false positive are the standard brute-force solutions of (1) making sure the gradients are small and the Hessian is positive definite (these are already checked internally); (2) trying different starting conditions, including re-starting at the current optimum; and (3) trying a different optimizer. We'll try option 3 and refit with the BFGS option from optim():

mod2_optim <- update(mod2,
                control=glmmTMBControl(optimizer=optim,
                                       optArgs=list(method="BFGS")))

## Warning in glmmTMB(formula = taeCPUE ~ Avg.Relief + (1 | site.name/reef.name),
## : non-integer counts in a nbinom1 model

BFGS doesn't give us any warnings. Comparing the parameter estimates:

(parcomp <- cbind(nlminb=unlist(fixef(mod2)),
                  optim=unlist(fixef(mod2_optim))))

##                       nlminb       optim
## cond.(Intercept)  -4.6422905  -4.6152081
## cond.Avg.Relief    0.1422406   0.1674755
## disp.(Intercept) -39.2185102 -11.4800820

The conditional-model parameters are practically identical. The dispersion parameters look quite different (-39.2 vs. -11.5), but if we back-transform from the log scale (via exp()) we can see that they are both extremely small (\(9.28\times 10^{-18}\) vs. \(1.03\times 10^{-5}\)).

Simplify the model by switching from NB1 to Poisson:

mod3 <- update(mod2, family=poisson)

summary(mod3)

diagnose(mod3)

You can also check directly whether the Hessian (curvature) of the model is OK by examining the pdHess ("positive-definite Hessian") component of the sdr ("standard deviation report") component of the model:

mod3$sdr$pdHess                       

In general models with non-positive definite Hessian matrices should be excluded from further consideration.

Model convergence problem: eigenvalue problems

m1 <- glmmTMB(count~spp + mined + (1|site), zi=~spp + mined, Salamanders, family=genpois)

diagnose(m1)

In this example, the fixed-effect covariance matrix is NaN. It may have to do with the generalized Poisson (genpois) distribution, which is known to have convergence problems; luckily, the negative binomial (nbinom1 and nbinom2) and/or Conway-Maxwell Poisson (compois) are good alternatives.

Models with convergence problems should be excluded from further consideration, in general.

In some cases, extreme eigenvalues may be caused by having predictor variables that are on very different scales: try rescaling, and centering, continuous predictors in the model.

NA/NaN function evaluation

Warning in nlminb(start = par, objective = fn, gradient = gr) : NA/NaN function evaluation

This warning occurs when the optimizer visits a region of parameter space that is invalid. It is not a problem as long as the optimizer has left that region of parameter space upon convergence, which is indicated by an absence of the model convergence warnings described above.

The following warnings indicate possibly-transient numerical problems with the fit, and can be treated in the same way (i.e. ignored if there are no errors or convergence warnings about the final fitted model).

Cholmod warning 'matrix not positive definite'

In older versions of R (< 3.6.0):

Warning in f(par, order = order, ...) : value out of range in 'lgamma'

false convergence

This warning:

false convergence: the gradient ∇f(x) may be computed incorrectly, the other stopping tolerances may be too tight, or either f or ∇f may be discontinuous near the current iterate x

comes from the nlminb optimizer used by default in glmmTMB. It's usually hard to diagnose the source of this warning (this Stack Overflow answer explains a bit more about what it means). Reasonable methods for making sure your model is OK are:

• restart the model at the estimated fitted values
• try using a different optimizer, e.g. control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS"))

and see if the results are sufficiently similar to the original fit.

NA/NaN gradient evaluation

dat1 = expand.grid(y=-1:1, rep=1:10)
m1 = glmmTMB(y~1, dat1, family=nbinom2)

FIXME: this is no longer a "gradient evaluation" error ...

The error occurs here because the negative binomial distribution is inappropriate for data with negative values.

If you see this error, check that the response variable meets the assumptions of the specified distribution.

gradient length

Error in nlminb(start = par, objective = fn, gradient = gr) : gradient function must return a numeric vector of length x

Error in optimHess(par.fixed, obj$fn, obj$gr): gradient in optim evaluated to length x

Try rescaling predictor variables. Try a simpler model and build up. (If you have a simple reproducible example of these errors, please post them to the issues list.)