`check_singularity` doesn't work for `glmmTMB` about performance HOT 9 CLOSED

Can you give a reproducible example?

from performance.

Probably there are different tolerance levels? @bbolker, could you point to the code in glmmTMB where singularity is checked? This is the code used in performance:

function (x, tolerance = 1e-05, ...) 
{
    insight::check_if_installed("lme4")
    vc <- .collapse_cond(lme4::VarCorr(x))
    any(sapply(vc, function(.x) any(abs(diag(.x)) < tolerance)))
}

from performance.

@raffaem You could try check_singularity(model, tolerance = 1e-4) or similar (like 1e-3), and that would detect "singularity" earlier.

from performance.

"singular convergence (7)" is a warning message from nlminb, and is (sigh) completely different from a singular fit to a mixed model; see e.g. here in the R code base, and the documentation linked here for details.

• the glmmTMB package at present has no code for checking singularity (it would generally manifest as a non-positive-definite Hessian)
• If I'm reading it correctly, the check_singularity() function from performance is incomplete. Especially for covariance matrices larger than 2×2, it's very easy to get a singular fit where none of the variances (diagonal elements) is close to zero. lme4::isSingular() will do a much better job, but only works as written for the way that lme4 is parameterized - it would have to be adapted to work for glmmTMB (and that might not be easy).

from performance.

Thanks for clarification! The code for performance::check_singularity.merMod() should be equivalent to lme4::isSingular():

function(x, tolerance = 1e-05, ...) {
  insight::check_if_installed("lme4")
  theta <- lme4::getME(x, "theta")
  diag.element <- lme4::getME(x, "lower") == 0
  any(abs(theta[diag.element]) < tolerance)
}

from performance.

I think lme4::getME(x, "lower") is not available for glmmTMB.

from performance.

theta represents the parameters underlying the random effects covariance matrix in both packages, but the two packages use completely different parameterizations, so the vector can't be treated the same way. (For details, see the Bates et al JSS paper/lmer vignette from lme4, or the "covariance structures" vignette from glmmTMB ...)

isSingular tests the elements of theta corresponding to the diagonal elements of the scaled Cholesky factor of the RE covariance matrix. It's possible that computing chol() for each element of VarCorr() and testing the diagonal elements would work for glmmTMB, but I'm not sure. It could certainly be done by brute force by calculating eigen() for each element of VarCorr() and testing the magnitude of the smallest eigenvalue ...

from performance.

Something like (I'm sure this can be coded more elegantly):

Set up example and show (1) that fitting with lmer gives a singular fit and is correctly diagnosed; (2) glmmTMB is not correctly diagnosed

library(glmmTMB)
library(lme4)
set.seed(101)
dd <- expand.grid(x=factor(1:6), f = factor(1:20), rep = 1:5)
dd$y <- simulate_new(~ 1 + (x|f),
                     newdata = dd,
                     newparam = list(beta = 0,
                                     theta = rep(0,21),
                                     betad = 0))[[1]]

m1 <- lmer(y ~ 1 + (x|f), data = dd, REML = FALSE)
isSingular(m1)  ## TRUE
any(abs(diag(VarCorr(m1)[[1]])) < 1e-3) ## FALSE

m2 <- glmmTMB(y ~ 1 + (x|f), data = dd, REML = FALSE)
## 'singular convergence' warning from nlminb ** AND ** it's really singular
performance::check_singularity(m2)  ## FALSE

And a quick implementation of a singularity check (with the exception of the "component" structure, something like this would also work for any other mixed model framework with a VarCorr method ...

eigs <- list()
vv <- VarCorr(m2)
for (component in c("cond", "zi")) {
    for (i in seq_along(vv[[component]])) {
        eigs <- c(eigs, list(eigen(vv[[component]][[i]], only.values = TRUE)$values))
    }
}
any(sapply(eigs, function(x) min(x)) < 1e-6) ## TRUE

from performance.

I could also add details about the nlminb codes to the troubleshooting vignette, if they're not there already ...

from performance.