Bias reduction for negative binomial regression models

Description

brnb is a function that fits negative binomial regression models using implicit and explicit bias reduction methods.

Usage

brnb(
  formula,
  data,
  subset,
  weights = NULL,
  offset = NULL,
  link = "log",
  start = NULL,
  etastart = NULL,
  mustart = NULL,
  control = list(...),
  na.action,
  model = TRUE,
  x = FALSE,
  y = TRUE,
  contrasts = NULL,
  intercept = TRUE,
  singular.ok = TRUE,
  ...
)

Arguments

Details

A detailed description of the fitting procedure is given in the iteration vignette (see, vignette("iteration", "brglm2") and Kosmidis et al, 2020). The number of iterations when estimating parameters are controlled by the maxit argument of brglmControl.

The type of score adjustment to be used is specified through the type argument (see brglmControl for details).

The available options are:

• type = "AS_mixed": the mixed bias-reducing score adjustments in Kosmidis et al (2020) that result in mean bias reduction for the regression parameters and median bias reduction for the dispersion parameter, if any; default.

• type = "AS_mean": the mean bias-reducing score adjustments in Firth (1993) and Kosmidis & Firth (2009).

• type = "AS_median": the median bias-reducing score adjustments in Kenne Pagui et al. (2017)

• type = "MPL_Jeffreys": maximum penalized likelihood with powers of the Jeffreys prior as penalty.

• type = "ML": maximum likelihood.

• type = "correction": asymptotic bias correction, as in Cordeiro & McCullagh (1991).

The choice of the parameterization for the dispersion is controlled by the transformation argument (see brglmControl for details). The default is "identity". Using transformation = "inverse" uses the dispersion parameterization that glm.nb uses.

Value

A fitted model object of class brnb inheriting from negbin and brglmFit. The object is similar to the output of brglmFit but contains four additional components: theta for the maximum likelihood estimate of the dispersion parameter as in glm.nb, vcov.mean for the estimated variance-covariance matrix of the regression coefficients, vcov.dispersion for the estimated variance of the dispersion parameter in the chosen parameterization (using the expected information), and twologlik for twice the log-likelihood function.

Author(s)

Euloge Clovis Kenne Pagui [ctb] kenne@stat.unipd.it, Ioannis Kosmidis [aut, cre] ioannis.kosmidis@warwick.ac.uk

References

Cordeiro GM & McCullagh, P (1991). Bias correction in generalized linear models. *Journal of the Royal Statistical Society. Series B(Methodological)*, **53**, 629-643.

Firth D (1993). Bias reduction of maximum likelihood estimates, Biometrika. **80**, 27-38.

Kenne Pagui EC, Salvan A, and Sartori N (2017). Median bias reduction of maximum likelihood estimates. *Biometrika*, **104**, 923–938

Kosmidis I, Kenne Pagui EC, Sartori N (2020). Mean and median bias reduction in generalized linear models. *Statistics and Computing*, **30** 43–59.

Kosmidis I and Firth D (2009). Bias reduction in exponential family nonlinear models. *Biometrika*, **96**, 793-804.

Examples

# Example in Saha, K., & Paul, S. (2005). Bias-corrected maximum
# likelihood estimator of the negative binomial dispersion
# parameter.  Biometrics, 61, 179--185.
#
# Number of revertant colonies of salmonella data
salmonella <- data.frame(freq = c(15, 16, 16, 27, 33, 20,
                                  21, 18, 26, 41, 38, 27,
                                  29, 21, 33, 60, 41, 42),
                         dose = rep(c(0, 10, 33, 100, 333, 1000), 3),
                         observation = rep(1:3, each = 6))

# Maximum likelihood fit with glm.nb of MASS
salmonella_fm <- freq ~ dose + log(dose + 10)
fitML_glmnb <- MASS::glm.nb(salmonella_fm, data = salmonella)

# Maximum likelihood fit with brnb
fitML <- brnb(salmonella_fm, data = salmonella,
              link = "log", transformation = "inverse", type = "ML")

# Mean bias-reduced fit
fitBR_mean <- update(fitML, type = "AS_mean")

# Median bias-reduced fit
fitBR_median <- update(fitML, type = "AS_median")

# Mixed bias-reduced fit
fitBR_mixed <- update(fitML, type = "AS_mixed")

# Mean bias-corrected fit
fitBC_mean <- update(fitML, type = "correction")

# Penalized likelihood with Jeffreys-prior penalty
fit_Jeffreys <- update(fitML, type = "MPL_Jeffreys")

# The parameter estimates from glm.nb and brnb with type = "ML" are
# numerically the same
all.equal(c(coef(fitML_glmnb), fitML_glmnb$theta),
            coef(fitML, model = "full"), check.attributes = FALSE)

# Because of the invariance properties of the maximum likelihood,
# median reduced-bias, and mixed reduced-bias estimators the
# estimate of a monotone function of the dispersion should be
# (numerically) the same as the function of the estimate of the
# dispersion:

# ML
coef(fitML, model = "dispersion")
1 / coef(update(fitML, transformation = "identity"), model = "dispersion")

# Median BR
coef(fitBR_median, model = "dispersion")
1 / coef(update(fitBR_median, transformation = "identity"), model = "dispersion")

# Mixed BR
coef(fitBR_mixed, model = "dispersion")
1 / coef(update(fitBR_mixed, transformation = "identity"), model = "dispersion")

## The same is not true for mean BR
coef(fitBR_mean, model = "dispersion")
1 / coef(update(fitBR_mean, transformation = "identity"), model = "dispersion")


## An example  from Venables & Ripley (2002, p.169).
data("quine", package = "MASS")
quineML <- brnb(Days ~ Sex/(Age + Eth*Lrn), link = "sqrt", transformation="inverse",
                data = quine, type="ML")
quineBR_mean <- update(quineML, type = "AS_mean")
quineBR_median <- update(quineML, type = "AS_median")
quineBR_mixed <- update(quineML, type = "AS_mixed")
quine_Jeffreys <- update(quineML, type = "MPL_Jeffreys")

fits <- list(ML = quineML,
             AS_mean = quineBR_mean,
             AS_median = quineBR_median,
             AS_mixed = quineBR_mixed,
             MPL_Jeffreys = quine_Jeffreys)
sapply(fits, coef, model = "full")