Troubleshooting with glmmTMB

2020-07-01

This vignette covers common problems that occur while using glmmTMB. The contents will expand with experience.

If your problem is not covered below, there’s a chance it has been solved in the development version; try updating to the latest version of glmmTMB on GitHub.

Warnings

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:

How do we diagnose the problem?

Example 1.

Consider this example:

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

First, see if any of the estimated coefficients are extreme. If you’re using a non-identity link function (e.g. log, logit), then parameter values with \(|\beta|>10\) are suspect (for a logit link, this implies probabilities very close to 0 or 1; for a log link, this implies mean counts that are close to 0 or extremely large).

Inspecting the fixed-effect estimates for this model:

fixef(zinbm0)

The zero-inflation intercept parameter is tiny (\(\approx -17\)): since the parameters are estimated on the logit scale, we back-transform with plogis(-17) to see the at the zero-inflation probability for the baseline level is about \(4 \times 10^{-8}\))). Many of the other ZI parameters are very large, compensating for the intercept: the estimated zero-inflation probabilities for all species are

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

Since the baseline probability is already effectively zero, making the intercept parameter larger or smaller will have very little effect - the likelihood is flat, which leads to the non-positive-definite warning.

Now that we suspect the problem is in the zero-inflation component, we can try to come up with ways of simplifying the model: for example, we could use a model that compared the first species (“GP”) to the rest:

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

This fits without a warning, although the GP zero-inflation parameter is still extreme:

fixef(zinbm0_A)[["zi"]]

Another possibility would be to fit the variation among species in the zero-inflation parameter as a random effect, rather than a fixed effect: this is slightly more parsimonious. This again fits without an error, although both the average level of zero-inflation and the among-species variation are estimated as very small:

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

The original analysis considered variation in zero-inflation by site status (mined or not mined) rather than by species - this simpler model only tries to estimate two parameters (mined + difference between mined and no-mining) rather than 7 (one per species) for the zero-inflation model.

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

This again fits without a warning, but we see that the zero-inflation is effectively zero in the unmined (“minedno”) condition (plogis(0.38-17.5) is approximately \(4 \times 10^{-8}\)). We can estimate the confidence interval, but it takes some extra work: the default Wald standard errors and confidence intervals are useless in this case.

## 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)

The lower CI is not defined; the upper CI is -2.08, i.e. we can state that the zero-inflation probability is less than plogis(-2.08) = 0.11.

More broadly, general inspection of the data (e.g., plotting the response against potential covariates) should help to diagnose overly complex models.

Example 2.

In some cases, scaling predictor variables may help. For example, in this example from @phisanti, the results of glm and glmmTMB applied to a scaled version of the data set agree, while glmmTMB applied to the raw data set gives a non-positive-definite Hessian warning.

## data taken from gamlss.data:plasma, originally
## http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/plasma.html
load(system.file("vignette_data","plasma.rda", package="glmmTMB"))
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):

load(system.file("vignette_data","troubleshooting.rda",package="glmmTMB"))

The first model gives the specified warning when it runs, as well as the other symptoms such as NA values for the likelihood:

summary(mod1)

We can immediately see that the dispersion is very small and that the zero-inflation parameter is strongly negative. However, we’ll develop some fancier machinery that checks the variance-covariance matrix or Hessian of the model, finds eigenvalues that are negative or close to zero, and identifies which model components contribute to those eigenvalues:

diagnose_vcov <- function(model, tol=1e-5, digits=2, analyze_hessian=FALSE) {
    vv <- vcov(model, full=TRUE)
    nn <- rownames(vv)
    if (!all(is.finite(vv))) {
        if (missing(analyze_hessian)) warning("analyzing Hessian, not vcov")
        if (!analyze_hessian) stop("can't analyze vcov")
        analyze_hessian <- TRUE
    }
    if (analyze_hessian) {
        par.fixed <- model$obj$env$last.par.best
        r <- model$obj$env$random
        if (!is.null(r)) par.fixed <- par.fixed[-r]
        vv <- optimHess(par.fixed, fn=model$obj$fn, gr=model$obj$gr)
        ## note vv is now HESSIAN, not vcov
    }
    ee <- eigen(vv)
    if (all(ee$values>tol)) {message("var-cov matrix OK"); return(invisible(NULL))}
    ## find negative or small-positive eigenvalues (flat/wrong curvature)
    bad_evals <- which(ee$values<tol)
    ## order worst to best
    bad_evals <- bad_evals[order(-ee$values[bad_evals])]
    ret <- lapply(bad_evals,
                  function(i) {
                      ## extract loadings
                      v <- setNames(ee$vectors[,i], nn)
                      ## order in decreasing magnitude & round
                      list(val=ee$values[i],vec=round(v[order(-abs(v))],digits))
                  })
    return(ret)
}

