[R] variance specification using glm and quasi

[Prof Brian Ripley]

On Mon, 12 Jun 2006, Jeff Miller wrote:

> Cameron and Trivedi in their 1998 Regression Analysis of Count Data refer to
> NB1 and NB2
>
> NB1 is the negative binomial model with variance = mu + (alpha * mu^1)
> yielding (1+alpha)*mu
>
> NB2 sets the power to 2; hence, variance = mu + (alpha*mu^2)
>
> I think that NB2 can be requested via
>
> negbin2<-glm(hhm~sex+age,family=quasi(var="mu^2",link="log"))
>
> Is that right?

No.  That is variance = phi*mu^2, not mu + (alpha*mu^2).

> If so, how I can get NB1? The quasi family appears to be very
> limited in variance specification options.

[Not so in the R-devel version of R: you can supply any variance 
function.]

You can use your own family in any version of R.  Package MASS has for 
many years supplied negative.binomial, which has mu + mu^2/theta, a 
parametrization of `NB2'.  It even provides glm.nb to estimate theta.

Note that (and this is explicitly in your reference) (1+alpha)*mu = phi*mu 
so that NB1 can be fitted as a quasipoisson GLM, although the 
quasilikelihood used is the *not* the likelihood of the model (which is 
not a GLM).  You could easily fit this model by maximum likelihood by 
direct maximization: p.445 of MASS provides a suitable template.

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595