1
$\begingroup$

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")

This is what the normal plot looks like

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
$\endgroup$
8
  • $\begingroup$ Why do you include observations where you know the fish should not occur? This seems to me it does nothing for your analysis. Would you, for example, want to run your model to predict how many fish are in an area where you know the correct answer is zero anyway? $\endgroup$
    – jbowman
    Commented Feb 1, 2023 at 23:37
  • $\begingroup$ It was recommended to me to look into zero-inflated models because there is an excess of zeros in the data. I already took out the sites in the geographic province where they do not exist. At the moment, I don't know which observations should be zero for the remaining sites, so I've been working with the zero-inflated negative binomial model. $\endgroup$
    – Nocomis
    Commented Feb 1, 2023 at 23:43
  • $\begingroup$ You actually don't explain what the data values represent. But if we assume the values are counts of fishes and the species you study is endangered, don't you expect to have many zeros? $\endgroup$
    – dipetkov
    Commented Feb 1, 2023 at 23:55
  • $\begingroup$ You are correct. Many zeros are expected. I'm trying to model those zeros to see if there has been any significant change over time. $\endgroup$
    – Nocomis
    Commented Feb 2, 2023 at 0:02
  • $\begingroup$ You provide few details about the data and the models. It might help to include more information. For example, how was the data collected? How you model time (to study changes over time)? Do the geographic sites have the same size? ... $\endgroup$
    – dipetkov
    Commented Feb 2, 2023 at 0:33

1 Answer 1

3
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Awesome! That makes a lot of sense and makes my analysis a lot easier. Thank you for all the help! $\endgroup$
    – Nocomis
    Commented Feb 3, 2023 at 17:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.