% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/brmultinom.R
\name{brmultinom}
\alias{brmultinom}
\title{Bias reduction for multinomial response models using the
Poisson trick.}
\usage{
brmultinom(
  formula,
  data,
  weights,
  subset,
  na.action,
  contrasts = NULL,
  ref = 1,
  model = TRUE,
  x = TRUE,
  control = list(...),
  ...
)
}
\arguments{
\item{formula}{
a formula expression as for regression models, of the form
\code{response ~ predictors}. The response should be a factor or a
matrix with K columns, which will be interpreted as counts for each of
K classes.
A log-linear model is fitted, with coefficients zero for the first
class. An offset can be included: it should be a numeric matrix with K columns
if the response is either a matrix with K columns or a factor with K >= 2
classes, or a numeric vector for a response factor with 2 levels.
See the documentation of \code{\link{formula}()} for other details.
}

\item{data}{
an optional data frame in which to interpret the variables occurring
in \code{formula}.
}

\item{weights}{
optional case weights in fitting.
}

\item{subset}{
expression saying which subset of the rows of the data should  be used
in the fit. All observations are included by default.
}

\item{na.action}{
a function to filter missing data.
}

\item{contrasts}{
a list of contrasts to be used for some or all of
the factors appearing as variables in the model formula.
}

\item{ref}{the reference category to use for multinomial
regression. Either an integer, in which case
\code{levels(response)[ref]} is used as a baseline, or a character
string. Default is 1.}

\item{model}{
  logical. If true, the model frame is saved as component \code{model}
  of the returned object.
}

\item{x}{should the model matrix be included with in the result
(default is \code{TRUE}).}

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

\item{...}{arguments to be used to form the default \code{control}
argument if it is not supplied directly.}
}
\description{
\code{\link[=brmultinom]{brmultinom()}} is a wrapper of \code{\link[=brglmFit]{brglmFit()}} that fits multinomial
regression models using implicit and explicit bias reduction
methods. See Kosmidis & Firth (2011) for details.
}
\details{
The models \code{\link[=brmultinom]{brmultinom()}} handles are also known as
baseline-category logit models (see, Agresti, 2002, Section 7.1),
because they model the log-odds of every category against a
baseline category. The user can control which baseline (or
reference) category is used via the \code{ref}. By default
\code{\link[=brmultinom]{brmultinom()}} uses the first category as reference.

The maximum likelihood estimates for the parameters of
baseline-category logit models have infinite components with
positive probability, which can result in problems in their
estimation and the use of inferential procedures (e.g. Wad
tests). Albert and Andreson (1984) have categorized the possible
data patterns for such models into the exclusive and exhaustive
categories of complete separation, quasi-complete separation and
overlap, and showed that infinite maximum likelihood estimates
result when complete or quasi-complete separation occurs.

The adjusted score approaches to bias reduction that \code{\link[=brmultinom]{brmultinom()}}
implements for \code{type = "AS_mean"} and \code{type = "AS_median"} are
alternatives to maximum likelihood that result in estimates with
smaller asymptotic mean and median bias, respectively, that are
also \emph{always} finite, even in cases of complete or quasi-complete
separation.

\code{\link[=brmultinom]{brmultinom()}} is a wrapper of \code{\link[=brglmFit]{brglmFit()}} that fits multinomial
logit regression models through the 'Poisson trick' (see, for
example, Palmgren, 1981; Kosmidis & Firth, 2011).

The implementation relies on the construction of an extended model
matrix for the log-linear model and constraints on the sums of the
Poisson means. Specifically, a log-linear model is fitted on a
\href{https://en.wikipedia.org/wiki/Kronecker_product}{Kronecker product} of the
original model matrix \code{X} implied by the formula, augmented by
\code{nrow(X)} dummy variables.

The extended model matrix is sparse, and the
\href{https://cran.r-project.org/package=Matrix}{\pkg{Matrix}} package is
used for its effective storage.

While \code{\link[=brmultinom]{brmultinom()}} can be used for analyses using multinomial
regression models, the current implementation is more of a proof of
concept and is not expected to scale well with either of \code{nrow(X)},
\code{ncol(X)} or the number of levels in the categorical response.
}
\examples{
## The housing data analysis from ?MASS::housing

data("housing", package = "MASS")

# Maximum likelihood using nnet::multinom
houseML1nnet <- nnet::multinom(Sat ~ Infl + Type + Cont, weights = Freq,
                               data = housing)
# Maximum likelihood using brmultinom with baseline category 'Low'
houseML1 <- brmultinom(Sat ~ Infl + Type + Cont, weights = Freq,
                       data = housing, type = "ML", ref = 1)
# The estimates are numerically the same as houseML0
all.equal(coef(houseML1nnet), coef(houseML1), tolerance = 1e-04)

# Maximum likelihood using brmultinom with 'High' as baseline
houseML3 <- brmultinom(Sat ~ Infl + Type + Cont, weights = Freq,
                      data = housing, type = "ML", ref = 3)
# The fitted values are the same as houseML1
all.equal(fitted(houseML3), fitted(houseML1), tolerance = 1e-10)

# Bias reduction
houseBR3 <- update(houseML3, type = "AS_mean")
# Bias correction
houseBC3 <- update(houseML3, type = "correction")

## Reproducing Bull et al. (2002, Table 3)
data("hepatitis", package = "brglm2")

# Construct a variable with the multinomial categories according to
# the HCV and nonABC columns
hepat <- hepatitis
hepat$type <- with(hepat, factor(1 - HCV * nonABC + HCV + 2 * nonABC))
hepat$type <- factor(hepat$type, labels = c("noDisease", "C", "nonABC"))
contrasts(hepat$type) <- contr.treatment(3, base = 1)

# Maximum likelihood estimation fails to converge because some estimates are infinite
hepML <- brmultinom(type ~ group * time, data = hepat, weights = counts, type = "ML", slowit = 0.1)

# Mean bias reduction returns finite estimates
hep_meanBR <- brmultinom(type ~ group * time, data = hepat, weights = counts, type = "AS_mean")
# The estimates in Bull et al. (2002, Table 3, DOI: 10.1016/S0167-9473(01)00048-2)
coef(hep_meanBR)

# Median bias reduction also returns finite estimates, which are a bit larger in absolute value
hep_medianBR <- brmultinom(type ~ group * time, data = hepat, weights = counts, type = "AS_median")
coef(hep_medianBR)

}
\references{
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}.

Agresti A (2002). \emph{Categorical data analysis} (2nd edition). Wiley
Series in Probability and Statistics. Wiley.

Albert A, Anderson J A (1984). On the Existence of Maximum
Likelihood Estimates in Logistic Regression Models. \emph{Biometrika},
\strong{71} 1--10. \doi{10.2307/2336390}.

Kosmidis I, Firth D (2011). Multinomial logit bias reduction
via the Poisson log-linear model. \emph{Biometrika}, \strong{98}, 755-759.
\doi{10.1093/biomet/asr026}.

Palmgren, J (1981). The Fisher Information Matrix for Log Linear
Models Arguing Conditionally on Observed Explanatory
Variables. \emph{Biometrika}, \strong{68}, 563-566.
\doi{10.1093/biomet/68.2.563}.
}
\seealso{
\code{\link[nnet:multinom]{nnet::multinom()}}, \code{\link[=bracl]{bracl()}} for adjacent category logit models for ordinal responses
}
\author{
Ioannis Kosmidis \verb{[aut, cre]} \email{ioannis.kosmidis@warwick.ac.uk}
}