Running the diagnostics on the model:

(d1 <- diagnose_vcov(mod1))

This model has a very bad eigenvalue that is mostly driven by the zero-inflation parameter, and a little bit by the dispersion parameter. Let’s try dropping the zero-inflation term:

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

We still get the warning, and the NA-valued likelihoods (and the very small dispersion parameter). Diagnose:

diagnose_vcov(mod2)

We can see that the dispersion parameter is still problematic. Simplify the model by switching from NB1 to Poisson:

mod3 <- update(mod2, family=poisson)
summary(mod3)

There are no warnings, the model looks OK now, and the diagnostic function agrees:

diagnose_vcov(mod3)

You can also check directly whether 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

(FIXME: add an accessor method for this?)

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)

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:

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

Errors

NA/NaN gradient evaluation

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

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

inspecting Hessian (experimental!)

Can we directly inspect the Hessian to find where the problems are?

##' @param fit a \code{glmmTMB} fit
##' @param h Hessian (if not provided, will be computed from the Jacobian of the gradient
##' @param eval.eps numeric tolerance for 'bad' eigenvalues
##' @param evec.eps numeric tolerance for 'bad' eigenvector elements
diagnose_hessian <- function(fit,h=NULL, eval.eps=1e-5,evec.eps=1e-2) {
    ## pull out the TMB object from the fit
    obj <- fit$obj
    ee <- environment(obj$fn)
    ## extract parameters
    pp <- ee$last.par[-ee$random]
    ## easiest way to get names corresponding to all of the parameters
    nn <- tryCatch(colnames(vcov(fit,full=TRUE)),
                   ## fall-back position
                   error = function(e) make.unique(names(pp)))
    ## fit$sdr$pdHess
    if ("sdr" %in% names(fit)) {
        cat("bad params according to sdreport:",
            paste(nn[!is.finite(suppressWarnings(sqrt(diag(fit$sdr$cov.fixed))))],
                  collapse=", "),"\n")
    }
    ## two ways to compute the Hessian
    ## (1) directly from the objective function, via finite difference+Richardson extrapolation
    ## h1 <- hessian(obj$fn, pp)
    ## (2) use the gradient and compute its Jacobian (faster and probably more stable)
    if (is.null(h)) {
        if (!require(numDeriv)) stop("need numDeriv package installed") 
        h <- jacobian(obj$gr, pp)
    }
    ## double-check we get the same answer (approximately)
    ## all.equal(h1,h,tolerance=1e-5)
    ## now investigate the Hessian
    eigs <- eigen(h)
    ## non-positive definite means some of the eigenvectors are <= 0
    bad <- which(eigs$values/max(eigs$values)<=eval.eps)
    if (length(bad)==0) {
        cat("Hessian seems OK\n")
        return(invisible(h))
    }
    cat(sprintf("max eigenvalue = %1.3g",eigs$values[1]),"\n")
    for (b in bad) {  ## there could be more than one 'bad' direction/eigenvector ..
        cat(sprintf("Hessian eigenvalue %d = %1.3g (relative val = %1.3g)",
                    b,eigs$values[b],eigs$values[b]/eigs$values[1]),"\n")
        bad_vec <- eigs$vectors[,b]
        bad_elements <- which(abs(bad_vec)>evec.eps)
        cat("   bad elements:",nn[bad_elements],"\n")
    }
    cat("SDs computed from sqrt(diag(solve(H))):",
        paste(suppressWarnings(sqrt(diag(solve(h)))), collapse=", "),"\n")
    return(invisible(h))
}
diagnose_hessian(zinbm0)