My data has overdispersion but the Hurdle model estimated theta is 0. What am I doing wrong?

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

• 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. Apr 17, 2019 at 22:07
• 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? Apr 18, 2019 at 5:25
• 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? Apr 19, 2019 at 13:17
• 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 Apr 19, 2019 at 14:20
• April why not edit that output into your original question? Apr 19, 2019 at 15:56