[R] glm.poisson.disp versus glm.nb
HankeA at mar.dfo-mpo.gc.ca
Mon Feb 2 16:13:50 CET 2004
This is a question about overdispersion and the ML estimates of the
parameters returned by the glm.poisson.disp (L. Scrucca) and glm.nb
(Venables and Ripley) functions. Both appear to assume a negative binomial
distribution for the response variable.
Paul and Banerjee (1998) developed C(alpha) tests for "interaction and main
effects, in an unbalanced two-way layout of counts involving two fixed
factors, when data are Poisson distributed, and when data are extra
dispersed". In R I coded their C(alpha) test statistic (called TNBI) for
interaction for the case where the counts are extra-dispersed, as well as
their test for extra-dispersion (called T_a). Using the Quine data set
(Quine, 1975) the authors collapse the orginal 4x2x2x2 study design into a
2x4 table and obtained estimates of TNBI=10.36 and T_a=90.81.
Using the dispersion estimate from glm.poisson.disp and the estimates for mu
I got exactly the same results for TNBI and T_a. This made me happy. Now I
thought to try the ML estimates from glm.nb to see if the results would be
the same but I am having difficulty relating the dispersion phi from
glm.poisson.disp to theta estimated by glm.nb.
According to the R help for glm.poisson.disp " Var(y_i) = mu_i(1+mu_i*phi)
". The help for glm.nb lead me to a book by V&R (1994) which indicates that
Var(y)=mu+mu^2/theta. From this I gathered that phi=1/theta but the
estimates do not relate to each other in this way unless one is in error. In
a document by L.P. Ammann he says a "negative binomial model can be
specified with mean mu and dispersion phi by taking theta=mu/(phi-1)". I had
a problem implementing this because in my mind mu is a vector whereas phi
and theta are scalars.
Consequently, I would like to know how to get phi from theta so that I can
compare the glm.poisson.disp and glm.nb methods for estimating dispersion.
Department of Fisheries and Oceans
St. Andrews Biological Station
531 Brandy Cove Road
St. Andrews, NB
More information about the R-help