The R Book

Negative binomial errors

Instead of using quasi-Poisson errors (as above) we could use a negative binomial model. This is in the MASS library and involves the function glm.nb. The modelling proceeds in exactly the same way as with a typical GLM:

model.nb1<-glm.nb(Days~Eth*Sex*Age*Lrn)
summary(model.nb1,cor=F)

Call:
glm.nb(formula = Days ~ Eth * Sex * Age * Lrn, init.theta =
1.92836014510701, link = log)

(DispersionparameterforNegativeBinomial(1.9284)family taken to be 1)

    Null deviance: 272.29 on 145 degrees of freedom
Residual deviance: 167.45 on 118 degrees of freedom
AIC: 1097.3

            Theta: 1.928
        Std. Err.: 0.269
2 x log-likelihood: -1039.324

The output is slightly different than a conventional GLM: you see the estimated negative binomial parameter (here called theta, but known to us as k and equal to 1.928) and its approximate standard error (0.269) and 2 times the log-likelihood (contrast this with the residual deviance from our quasi-Poisson model, which was 1301.1; see above). Note that the residual deviance in the negative binomial model (167.45) is not 2 times the log-likelihood.

An advantage of the negative binomial model over the quasi-Poisson is that we can automate the model simplification with stepAIC:

model.nb2<-stepAIC(model.nb1)
summary(model.nb2,cor=F) Coefficients: (3 not defined because of singularities) Estimate Std. Error z value Pr(>|z|) (Intercept) 3.1693 0.3411 9.292 < 2e-16 *** EthN -0.3560 0.4210 -0.845 0.397848 SexM -0.6920 0.4138 -1.672 0.094459 . AgeF1 ...