Using glmmTMB
I have (so far) found 4 convergence scenarios:
model$fit$convergence == 0
model$fit$message == "relative convergence (4)"
warning: "model convergence problem; singular convergence (7)"
model$fit$convergence: 1
model$fit$message : "singular convergence (7)"
warning: "Model convergence problem; non-positive-definite Hessian matrix"
model$fit$convergence: 0
model$fit$message: "relative convergence (4)"
warning 1: "Model convergence problem; non-positive-definite Hessian matrix."
warning 2: " Model convergence problem; false convergence (8)"
model$fit$convergence: 1
model$fit$message: "false convergence (8)"
My problem is that I need to check the status programmatically. So far I have not found how to extract the warnings from the fitted object, so I have been working with the output just model$fit$convergence
and model$fit$message
, and as you can see, these are identical for 1) and 3) above.
So my questions:
If anyone wants to try this for themselves, then obviously you need data and code which produces these scenarios, so here it is:
# 1. Normal convergence:
set.seed(1)
data <- data.frame(y = rnorm(100), x = rnorm(100), group = rep(1:10, each = 10))
model <- glmmTMB(y ~ x + (1 | group), data = data)
model$fit$convergence # == 0
model$fit$message # "relative convergence (4)"
# 2. Singular fit
set.seed(3)
data <- data.frame(y = rnorm(100), x = rnorm(100), group = rep(1:10, each = 10))
model <- glmmTMB(y ~ x + (1 + x | group), data = data)
model$fit$convergence # == 1
model$fit$message # "singular convergence (7)"
# 3. Possibly a singular fit:
set.seed(5)
data <- data.frame(y = rnorm(100), x = rnorm(100), group = rep(1:10, each = 10))
model <- glmmTMB(y ~ x + (1 + x^2 | group), data = data)
model$fit$convergence # == 0
model$fit$message # "relative convergence (4)"
# 4. False Convergence
set.seed(1)
data <- data.frame(y = rnorm(100), x = rnorm(100), group = rep(1:10, each = 10))
data$y1 = data$y < -2
model <- glmmTMB(y1 ~ x + (1 + x | group), family = binomial, data = data)
model$fit$convergence # == 1
model$fit$message # "false convergence (8)"
tl;dr model should be OK if model$fit$convergence == 0 && model$sdr$pdHess
is satisfied.
These warnings come from two different places.
model$fit$convergence
and model$fit$message
come from whatever nonlinear optimizer was used to minimize the negative log-likelihood (by default,nlminb
). To understand the meanings of these convergence codes and messages, see the nlminb
documentation (or Stack Overflow, in the notoriously poorly documented cases of "singular convergence" or "false convergence" ...). R's interface to nlminb
returns a convergence code of 0 if the fit is OK (i.e., nlminb return codes 3 to 6) or 1 otherwise.
Note that "singular convergence"1 is not the same as a "singular fit" in the mixed-model sense (where the random-effects covariance matrix is singular, i.e. variances estimated as 0 or non-positive-definite correlation matrices ...)
More generally, if you're using an alternate optimizer (other than nlminb
), the specific messages/convergence codes will differ, but 0 is universally the convergence code indicating success.
non-positive-definite Hessian matrices occur when the covariance of "top-level" parameters (fixed effect coefficients and coefficients determining the random effects covariance matrix, along with coefficients governing zero inflation and dispersion) is not positive definite. This may be the result of a LMM-type singular fit (in which case the random effects variance will go to zero, the RE log-standard deviation [the scale on which glmmTMB actually fits the model] goes to -∞, and the likelihood surface becomes flat) or of some other numerical instability/unidentifiability/etc. In general you can check these by examining model$sdr$pdHess
(generally TRUE if the model is OK, although see here for an edge case)
nlminb
does not list the possible convergence codes. Here is a subset of "bad" codes from the PORT documentation:
7 — singular convergence: x may have too many free components. See §5.
8 — 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.
9 — function evaluation limit: no convergence after IV(MXFCAL) = IV(17) evaluations of f(x).
10 — iteration limit: no convergence after IV(MXITER) = IV(18) iterations.
11 — STOPX returned .TRUE.: you supplied a system-dependent STOPX (see §12) routine and hit the BREAK key.
12–13 — impossible: these are input IV(1) values only. (12 means allocate storage within IV and V and start the algorithm; this is the default IV(1) value supplied by [D]IVSET. 13 means just allocate storage and return. See §4a for an example.)
14 — storage has been allocated (after a call with IV(1) = 13 — see, for example, §4a below)
1 From the PORT documentation:
V(37) is the singular-convergence tolerance. A return with IV(1) = 7 occurs if a more favorable stopping test is not satisfied and if the algorithm thinks
f (x) − min {f(y) :||D(y-x)|| ≤ V(LMAXS) } < V(SCTOL) .|f(x)|,
where
D
is given by (4.1). When this test is satisfied, it appears that x has too many degrees of freedom — and you should ponder whether f was properly formulated. Default = max {10^{-10}, MACHEP^{2/3} }.