Zero-inflated negative binomial model still underfitting zeros

Question

I use R studio. I am trying to model an endangered fish by geographic province and count. Data is collected at the same exact sites, every year. All fish species are counted. For the species I am looking at, I summed the total number at each site and plotted it. Some sites exist where the fish should not occur (true zeroes), and other sites may have lost the species over time (happen to be zero). I already took out the sites in the geographic province where I know they do not exist. At the moment, I don't know which observations should be zero for the remaining sites.

ggplot ( sentinel, aes ( COUNT )) + geom_histogram ( binwidth = 5, position = "dodge" ) + xlab ( "Abundance of Brook Trout" ) + theme_bw () + ylab ( "Count")

I modeled the Poisson and negative binomial regression and then tested for over/under dispersion, which exists in both models. Then I modeled zero-inflated Poisson and zero-inflated negative binomial regressions.

Based on AIC values, the zero-inflated negative binomial is the best option. I checked for over/under dispersion in the zero-inflated negbinom model, and the value is close to 1 (which is good). I ran the likelihood of ratio test to compare the zero-inflated model with the normal model and it favors the zero-inflated. However, when I use the check_zeroinflation function, the zero-inflated negative binomial is still "underfitting zeroes (probable zero-inflation)." I'm not sure what the next step is for my analysis. Not sure if I've provided enough information, I am still learning. If this isn't the right approach, I'm open to ideas. Thank you.

# Main dataframe (FIBISTRATA refers to geographic province)
sentinel <- brook_done [ brook_done $ Type == "Sentinel", ]
 
# Normal poisson regression
brook_poisson1 <- glm ( COUNT ~ FIBISTRATA + Round,
                   family = "poisson", data = sentinel)
summary ( brook_poisson1 )

# Check for overdispersion
dispersiontest ( brook_poisson1 )
# Overdispersion exists

# Normal negative binomial model
brook_nb1 <- glm.nb ( COUNT ~ FIBISTRATA + Round, data = sentinel )
summary ( brook_nb1 )  

# Check models for zero-inflation
check_zeroinflation (brook_poisson1)
check_zeroinflation (brook_nb1)
# Zero-inflation probable in both models 

# Zero-inflated Poisson model
brook_zip1 <- zeroinfl ( COUNT ~ FIBISTRATA + Round, data = sentinel )
summary ( brook_zip1 )

# Test for over/under dispersion
zip_disp <- resid ( brook_zip1, type = "pearson")
N <- nrow ( sentinel )
p <- length ( coef ( brook_zip1 ))
sum ( zip_disp^2 ) / ( N - p )
# Overdispersion still exists

# Zero-inflated negative binomial
brook_zinb1 <- zeroinfl ( COUNT ~ FIBISTRATA + Round, link = "logit",
                     dist = "negbin", data = sentinel )
summary ( brook_zinb1 )

# Check for over/under dispersion
negbin_disp <- resid ( brook_zinb1, type = "pearson" )
N <- nrow (sentinel)
p <- length ( coef ( brook_zinb1 )) + 1
sum ( negbin_disp^2 ) / ( N - p )
# Value is close to 1 (1.06)

# Compare zero-inflated models
lrtest (brook_zip1, brook_zinb1)
# Zero-inflated binomial is better

# Compare models
AIC (brook_zip1, brook_zinb1)
# Zero-inflated negative binomial is far better

# Check zero-inflation
check_zeroinflation (brook_zinb1)
# Probable zero-inflation

Answer 1

I looked at the code, and the check_zeroinflation function ignores the zero inflation part of the model - it only uses a negative binomial model (and not the one fitted using the zero inflation model) - when calculating the predicted number of zeros. The relevant lines of code inside check_zeroinflation are:

   mu <- stats::fitted(x)
      ...
   if (model_info$is_negbin) {
      ...   
        else {
            theta <- x$theta
        }
    }

      ...

    if (!is.null(theta)) {
       pred.zero <- round(sum(stats::dnbinom(x = 0, size = theta, mu = mu)))
    }
    else {
        pred.zero <- round(sum(stats::dpois(x = 0, lambda = mu)))
    }
    structure(class = "check_zi", list(predicted.zeros = pred.zero, 
        observed.zeros = obs.zero, ratio = pred.zero/obs.zero, 
        tolerance = tolerance))

where theta is pulled directly from the model object and mu is the mean of the fitted model (which naturally includes the zeros from the zero inflation factor!)

So... you should ignore it when fitting a zero-inflation model of any type.