Mixed effect zero inflated negative binomial model: "the leading minor of order 1 is not positive definite"

Question

I am having trouble fitting a mixed effect zero inflated negative binomial model to my data using the GLMMadaptive package:

mixed_model(abun ~ TREATMENT + Size + R , random = ~1 | PLOT.TR,
              data = abunIP, family = zi.negative.binomial(), 
              zi_fixed = ~grass.cut)

returns: Error in chol.default(X[[i]], ...) : the leading minor of order 1 is not positive definite I can't find anything about this particular error message for these models.

My data has 978 rows. The resonse abun is count data (flowers). TREATMENT and PLOT.TR are factors indicating plot type and a unique plot ID, respectively. Size, and R are continuous variables. grass.cut indicates whether weeding occurred before a survey, which could reduce the abun value, often to 0, which is why I am trying to use a zero-inflated model.

UPDATE

Following Dimitris' answer below, I found that the combination of solutions 1 and 2 worked for me, though I made the change list(betas = rep(0,7)) because TREATMENT is a factor with 5 levels.

By only adding the initial_values argument, I got the error message

Error in mixed_fit(y, X, Z, X_zi, Z_zi, id, offset, offset_zi, family,  : 
  A large coefficient value has been detected during the optimization.
 Please re-scale you covariates. Alternatively, this may due to a
 divergence of the optimization algorithm, indicating that an overly
 complex model is fitted to the data. For example, this could be
 caused when including random-effects terms (e.g., in the
 zero-inflated part) that you do not need. Otherwise, adjust the
 'max_coef_value' control argument.

However, specifying iter_EM = 0 in addition resulted in no error messages.

I have a follow-up question: Conceptually, what is iter_EM = 0 doing, and is it addressing the previous error message I get by specifying initial_values appropriately?

Answer 1

It is difficult to see what may be causing the problem without a reproducible example, but some things to check are:

1. Change initial values for the fixed effects coefficients using the initial_values argument of mixed_model(), e.g.,

   mixed_model(abun ~ TREATMENT + Size + R , random = ~ 1 | PLOT.TR,
               data = abunIP, family = zi.negative.binomial(), 
               zi_fixed = ~ grass.cut,
               initial_values = list(betas = rep(0, 3)))
   

   potentially you could also change the initial value for the variance of the random intercepts using initial_values = list(betas = rep(0, 3)).

2. You could try setting the number of EM iterations to zero using mixed_model(..., iter_EM = 0) to only use the quasi-Newton part of the optimization procedure.

3. If you have a (complete) separation issue, you could try invoking the penalized argument of mixed_model() that will place a Student's t penalty for the fixed effects.