Multinomial logistic regression using brglm2

Ioannis Kosmidis

01 July 2017

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 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 and Firth (2011).

This vignettes illustrates the use of brmultinom and of the associated methods, using the alligator food choice example in Agresti (2002, Section 7.1)

Alligator data

The alligator data set ships with brglm2. Agresti (2002, 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 fit@agresti:02[, Table 7.4]. Note that in order to get the estimates and standard errors reported in latter table, we have to explicitly specify the contrasts that Agresti (2002) uses.

## Warning in brmultinom(foodchoice ~ size + lake, weights = freq, data =
## alligators, : Observations with non-positive weights have been omited from
## the computations
##              (Intercept) size<=2.3 lakeHancock lakeOklawaha lakeTrafford
## Invertebrate       -1.55      1.46       -1.66         0.94         1.12
## Reptile            -3.31     -0.35        1.24         2.46         2.94
## Bird               -2.09     -0.63        0.70        -0.65         1.09
## Other              -1.90      0.33        0.83         0.01         1.52
##              (Intercept) size<=2.3 lakeHancock lakeOklawaha lakeTrafford
## Invertebrate        0.42      0.40        0.61         0.47         0.49
## Reptile             1.05      0.58        1.19         1.12         1.12
## Bird                0.66      0.64        0.78         1.20         0.84
## Other               0.53      0.45        0.56         0.78         0.62

Mean and median bias reduction

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

## Warning in brmultinom(formula = foodchoice ~ size + lake, data =
## alligators, : Observations with non-positive weights have been omited from
## the computations
## Call:
## brmultinom(formula = foodchoice ~ size + lake, data = alligators, 
##     weights = freq, contrasts = agresti_contrasts, ref = 1, type = "AS_mean")
## 
## Coefficients:
##              (Intercept)  size<=2.3 lakeHancock lakeOklawaha lakeTrafford
## Invertebrate   -1.490818  1.4019624  -1.5600445   0.90109748     1.081322
## Reptile        -2.908819 -0.3209487   0.9800182   2.10269205     2.559008
## Bird           -1.940190 -0.5828768   0.6213169  -0.41459453     1.026221
## Other          -1.815225  0.3147753   0.7792789   0.06001729     1.450717
## 
## Std. Errors:
##              (Intercept) size<=2.3 lakeHancock lakeOklawaha lakeTrafford
## Invertebrate   0.4250711 0.3962267   0.6003569    0.4742043    0.4924905
## Reptile        0.8864071 0.5565974   1.0218118    0.9619417    0.9604963
## Bird           0.6276892 0.6084642   0.7436454    1.0458118    0.8031649
## Other          0.5177641 0.4443760   0.5521476    0.7489610    0.6154619
## 
## Residual Deviance: 540.8142 
## AIC: 540.8142

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

## Warning in brmultinom(formula = foodchoice ~ size + lake, data =
## alligators, : Observations with non-positive weights have been omited from
## the computations
## Call:
## brmultinom(formula = foodchoice ~ size + lake, data = alligators, 
##     weights = freq, contrasts = agresti_contrasts, ref = 1, type = "AS_median")
## 
## Coefficients:
##              (Intercept)  size<=2.3 lakeHancock lakeOklawaha lakeTrafford
## Invertebrate   -1.508934  1.4134512  -1.6080461   0.90952486     1.090687
## Reptile        -3.128931 -0.3379387   1.1208787   2.28462590     2.743766
## Bird           -2.018328 -0.6036652   0.6554327  -0.55207076     1.044909
## Other          -1.854740  0.3201676   0.8002874   0.02330188     1.469440
## 
## Std. Errors:
##              (Intercept) size<=2.3 lakeHancock lakeOklawaha lakeTrafford
## Invertebrate   0.4238148 0.3949248   0.6066348    0.4720116    0.4901870
## Reptile        0.9729317 0.5699878   1.1068113    1.0427073    1.0411046
## Bird           0.6449317 0.6254237   0.7621085    1.1250548    0.8221610
## Other          0.5208087 0.4456384   0.5538846    0.7614315    0.6170972
## 
## Residual Deviance: 540.2445 
## AIC: 540.2445

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 and Anderson 1984).

all_ml_sparse <- update(all_ml, weights = round(freq/3), slowit = 0.2)
## Warning in brmultinom(formula = foodchoice ~ size + lake, data =
## alligators, : Observations with non-positive weights have been omited from
## the computations
## Warning: brglmFit: algorithm did not converge
summary(all_ml_sparse)
## Call:
## brmultinom(formula = foodchoice ~ size + lake, data = alligators, 
##     weights = round(freq/3), contrasts = agresti_contrasts, ref = 1, 
##     type = "ML", slowit = 0.2)
## 
## Coefficients:
##              (Intercept)   size<=2.3 lakeHancock lakeOklawaha lakeTrafford
## Invertebrate   -1.892117  1.74080751  -1.7936496   0.96047783    1.0747116
## Reptile       -22.030425 -1.05673342  20.4794861  21.17196942   21.3070760
## Bird           -2.225759 -0.33405988   0.9538709 -19.71099992    0.6957762
## Other          -1.730365  0.04554622   1.1091441  -0.07652818    1.2070845
## 
## Std. Errors:
##               (Intercept) size<=2.3  lakeHancock lakeOklawaha lakeTrafford
## Invertebrate 8.042604e-01 0.7409450 1.188368e+00 8.431835e-01 8.917234e-01
## Reptile      2.297088e+04 1.2808647 2.297088e+04 2.297088e+04 2.297088e+04
## Bird         1.184975e+00 1.1598538 1.324354e+00 2.486631e+04 1.545522e+00
## Other        8.869382e-01 0.7815002 9.589614e-01 1.337964e+00 1.084437e+00
## 
## Residual Deviance: 161.913 
## AIC: 161.913

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 (see, also the separation vignette).

