% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/brnb.R
\name{brnb}
\alias{brnb}
\title{Bias reduction for negative binomial regression models}
\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{
\item{formula}{an object of class \code{"\link[stats]{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 \sQuote{Details}.}

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

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

\item{weights}{an optional vector of \sQuote{prior weights} to be used
    in the fitting process.  Should be \code{NULL} or a numeric vector.}

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

\item{link}{The link function. Currently must be one of \code{"log"},
\code{"sqrt"} or \code{"identity"}.}

\item{start}{starting values for the parameters in the linear predictor.}

\item{etastart}{starting values for the linear predictor.}

\item{mustart}{starting values for the vector of means.}

\item{control}{a list of parameters for controlling the fitting
process. See \code{\link[=brglmControl]{brglmControl()}} for details.}

\item{na.action}{a function which indicates what should happen
    when the data contain \code{NA}s.  The default is set by
    the \code{na.action} setting of \code{\link{options}}, and is
    \code{\link[stats]{na.fail}} if that is unset.  The \sQuote{factory-fresh}
    default is \code{\link[stats]{na.omit}}.  Another possible value is
    \code{NULL}, no action.  Value \code{\link[stats]{na.exclude}} can be useful.}

\item{model}{a logical value indicating whether \emph{model frame}
    should be included as a component of the returned value.}

\item{x, y}{For \code{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 \code{glm.fit}: \code{x} is a design matrix of dimension
    \code{n * p}, and \code{y} is a vector of observations of length
    \code{n}.
  }

\item{contrasts}{an optional list. See the \code{contrasts.arg}
    of \code{model.matrix.default}.}

\item{intercept}{logical. Should an intercept be included in the
    \emph{null} model?}

\item{singular.ok}{logical; if \code{FALSE} a singular fit is an
    error.}

\item{...}{
    For \code{glm}: arguments to be used to form the default
    \code{control} argument if it is not supplied directly.

    For \code{weights}: further arguments passed to or from other methods.
  }
}
\value{
A fitted model object of class \code{\link[=brnb]{"brnb"}} inheriting
from \code{\link[=negbin]{"negbin"}} and \code{\link[=brglmFit]{"brglmFit"}}. The
object is similar to the output of \code{\link[=brglmFit]{brglmFit()}} but contains
four additional components: \code{theta} for the maximum likelihood
estimate of the dispersion parameter as in \code{\link[MASS:glm.nb]{MASS::glm.nb()}},
\code{vcov.mean} for the estimated variance-covariance matrix of the
regression coefficients, \code{vcov.dispersion} for the estimated
variance of the dispersion parameter in the chosen
parameterization (using the expected information), and
\code{twologlik} for twice the log-likelihood function.
}
\description{
\code{\link[=brnb]{brnb()}} is a function that fits negative binomial regression
models using implicit and explicit bias reduction methods.
}
\details{
A detailed description of the fitting procedure is given in the
iteration vignette (see, \code{vignette("iteration", "brglm2")} and
Kosmidis et al, 2020). The number of iterations when estimating
parameters are controlled by the \code{maxit} argument of
\code{\link[=brglmControl]{brglmControl()}}.

The type of score adjustment to be used is specified through the
\code{type} argument (see \code{\link[=brglmControl]{brglmControl()}} for details).

The available options are:
\itemize{
\item \code{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.
\item \code{type = "AS_mean"}: the mean bias-reducing score
adjustments in Firth (1993) and Kosmidis & Firth (2009).
\item \code{type = "AS_median"}: the median bias-reducing score
adjustments in Kenne Pagui et al. (2017)
\item \code{type = "MPL_Jeffreys"}: maximum penalized likelihood
with powers of the Jeffreys prior as penalty.
\item \code{type = "ML"}: maximum likelihood.
\item \code{type = "correction"}: asymptotic bias correction, as in
Cordeiro & McCullagh (1991).
}

The choice of the parameterization for the dispersion is controlled
by the \code{transformation} argument (see \code{\link[=brglmControl]{brglmControl()}} for
details).  The default is \code{"identity"}. Using \code{transformation = "inverse"} uses the dispersion parameterization that
\code{\link[MASS:glm.nb]{MASS::glm.nb()}} uses.
}
\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")
}

}
\references{
Cordeiro G M, McCullagh P (1991). Bias correction in generalized
linear models. \emph{Journal of the Royal Statistical Society. Series B
(Methodological)}, \strong{53}, 629-643. \doi{10.1111/j.2517-6161.1991.tb01852.x}.

Firth D (1993). Bias reduction of maximum likelihood estimates.
\emph{Biometrika}. \strong{80}, 27-38. \doi{10.2307/2336755}.

Kenne Pagui E C, Salvan A, Sartori N (2017). Median bias
reduction of maximum likelihood estimates. \emph{Biometrika}, \strong{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. \emph{Statistics and Computing},
\strong{30}, 43-59. \doi{10.1007/s11222-019-09860-6}.

Kosmidis I, Firth D (2009). Bias reduction in exponential family
nonlinear models. \emph{Biometrika}, \strong{96}, 793-804. \doi{10.1093/biomet/asp055}.
}
\author{
Euloge Clovis Kenne Pagui \verb{[aut]} \email{kenne@stat.unipd.it}, Ioannis Kosmidis \verb{[aut, cre]} \email{ioannis.kosmidis@warwick.ac.uk}
}