Evaluating dispersion parameter and goodness of fit for binomial regression in r

First, I have looked at other posts but none address this specific question. I am trying to evaluate goodness of fit for a negative binomial regression, and what the dispersion parameter indicates in this case. First I tried Poisson (and quasipoisson) glm, but the goodness of fit test resulted in a rejection of fit.

Here are my data:

season  e.total hunting gosh    o.traps
2012    93      7       0       0
2013    61      10      5       14
2014    4       20      58      9
2015    10      10      27      15
2016    3       5       64      1
2017    0       4       46      0


Briefly, e.total is the number of eggs of an endangered species that were depredated due to a predator. Hunting, gosh, and o.traps are the number of those predators removed via human control methods. I am trying to determine which predictor has strongest relationship to e.total.

First I tried poisson and quasiposisson with all variations of the 3 predictors, but whenever I tested for goodness of fit, I would always reject.

So I try negative binomial:

nb.1 <- glm.nb(e.total ~ hunting + gosh + o.traps, data = pred.data)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.77051    0.86806   4.344   0.0491 *
hunting      0.10873    0.12086   0.900   0.4633
gosh        -0.07331    0.01841  -3.983   0.0576 .
o.traps     -0.02903    0.03510  -0.827   0.4951


No great result. Also I should note I try everything up to saturated model and get similar results. So I tried collapsing the hunting and o.traps data into one variable and call it non.gosh:

nb.5 <- glm.nb(e.total ~ gosh + non.gosh, data = pred.data)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  4.517387   0.299408  15.088 0.000632 ***
gosh        -0.069117   0.013874  -4.982 0.015547 *
non.gosh    -0.002412   0.018570  -0.130 0.904869
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(80519.82) family taken to be 3.990911)


Finally to my explicit question: Can I test goodness of fit here using:

1 - pchisq(summary(nb.5)$deviance, summary(nb.5)$df.residual
)


and the dispersion parameter seems way too large here. How do I interpret 80519.82? Am I correctly going about model selection?

• How bad is the goodness of fit, really? You have a nice, nearly linear relationship between e.total and log(1+gosh), suggesting the simple Poisson model fit with glm(total ~ I(log(1 + gosh)), family="poisson", data=pred.data) might work quite well. Your models use so many variables--five parameters for six observations!--that they're unlikely to be useful. A scatterplot matrix of your data makes it abundantly clear that gosh has the strongest relationship with e.total, regardless: see with(data.pred, pairs(cbind(e.total, hunting, o.traps, log.gosh=log(1+gosh)))).
– whuber
Commented Dec 13, 2017 at 23:59
• So do you think the best way to report this is to present the scatter plot matrix (and R^2 values?) and forget about the more complicated models? I guess I was just overthinking this. Commented Dec 14, 2017 at 1:16

I agree with whuber that fitting more than one parameter when you have only 6 samples is not a great idea. If that was just an example and your data has more rows, you may want to consider the liklihood ratio test to compare models, and determine if nb should be used instead of poisson. This post may help: Compare poisson and negative binomial regression with LR test

• Nope, that is the full dataset. Since that is the case, do you agree that reporting the difference in linear associations (log transformed in the case of gosh) is sufficient to illustrate the point? Commented Dec 14, 2017 at 1:41