brnb | R Documentation |
brnb()
is a function that fits negative binomial regression
models using implicit and explicit bias reduction methods.
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,
...
)
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 |
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.
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.
Euloge Clovis Kenne Pagui [aut]
kenne@stat.unipd.it, Ioannis Kosmidis [aut, cre]
ioannis.kosmidis@warwick.ac.uk
Cordeiro G M, McCullagh P (1991). Bias correction in generalized linear models. Journal of the Royal Statistical Society. Series B (Methodological), 53, 629-643. \Sexpr[results=rd]{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. \Sexpr[results=rd]{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. \Sexpr[results=rd]{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. \Sexpr[results=rd]{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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1093/biomet/asp055")}.
## 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")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.