1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
|
---
title: "Multinomial logistic regression using **brglm2**"
author: "[Ioannis Kosmidis](https://www.ikosmidis.com)"
date: "01 July 2017"
output: rmarkdown::html_vignette
bibliography: brglm2.bib
vignette: >
%\VignetteIndexEntry{Multinomial logistic regression using **brglm2**}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width = 6,
fig.height = 6
)
```
# **brmultinom**
The [**brglm2**](https://github.com/ikosmidis/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.
```{r, echo = TRUE}
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.
```{r, echo = TRUE}
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
```{r, echo = TRUE}
all_mean <- update(all_ml, type = "AS_mean")
summary(all_mean)
```
The corresponding fit using median-bias reducing adjusted score equations is
```{r, echo = TRUE}
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].
```{r, echo = TRUE, error = TRUE}
all_ml_sparse <- update(all_ml, weights = round(freq/3), slowit = 0.1)
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.
```{r, echo = TRUE}
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
```{r, echo = TRUE}
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`](https://cran.r-project.org/package=brglm2/brglm2.pdf) 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
|
