# Convergence warning when using mixed effect model

I carried out field observations where I counted the number of birds in a plot, I repeated the observation 4 times for different months. The objective of this is to see if land-use influences the number of birds observed.

My formula is bird.count ~ lu + (1 | plot)

bird.count a count data with a lot of Zero observations and mean of count = 2 with variance = 5.9, lu is a factor with six levels and plot is a factor with 36 levels the where I carried out the repeated observations. In total I've got 144 observations.

Based on all of this I carried out a negative binomial mixed-effect regression, but got warnings regarding failure to converge and the Std.Error of the coefficients are very high and all the same.

I then tried averaging the bird.count of each plot and then had a total of 36 observations. I did this so I wouldn't have to use a mixed effect model. I ran a glm assuming a Gaussian distribution but the Std.Error of the coefficients are very high and all the same.

The things I have tried,

fit.x10 <- glmer.nb(bird.count ~ 1 + (1 | plot),
data = birds.data)

Error in negative.binomial(theta = 2543.33244700468) :
unused argument (theta = 2543.33244700468)
In theta.ml(Y, mu, weights = object@resp$weights, limit = limit, : iteration limit reached  fit.x11 <- glmer.nb(bird.count ~ lu + (1 | plot), data = birds.data)  Error in negative.binomial(theta = 1947.61685250933) : unused argument (theta = 1947.61685250933) In addition: Warning messages: 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00180515 (tol = 0.001, component 1) 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio - Rescale variables? 3: In theta.ml(Y, mu, weights = object@resp$weights, limit = limit,  :
iteration limit reached

fit.y10 <- glmmTMB(bird.count ~ (1 | plot), data = birds.data,
family = nbinom2, ziformula = ~ 0, se = TRUE,
verbose = FALSE, doFit = TRUE)

Family: nbinom2  ( log )
Formula:          bird.count ~ (1 | plot)
Data: birds.data

AIC      BIC   logLik deviance df.resid
501.6    510.5   -247.8    495.6      141

Random effects:

Conditional model:
Groups Name        Variance Std.Dev.
plot   (Intercept) 1.434    1.198
Number of obs: 144, groups:  plot, 36

Overdispersion parameter for nbinom2 family (): 14.2

Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.1303     0.2333   0.558    0.577

fit.y11 <- glmmTMB(bird.count ~ lu + (1 | plot),
data = birds.data, family = nbinom2,
ziformula = ~ 0, se = TRUE, verbose = FALSE,
doFit = TRUE)

> fit.y11 <- glmmTMB(bird.count~ lu+  (1|plot), data= birds.data, family = nbinom2,
+                    ziformula = ~0, se = TRUE, verbose = FALSE, doFit = TRUE)
> summary(fit.y11)
Family: nbinom2  ( log )
Formula:          bird.count ~ lu + (1 | plot)
Data: birds.data

AIC      BIC   logLik deviance df.resid
454.1    477.9   -219.1    438.1      136

Random effects:

Conditional model:
Groups Name        Variance Std.Dev.
plot   (Intercept) 0.1946   0.4412
Number of obs: 144, groups:  plot, 36

Overdispersion parameter for nbinom2 family ():   14

Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept)           -20.51    5531.89  -0.004    0.997
lureserve              21.94    5531.89   0.004    0.997
lunational park        21.48    5531.89   0.004    0.997
lunational park.set    20.62    5531.89   0.004    0.997
luplantation.ns        20.11    5531.89   0.004    0.997
luplantation.cv        21.33    5531.89   0.004    0.997

fit.z11 <- glm(bird.count ~ lu, data = birds.data.avg,
family = gaussian())

> summary(fit.z11)

Call:
glm(formula = bird.count ~ lu, family = gaussian, data = birds.data.avg)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-3.2500  -0.5625   0.0000   0.7917   3.7500

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)         4.441e-16  5.886e-01   0.000  1.00000
lureserve           4.750e+00  8.324e-01   5.707 3.17e-06 ***
lunational park     2.792e+00  8.324e-01   3.354  0.00217 **
lunational park.set 1.250e+00  8.324e-01   1.502  0.14361
luplantation.ns     7.500e-01  8.324e-01   0.901  0.37474
luplantation.cv     2.458e+00  8.324e-01   2.953  0.00606 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 2.078472)

Null deviance: 149.500  on 35  degrees of freedom
Residual deviance:  62.354  on 30  degrees of freedom
AIC: 135.94

Number of Fisher Scoring iterations: 2


Any advice on how to diagnose the problem or am I doing something wrong?

EDIT Based on the feed back, yes there is a separation issue, I tried using the GLMMadaptive package but it did not help.

mm1 <- mixed_model(fixed = bird.count ~ lu , random = ~ 1 | plot, data = birds.data , family = negative.binomial() )

This is the error I got,

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 and/or try setting the control argument
'iter_EM = 0'. 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


The alternative solution found is to remove the level from 'lu' variable which has zero observations or removing 'plot' the random effect.

• If land use is a factor and you have no other covariates, the standard errors can only be expected to be the same for each level. – Frans Rodenburg Apr 26 '19 at 4:18
• could you please elaborate, cause that doesn't make sense to me – Leon D Apr 26 '19 at 5:49
• @FransRodenburg, I think that's only true for linear models with balanced designs. – Ben Bolker Apr 26 '19 at 10:45
• not enough time for a proper answer, but: (1) you do indeed have a complete separation issue (all the values in your first group are zero, I think); you can solve this with GLMMadaptive or blme, see e.g. complete separation examples here. (2) do you have an alternative negative.binomial() function defined? What do you get from find("negative.binomial") ? – Ben Bolker Apr 26 '19 at 10:50

You could also give a try to the GLMMadaptive package that can fit the zero-inflated negative binomial mixed model. For example, see here. In addition, if you have a separation problem you could invoke a penalized estimation using the penalized argument; for an example see here and here.