Mixed effects model with zero inflated data - Error message using zeroinfl() from 'pscl' package

Question

I am having issues with the function zeroinfl() from the 'pscl' package. Here is an overview of my situations:

I am trying to find out if the non-native stem density in a plot is influenced by the focal species at that plot. I am using a mixed effects model with the random effect being the site (I collected data at 6 different sites). This is what the mixed model looks like using the glmer() function from the 'lme4' package:

nonstem.model <- glmer(nonstem ~ focalspecies + (1|site), family = "poisson", data = data, na.action=na.omit)

The problem is that my data is zero inflated, meaning that there were a lot of plots with no non-native species present. So I tried using the zeroinfl() function from the 'pscl' package.

nonstem.ZIP = zeroinfl(nonstem ~ focalspecies + (1|site), dist="poisson", link = "logit", data = data)

But I received an error message:

Error in contrasts<-(*tmp*, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels In addition: Warning message: In Ops.factor(1, site) : ‘|’ not meaningful for factors

So I gathered that I cannot have a random effect here, and so I changed it to a fixed effect.

nonstem.ZIP = zeroinfl(nonstem ~ focalspecies + site, dist="poisson", link = "logit", data = data)

However, now I am getting this error message:

Error in solve.default(as.matrix(fit$hessian)) : system is computationally singular: reciprocal condition number = 4.70937e-36

If I change the distribution from "poisson" to "negbin" then a similar error message comes up:

Error in solve.default(as.matrix(fit$hessian)) : system is computationally singular: reciprocal condition number = 2.92265e-19

Does anyone know what this error message means? Or if there is a different package/function that I can use? Any help is much appreciated.

Answer 1

Concerning your comment, you can use base R's stats::anova to compare two models.

Here is a reproducible example using the Salamanders sample data from glmmTMB

library(glmmTMB);
fit1 = glmmTMB(
    count ~ spp * mined + (1|site),
    zi = ~ spp * mined,
    data = Salamanders,
    family = poisson);
fit2 = glmmTMB(
    count ~ spp + mined + (1|site),
    zi = ~ spp + mined,
    data = Salamanders,
    family = poisson);

anova(fit1, fit2)
#Data: Salamanders
#Models:
#fit2: count ~ spp + mined + (1 | site), zi=~spp + mined, disp=~1
#fit1: count ~ spp * mined + (1 | site), zi=~spp * mined, disp=~1
#     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
#fit2 17 1785.5 1861.4 -875.75   1751.5
#fit1 29 1778.1 1907.7 -860.04   1720.1 31.405     12   0.001708 **
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Credit where credit is due: The example was adjusted from here.