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

Value

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

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.

References

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

Examples

Run this code

# NOT RUN {
  # 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))) 
  
# }

Run the code above in your browser using DataCamp Workspace