negbin {aod}

Negative-Binomial Model for Counts

Description

The function fits a negative-binomial log linear model accounting for overdispersion in counts y.

Usage

  negbin(formula, random, data, phi.ini = NULL, warnings = FALSE, 
         na.action = na.omit, fixpar = list(),
         hessian = TRUE, control = list(maxit = 2000), ...)
  

Arguments

formula

A formula for the fixed effects. The left-hand side of the formula must be the counts y i.e., positive integers (y >= 0). The right-hand side can involve an offset term.

random

A right-hand formula for the overdispersion parameter(s) \phi.

data

A data frame containing the response (y) and explanatory variable(s).

phi.ini

Initial values for the overdispersion parameter(s) \phi. Default to 0.1.

warnings

Logical to control printing of warnings occurring during log-likelihood maximization. Default to FALSE (no printing).

na.action

A function name. Indicates which action should be taken in the case of missing value(s).

fixpar

A list with 2 components (scalars or vectors) of the same size, indicating which parameters are fixed (i.e., not optimized) in the global parameter vector (b, \phi) and the corresponding fixed values.
For example, fixpar = list(c(4, 5), c(0, 0)) means that 4th and 5th parameters of the model are set to 0.

hessian

A logical. When set to FALSE, the hessian and the variances-covariances matrices of the parameters are not computed.

control

A list to control the optimization parameters. See optim. By default, set the maximum number of iterations to 2000.

...

Further arguments passed to optim.

Details

For a given count y, the model is:

y~|~\lambda \sim Poisson(~\lambda)

with \lambda following a Gamma distribution Gamma(r,~\theta).
If G denote the gamma function, then:

P(\lambda) = r^{-\theta} * \lambda^{\theta - 1} * \frac{exp(-\frac{\lambda}{r})}{G(\theta)}

E[\lambda] = \theta * r

Var[\lambda] = \theta * r^2

The marginal negative-binomial distribution is:

P(y) = G(y + \theta) * \left(\frac{1}{1 + r}\right)^\theta * \frac{(\frac{r}{1 + r})^y}{y! * G(\theta)}

The function uses the parameterization \mu = \theta * r = exp(X b) = exp(\eta) and \phi = 1 / \theta, where X is a design-matrix, b is a vector of fixed effects, \eta = X b is the linear predictor and \phi the overdispersion parameter.
The marginal mean and variance are:

E[y] = \mu

Var[y] = \mu + \phi * \mu^2

The parameters b and \phi are estimated by maximizing the log-likelihood of the marginal model (using the function optim()). Several explanatory variables are allowed in b. Only one is allowed in \phi.
An offset can be specified in the formula argument to model rates y/T. The offset and the marginal mean are log(T) and \mu = exp(log(T) + \eta), respectively.

Value

An object of formal class “glimML”: see glimML-class for details.

Author(s)

Matthieu Lesnoff [e-mail], Renaud Lancelot [e-mail]

References

Lawless, J.F., 1987. Negative binomial and mixed Poisson regression. The Canadian Journal of Statistics, 15(3): 209-225.

See Also

glimML-class, glm and optim,
glm.nb in the recommended package MASS,
gnlr in package gnlm available at https://www.commanster.eu/rcode.html.

Examples

  # without offset
  data(salmonella)
  negbin(y ~ log(dose + 10) + dose, ~ 1, salmonella)
  library(MASS) # function glm.nb in MASS
  fm.nb <- glm.nb(y ~ log(dose + 10) + dose,
                  link = log, data = salmonella)
  coef(fm.nb)
  1 / fm.nb$theta # theta = 1 / phi
  c(logLik(fm.nb), AIC(fm.nb))
  # with offset
  data(dja)
  negbin(y ~ group + offset(log(trisk)), ~ group, dja)
  # phi fixed to zero in group TREAT
  negbin(y ~ group + offset(log(trisk)), ~ group, dja,
    fixpar = list(4, 0))
  # glim without overdispersion
  summary(glm(y ~ group + offset(log(trisk)),
    family = poisson, data = dja))
  # phi fixed to zero in both groups
  negbin(y ~ group + offset(log(trisk)), ~ group, dja,
    fixpar = list(c(3, 4), c(0, 0)))