Description Usage Arguments Details Value Author(s) References Examples
brnb
is a function that fits negative binomial regression
models using implicit and explicit bias reduction methods.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
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 
For For 
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 
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 glm.nb
uses.
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 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.
Euloge Clovis Kenne Pagui [ctb] kenne@stat.unipd.it, Ioannis Kosmidis [aut, cre] ioannis.kosmidis@warwick.ac.uk
Cordeiro GM & McCullagh, P (1991). Bias correction in generalized linear models. *Journal of the Royal Statistical Society. Series B(Methodological)*, **53**, 629643.
Firth D (1993). Bias reduction of maximum likelihood estimates, Biometrika. **80**, 2738.
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**, 793804.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77  # 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")

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.