# Multinomial logistic regression using **brglm2** In brglm2: Bias Reduction in Generalized Linear Models

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

# brmultinom

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]

# Alligator data

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")
```

## Maximum likelihood estimation

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)
```

## Mean and median bias reduction

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.

# Infinite estimates and multinomial logistic regression

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)
```

# Relevant resources

`?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.

# Citation

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")`.

# References

## Try the brglm2 package in your browser

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

brglm2 documentation built on Nov. 21, 2021, 5:08 p.m.