se_ratios <- check_infinite_estimates(all_ml_sparse)
matplot(se_ratios, type = "l", lty = 1, ylim = c(0.5, 1.5), xlab = "Iteration")

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 and Albert 1989).

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

all_mean_sparse <- update(all_ml_sparse, type = "AS_mean")
## Warning in brmultinom(formula = foodchoice ~ size + lake, data =
## alligators, : Observations with non-positive weights have been omited from
## the computations
summary(all_mean_sparse)
## Call:
## brmultinom(formula = foodchoice ~ size + lake, data = alligators, 
##     weights = round(freq/3), contrasts = agresti_contrasts, ref = 1, 
##     type = "AS_mean", slowit = 0.2)
## 
## Coefficients:
##              (Intercept)   size<=2.3 lakeHancock lakeOklawaha lakeTrafford
## Invertebrate   -1.679168  1.53865702  -1.4524526   0.85785527    0.9583239
## Reptile        -2.690142 -0.76143882   1.4034489   1.97464011    2.0975458
## Bird           -1.820474 -0.25362492   0.7182343  -0.59602306    0.6373375
## Other          -1.517519  0.04996481   0.9501195   0.06500802    1.0685650
## 
## Std. Errors:
##              (Intercept) size<=2.3 lakeHancock lakeOklawaha lakeTrafford
## Invertebrate   0.7986896 0.7384928   1.0933484    0.8526692    0.9002978
## Reptile        1.5495133 1.0465737   1.7628485    1.7032986    1.7165517
## Bird           1.0387358 1.0051365   1.1808982    1.8023206    1.3550407
## Other          0.8606607 0.7662529   0.9397838    1.2305934    1.0633023
## 
## Residual Deviance: 164.878 
## AIC: 164.878
all_median_sparse <- update(all_ml_sparse, type = "AS_median")
## Warning in brmultinom(formula = foodchoice ~ size + lake, data =
## alligators, : Observations with non-positive weights have been omited from
## the computations
## Warning: brglmFit: algorithm did not converge
summary(all_median_sparse)
## Call:
## brmultinom(formula = foodchoice ~ size + lake, data = alligators, 
##     weights = round(freq/3), contrasts = agresti_contrasts, ref = 1, 
##     type = "AS_median", slowit = 0.2)
## 
## Coefficients:
##              (Intercept)   size<=2.3 lakeHancock lakeOklawaha lakeTrafford
## Invertebrate   -1.738085  1.57656366  -1.6099398   0.87918613    0.9772025
## Reptile        -3.665433 -0.89728841   2.1898217   2.82443525    2.9355875
## Bird           -1.991952 -0.28581026   0.8094639  -1.37237204    0.6272294
## Other          -1.613288  0.04720733   1.0087178  -0.01733333    1.1038823
## 
## Std. Errors:
##              (Intercept) size<=2.3 lakeHancock lakeOklawaha lakeTrafford
## Invertebrate   0.7922381 0.7310924   1.1339038    0.8411221    0.8888212
## Reptile        2.4236130 1.1696370   2.5995650    2.5308831    2.5413662
## Bird           1.0941192 1.0683125   1.2317627    2.5549132    1.4277852
## Other          0.8695131 0.7700378   0.9449311    1.2711051    1.0676074
## 
## Residual Deviance: 162.9522 
## AIC: 162.9522

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 using

citation("brglm2")
## Warning in citation("brglm2"): no date field in DESCRIPTION file of package
## 'brglm2'
## Warning in citation("brglm2"): could not determine year for 'brglm2' from
## package DESCRIPTION file
## 
## To cite package 'brglm2' in publications use:
## 
##   Ioannis Kosmidis (NA). brglm2: Bias Reduction in Generalized
##   Linear Models. R package version 0.1.7.
##   https://github.com/ikosmidis/brglm2
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {brglm2: Bias Reduction in Generalized Linear Models},
##     author = {Ioannis Kosmidis},
##     note = {R package version 0.1.7},
##     url = {https://github.com/ikosmidis/brglm2},
##   }

References

Agresti, Alan. 2002. Categorical Data Analysis. New York: Wiley.

Albert, A., and J. A. Anderson. 1984. “On the Existence of Maximum Likelihood Estimates in Logistic Regression Models.” Biometrika 71 (1):1–10.

Kosmidis, I., and D. Firth. 2011. “Multinomial Logit Bias Reduction via the Poisson Log-Linear Model.” Biometrika 98 (3):755–59.

Lesaffre, E., and A. Albert. 1989. “Partial Separation in Logistic Discrimination.” Journal of the Royal Statistical Society. Series B (Methodological) 51 (1). [Royal Statistical Society, Wiley]:109–16. http://www.jstor.org/stable/2345845.