brglm2 (version 0.9.2)

brnb: 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,
  ...
)

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

Arguments

formula

an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. The details of model specification are given under ‘Details’.

data

an optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which glm is called.

subset

an optional vector specifying a subset of observations to be used in the fitting process.

weights

an optional vector of ‘prior weights’ to be used in the fitting process. Should be NULL or a numeric vector.

offset

this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. One or more offset terms can be included in the formula instead or as well, and if more than one is specified their sum is used. See model.offset.

link

The link function. Currently must be one of "log", "sqrt" or "identity".

start

starting values for the parameters in the linear predictor.

etastart

starting values for the linear predictor.

mustart

starting values for the vector of means.

control

a list of parameters for controlling the fitting process. See brglmControl() for details.

na.action

a function which indicates what should happen when the data contain NAs. The default is set by the na.action setting of options, and is na.fail if that is unset. The ‘factory-fresh’ default is na.omit. Another possible value is NULL, no action. Value na.exclude can be useful.

model

a logical value indicating whether model frame should be included as a component of the returned value.

x, y

For glm: logical values indicating whether the response vector and model matrix used in the fitting process should be returned as components of the returned value.

For glm.fit: x is a design matrix of dimension n * p, and y is a vector of observations of length n.

contrasts

an optional list. See the contrasts.arg of model.matrix.default.

intercept

logical. Should an intercept be included in the null model?

singular.ok

logical; if FALSE a singular fit is an error.

...

For glm: arguments to be used to form the default control argument if it is not supplied directly.

For weights: further arguments passed to or from other methods.

Author

Euloge Clovis Kenne Pagui [aut] , Ioannis Kosmidis [aut, cre]

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 MASS::glm.nb() uses.

References

Cordeiro G M, McCullagh P (1991). Bias correction in generalized linear models. Journal of the Royal Statistical Society. Series B (Methodological), 53, 629-643. tools:::Rd_expr_doi("10.1111/j.2517-6161.1991.tb01852.x").

Firth D (1993). Bias reduction of maximum likelihood estimates. Biometrika. 80, 27-38. tools:::Rd_expr_doi("10.2307/2336755").

Kenne Pagui E C, Salvan A, Sartori N (2017). Median bias reduction of maximum likelihood estimates. Biometrika, 104, 923–938. tools:::Rd_expr_doi("10.1093/biomet/asx046").

Kosmidis I, Kenne Pagui E C, Sartori N (2020). Mean and median bias reduction in generalized linear models. Statistics and Computing, 30, 43-59. tools:::Rd_expr_doi("10.1007/s11222-019-09860-6").

Kosmidis I, Firth D (2009). Bias reduction in exponential family nonlinear models. Biometrika, 96, 793-804. tools:::Rd_expr_doi("10.1093/biomet/asp055").

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")

# \donttest{
## 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")
# }