2
$\begingroup$

I am confused by the dispersion parameter from my model. My data fails the overdispersion test. It's mean is 28.7, the variance is 18655.27. N=2916 of which 32% are zeros. How can theta equal 0 in the model output? I thought zero means no overdispersion, use a poisson regression instead. What am I not seeing? Or am I using the wrong model? Thanks!

subset:
       structure(list(linm = c(0, 1, 0, 0, 2, 0, 4, 0, 12, 1, 0, 0, 
    0, 1, 2, 0, 0, 1, 42, 86, 4, 7, 2, 2, 18, 6, 18, 1, 3, 4, 878, 
    1, 68, 2, 70, 46, 13, 1, 5, 3), NTU_log = c(NA, NA, NA, NA, NA, 
    NA, 2.29253475714054, NA, 2.17475172148416, NA, NA, NA, NA, 3.75653810258775, 
    NA, NA, NA, NA, 2.87356463957978, 3.79997350161952, NA, 3.42100000895834, 
    NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 2.66025953726586, 3.13983261752775, 
    NA, NA, NA, NA, NA, NA), SAL_log = c(1.03641529246456, 0.459971279972899, 
    0.200380785351367, 0.270603376498819, 0.166450377468644, 0.227787936640467, 
    2.37138015942553, 0.141277083534084, 2.65589659780768, 0.206688663567527, 
    0.101274462601439, 0.265321561197923, 0.101274462601439, 2.22437811492075, 
    0.18609550896943, 0.164433066717692, 0.18010600593939, 1.14706213493597, 
    2.69126054093852, 2.11723007602852, 1.36181171813436, 1.54308326206508, 
    1.25149036708158, 1.49508467518896, 2.0967530929631, 2.01559416591294, 
    2.26755758214694, 0.654541643507039, 1.80658341904009, 1.82615526195206, 
    3.22972624556844, 2.03744983544891, 2.57538655454381, 1.96414186636755, 
    3.06039512696199, 2.57439664177754, 2.44533611665484, 0.589573512538608, 
    0.436014303243863, 0.127745166659952), TEMP = c(19.4, 20.2, 19.9, 
    21.4, 20.1, 22.5, 17.2, 20.4, 13.4, 23.1, 23.8, 20.9, 24.5, 16.9, 
    26.5, 21.8, 23.2, 20, 14.5, 11.8, 19.6, 11.8, 20, 19.5, 18.8, 
    20.1, 20.5, 20.3, 19.9, 19.7, 19.5, 20, 10.6, 17.7, 19.8, 19.3, 
    20.5, 21, 21, 21.4), x2_grp = structure(c(2L, 2L, 2L, 2L, 2L, 
    2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 
    1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 
    2L, 2L, 2L), .Label = c("0", "1"), class = "factor")), row.names = c(2L, 
    6L, 8L, 13L, 16L, 18L, 23L, 26L, 28L, 31L, 34L, 37L, 39L, 41L, 
    44L, 47L, 49L, 52L, 53L, 54L, 57L, 59L, 60L, 64L, 66L, 71L, 73L, 
    76L, 81L, 84L, 87L, 89L, 92L, 96L, 97L, 105L, 110L, 112L, 114L, 
    117L), class = "data.frame")

My Code:

AER::dispersiontest(disp_mod_linm.b,  trafo=1)
 summary(hurdle(linm~ x2_grp+TEMP+SAL_log+NTU_log|x2_grp+TEMP+SAL_log+NTU_log, 
                 dist = c("negbin"), data = dt_b))

Output:

    Overdispersion test

data:  disp_mod_linm.b
z = 3.9909, p-value = 3.291e-05
alternative hypothesis: true alpha is greater than 0
sample estimates:
   alpha 
649.0598

Call:
hurdle(formula = substitute(i ~ x2_grp + TEMP + SAL_log + NTU_log | x2_grp + TEMP + SAL_log + NTU_log, list(i = as.name(x))), 
    data = dt_b, dist = c("negbin"))

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-0.5184 -0.3950 -0.3097 -0.2131 42.1207 

Count model coefficients (truncated negbin with log link):
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -8.77634   38.94351  -0.225 0.821699    
x2_grp1       0.03114    0.14090   0.221 0.825066    
TEMP         -0.07034    0.01864  -3.773 0.000162 ***
SAL_log       0.50837    0.05625   9.038  < 2e-16 ***
NTU_log       0.46750    0.06879   6.796 1.07e-11 ***
Log(theta)  -13.24934   38.94100  -0.340 0.733675    
Zero hurdle model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.57340    0.41341  -1.387    0.165    
x2_grp1      0.78623    0.11831   6.645 3.02e-11 ***
TEMP        -0.12321    0.01557  -7.915 2.47e-15 ***
SAL_log      3.30413    0.28177  11.726  < 2e-16 ***
NTU_log      0.86916    0.08698   9.993  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta: count = 0
Number of iterations in BFGS optimization: 41 
Log-likelihood: -8149 on 11 Df
$\endgroup$
  • $\begingroup$ I'm not sure, but it looks like a pretty complicated model for only 40 observations. It looks like log(theta) = -13 +/- 38, so that doesn't really rule out overdispersion. $\endgroup$ – HStamper Apr 17 at 22:07
  • $\begingroup$ HStamper, N=2916. 40 is just a subset of the data. Mean is 28.7, the variance is 18655.27, so mean and variance are definitely not equal. What am I not seeing? $\endgroup$ – April Apr 18 at 5:25
  • $\begingroup$ I would definitely start with a simpler model and add terms incrementally. What happens if you just remove all the predictors from the logit/zeros part of the formula? $\endgroup$ – HStamper Apr 19 at 13:17
  • $\begingroup$ Good question HStamper! Still theta =0 though. Part of the model output (to meet character limits here): Zero hurdle model coefficients (binomial with logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.71730 0.03157 22.72 <2e-16 *** --- Signif. codes: 0 '' 0.001 '' 0.01 '' 0.05 '.' 0.1 ' ' 1 Theta: count = 0 Number of iterations in BFGS optimization: 41 Log-likelihood: -1.454e+04 on 7 Df $\endgroup$ – April Apr 19 at 14:20
  • $\begingroup$ April why not edit that output into your original question? $\endgroup$ – mdewey Apr 19 at 15:56
0
$\begingroup$

I already wrote some of this on Facebook, but I figured I would fill it out a bit more here. First, I would make sure that the NAs have been properly dealt with. I don't think that is affecting your hurdle model, and I'm not positive what the default method for pscl is when it encounters NA values, but it would be best to impute them. I would recommend the mice package and a nonparametric model that would be able to handle the zeros in the response variable as well as the relatively high number of NAs (maybe random forest?).

Secondly, I think you should find the best model using whatever criteria you planned to use, getting rid of uninformative variables with the truncated negbin model. Then, I would implement the same model with a truncated Poisson distribution, and compare those two models to determine which is better using the log likelihood ratio test. I wouldn't stress about what the estimate of Theta is, just determine which distribution fits the data better. Is there a biological reason you are using a hurdle model instead of a zero-inflated model?

$\endgroup$
  • $\begingroup$ bholly thanks so much for your concise response. The NA's are handled with na.omit, na.rm isn't an accepted argument in hurdles, unfortunately. The NA's are not at random so I do not want to imput them. I have resorted to using GAM's which I dont know why I feel like is still incorrect, but I do. The output is way less wonky with GAM's. $\endgroup$ – April Apr 19 at 14:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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