# hurdle model with negative binomial distribution of counts - error message and model selection

I'm working with over-dispersed count data, which is zero inflated (~2/3 zeros). I've fit a hurdle model using hurdle from {pscl}.

I have three related questions:

1) When fitting a hurdle model such as:

hurdle(count ~ factor1*factor2, data=data, dist="negbin")


I'm receiving a warning message that I'm having a hard time interpreting:

Warning message:
In sqrt(diag(vc_count)[kx + 1]) : NaNs produced


I also tried using dist="poisson", which does not return a warning.

Any idea what the warning message means?

2) When fitting both the Poisson and negative binomial GLMs to the counts, I wanted to analyze the model. I did a waldtest() and a lrtest() (Wald Test and Likelihood Ratio test respectively), comparing the fitted model to a null model with only an intercept. The comparison showed that the model with a negative binomial distribution for counts was not significant (potentially due to the warning?) and that the model with a Poisson distribution was significant. Is this an appropriate way to assess the model significance?

3) I wasn't sure how to assess the fits of the models (usually looking at plot(glm.object), so I subsetted the count>0 values and fit the GLM and looked at the residual-fitted values plot as well as the Q-Q plot, to assess the fit. Is that an acceptable method? Note that objects of which are class hurdle cannot be used in plot().

Any help or insights are much appreciated!

(1) It is hard to tell what exactly goes wrong here without a reproducible example. One possibility could be that there are no non-zero observations for certain combinations of factor levels. Another possibility could be that the theta estimate in the NB version of the model degenerates either towards zero or towards infinity and hence leads to numeric problems. It could also be something else, though...

(2) I wouldn't start testing the model that generated the warning before I figured out what went wrong in (1). I wouldn't just ignore the warning.

(3) I would recommend constructing the plots by hand rather than estimating a poorly fitting GLM. You can get fitted() and residuals() from the hurdle model and then call plot functions for scatter and QQ plots respectively.

In response to (1):

NaN is a reserved term in R. The fact that you're getting NaN means the function sqrt(diag(vc_count)[kx + 1]) is returning something that is not any element of an integer, logical, or raw vector. One possible explanation for NaN, per the documentation, is that "a complex number is regarded as NaN if either the real or imaginary part is NaN but not NA." The square root of a negative number contains an imaginary number and returns NaN in R:

sqrt(-1)
[1] NaN
Warning message:
In sqrt(-1) : NaNs produced


Given that the function that is causing you problems is sqrt() (inside the hurdle() environment), this seems like a likely explanation. I recommend checking to see if diag(vc_count)[kx + 1] returns a negative number.