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 ...