brnb {brglm2}  R Documentation 
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
formula 
an object of class 
data 
an optional data frame, list or environment (or object
coercible by 
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 
offset 
this can be used to specify an a priori known
component to be included in the linear predictor during fitting.
This should be 
link 
The link function. Currently must be one of 
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 
na.action 
a function which indicates what should happen
when the data contain 
model 
a logical value indicating whether model frame should be included as a component of the returned value. 
x , y 
For For 
contrasts 
an optional list. See the 
intercept 
logical. Should an intercept be included in the null model? 
singular.ok 
logical; if 
... 
For For 
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 biasreducing 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 biasreducing score adjustments in Firth (1993) and Kosmidis & Firth (2009). 
type = "AS_median"
: the median biasreducing 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.
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 variancecovariance 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 loglikelihood function.
Author(s)
Euloge Clovis Kenne Pagui [aut]
kenne@stat.unipd.it, Ioannis Kosmidis [aut, cre]
ioannis.kosmidis@warwick.ac.uk
References
Cordeiro G M, McCullagh P (1991). Bias correction in generalized linear models. Journal of the Royal Statistical Society. Series B (Methodological), 53, 629643. doi:10.1111/j.25176161.1991.tb01852.x.
Firth D (1993). Bias reduction of maximum likelihood estimates. Biometrika. 80, 2738. doi:10.2307/2336755.
Kenne Pagui E C, Salvan A, Sartori N (2017). Median bias reduction of maximum likelihood estimates. Biometrika, 104, 923–938. 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, 4359. doi:10.1007/s11222019098606.
Kosmidis I, Firth D (2009). Bias reduction in exponential family nonlinear models. Biometrika, 96, 793804. doi:10.1093/biomet/asp055.
Examples
## Example in Saha, K., & Paul, S. (2005). Biascorrected maximum
## likelihood estimator of the negative binomial dispersion
## parameter. Biometrics, 61, 179185.
#
# 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 biasreduced fit
fitBR_mean < update(fitML, type = "AS_mean")
# Median biasreduced fit
fitBR_median < update(fitML, type = "AS_median")
# Mixed biasreduced fit
fitBR_mixed < update(fitML, type = "AS_mixed")
# Mean biascorrected fit
fitBC_mean < update(fitML, type = "correction")
# Penalized likelihood with Jeffreysprior 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 reducedbias, and mixed reducedbias 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")