knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 6 )

The **brglm2** R package provides `brmultinom()`

which is a wrapper of `brglmFit`

for fitting multinomial logistic regression models (a.k.a. baseline category logit models) using either maximum likelihood or maximum penalized likelihood or any of the various bias reduction methods described in `brglmFit()`

. `brmultinom()`

uses the equivalent Poisson log-linear model, by appropriately re-scaling the Poisson means to match the multinomial totals (a.k.a. the "Poisson trick"). The mathematical details and algorithm on using the Poisson trick for mean-bias reduction are given in @kosmidis:11.

This vignettes illustrates the use of `brmultinom()`

and of the associated methods, using the alligator food choice example in @agresti:02[, Section 7.1]

The alligator data set ships with **brglm2**. @agresti:02[, Section 7.1] provides a detailed description of the variables recorded in the data set.

library("brglm2") data("alligators", package = "brglm2")

The following chunk of code reproduces @agresti:02[, Table 7.4]. Note that in order to get the estimates and standard errors reported in the latter table, we have to explicitly specify the contrasts that @agresti:02 uses.

agresti_contrasts <- list(lake = contr.treatment(levels(alligators$lake), base = 4), size = contr.treatment(levels(alligators$size), base = 2)) all_ml <- brmultinom(foodchoice ~ size + lake , weights = freq, data = alligators, contrasts = agresti_contrasts, ref = 1, type = "ML") all_ml_summary <- summary(all_ml) ## Estimated regression parameters round(all_ml_summary$coefficients, 2) ## Estimated standard errors round(all_ml_summary$standard.errors, 2)

Fitting the model using mean-bias reducing adjusted score equations gives

all_mean <- update(all_ml, type = "AS_mean") summary(all_mean)

The corresponding fit using median-bias reducing adjusted score equations is

all_median <- update(all_ml, type = "AS_median") summary(all_median)

The estimates and the estimated standard errors from bias reduction are close to those for maximum likelihood. As a result, it is unlikely that either mean or median bias is of any real consequence for this particular model and data combination.

Let's scale the frequencies in `alligators`

by 3 in order to get a sparser data set. The differences between maximum likelihood and mean and median bias reduction should be more apparent on the resulting data set. Here we have to "slow-down" the Fisher scoring iteration (by scaling the step-size), because otherwise the Fisher information matrix quickly gets numerically rank-deficient. The reason is data separation [@albert:84].

all_ml_sparse <- update(all_ml, weights = round(freq/3), slowit = 0.2) summary(all_ml_sparse)

Specifically, judging from the estimated standard errors, the estimates for `(Intercept)`

, `lakeHancock`

, `lakeOklawaha`

and `lakeTrafford`

for `Reptile`

and `lakeHancock`

for `Bird`

seem to be infinite.

To quickly check if that's indeed the case we can use the `check_infinite_estimates()`

method of the [**detectseparation**][https://cran.r-project.org/package=detectseparation] R package.

library("detectseparation") se_ratios <- check_infinite_estimates(all_ml_sparse) plot(se_ratios)

Some of the estimated standard errors diverge as the number of Fisher scoring iterations increases, which is evidence of complete or quasi-complete separation [@lesaffre:89].

In contrast, both mean and median bias reduction result in finite estimates

all_mean_sparse <- update(all_ml_sparse, type = "AS_mean") summary(all_mean_sparse) all_median_sparse <- update(all_ml_sparse, type = "AS_median") summary(all_median_sparse)

`?brglmFit`

and `?brglm_control`

contain quick descriptions of the various bias reduction methods supported in **brglm2**. The `iteration`

vignette describes the iteration and gives the mathematical details for the bias-reducing adjustments to the score functions for generalized linear models.

If you found this vignette or **brglm2**, in general, useful, please consider citing **brglm2** and the associated paper. You can find information on how to do this by typing `citation("brglm2")`

.

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.