Zero Inflated Negative Binomial not producing standard errors - R

Question

I've recently learned about using ZI-Negative Binomial in R, but I haven't been able to figure out why my model results are not calculating a standard error.

fit.zinb<-zeroinfl(comp_counts~bc,dist="negbin",link="logit",data=master_region)
summary(fit.zinb)

Output:

Call:
zeroinfl(formula = comp_counts ~ bc, data = master_region, dist = "negbin", link = "logit")

Pearson residuals:
     Min       1Q   Median       3Q      Max 
 -0.2402  -0.2385  -0.2385  -0.2385 461.0861 

Count model coefficients (negbin with log link):
              Estimate Std. Error z value Pr(>|z|)
(Intercept)  3.534e-01         NA      NA       NA
bc           2.301e-09         NA      NA       NA
Log(theta)  -2.826e+00         NA      NA       NA

Zero-inflation model coefficients (binomial with logit link):
             Estimate Std. Error z value Pr(>|z|)
(Intercept) -12.95553         NA      NA       NA
bc           -0.00567         NA      NA       NA

Theta = 0.0593 
Number of iterations in BFGS optimization: 48 
Log-likelihood: -3.431e+05 on 5 Df

The dependent variable comp_counts is indeed inflated with 0's (83% of observations for the dependent variable are actually 0's).

To provide additional details on the variables, comp_counts is a count, ranging from 0 to 2754, and bc is a continuous variable. However, the independent variable bc also has several 0's (99.4% of observations are 0). Could this be causing the problem in the model?

Thanks for your help! Apologies for the novice question. Happy to provide the data if it is helpful.

Edit: the same issue (no standard errors) also occurs when I use the zero inflated poisson model. No warnings are produced, so not sure what's happening

Edit 2: Performing the regression models by reversing the dependent and independent variable does produce results with standard errors (ie, the dependent is now bc and independent comp_counts). But this is not the direction of causality that I seek to explore

Call:
zeroinfl(formula = bc ~ comp_counts, data = master_region, dist = "poisson", link = "logit")

Pearson residuals:
       Min         1Q     Median         3Q        Max 
  -0.09871   -0.07143   -0.07143   -0.07143 4732.55835 

Count model coefficients (poisson with log link):
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) 1.321e+01  3.179e-05  415741   <2e-16 ***
comp_counts 7.672e-03  3.626e-06    2116   <2e-16 ***

Zero-inflation model coefficients (binomial with logit link):
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  5.2781149  0.0232286 227.224   <2e-16 ***
comp_counts -0.0002349  0.0010703  -0.219    0.826    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

EDIT 3: The data is attached here

Answer 1

Always a good idea to plot the data first (especially if like me you're not familiar with it): I like ggplot for its flexibility, but it doesn't really matter what tool you use. I arrived at this after a few iterations:

library(ggplot2); theme_set(theme_bw())
(ggplot(master_region,aes(bc,comp_counts)) 
    ## show overlapping/replicated points as circles (semi-transparent)
    + stat_sum(alpha=0.5) 
    ## cosmetic: make minimum size (non-replicated) points larger             
    + scale_size(range=c(2,8))
    ## put both x- and y-axes on a log(1+x) scale
    + scale_x_continuous(trans="log1p")
    + scale_y_continuous(trans="log1p")
    ## fit and display a flexible smooth curve for point data, **not** 
    ## including points at x=0
    + geom_smooth(data=subset(master_region,bc>0),
                  method="gam",method.args=list(family="quasipoisson"))
)

Additionally, with(master_region,table(bc=bc>0, cc=comp_counts>0)) gives:

       cc
bc       FALSE   TRUE
  FALSE 302109  64678
  TRUE    1420    452

So far we know:

• most of your data fall at bc==0, comp_counts==0
• the data with bc==0, comp_counts>0 have an extremely wide range in comp_counts (0-2754)
• the data have an extremely wide range in bc (0-1.84e8), although the range of comp_counts is narrower for bc>0 (only up to 108)

This is going to be extremely challenging to fit from a computational standpoint. The problems with your standard deviations stem ultimately from the fact that R is having trouble computing the Hessian, i.e. the matrix that describes the curvature of the log-likelihood surface at the estimated parameter values.

One very simple solution is to scale your predictor variable: this sets the mean of the predictor to 0 and the SD to 1, which generally improves performance.

fit.zinb0 <- zeroinfl(comp_counts~scale(bc),dist="negbin",link="logit",
                     data=master_region)

You'll have to back-scale your parameters if you want to know the values on the original scale (i.e. multiply the bc coefficient by sd(bc) ... you can also back-transform the intercept, I'm too lazy to post this at the moment). But everything else about the model (Z-scores/p-values of predictors, overall goodness of fit, predictions, etc.) will be unchanged.

Alternately, you could consider fitting on the log(1+bc) scale:

fit.zinb <- zeroinfl(comp_counts~log10(1+bc),dist="negbin",link="logit",
                     data=master_region)

Finally, independently of these I would consider whether you want to model the behaviour at bc==0 separately from the behaviour for bc>0.