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 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 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 Agresti (2002, 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 (2002) 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)
#>              (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
## Estimated standard errors
round(all_ml_summary$standard.errors, 2)
#>              (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

all_mean <- update(all_ml, type = "AS_mean")
summary(all_mean)
#> 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 
#> Log-likelihood: -270.4071 
#> AIC: 580.8142

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

all_median <- update(all_ml, type = "AS_median")
summary(all_median)
#> 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 
#> Log-likelihood: -270.1223 
#> AIC: 580.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: 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 
#> Log-likelihood: -80.95651 
#> AIC: 201.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)
#> Warning: 'check_infinite_estimates' will be removed from 'brglm2' at version
#> 0.8. A new version of 'check_infinite_estimates' is now maintained in the
#> 'detectseparation' package.
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 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)
#> 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.53865701  -1.4524525   0.85785527    0.9583239
#> Reptile        -2.690142 -0.76143865   1.4034488   1.97464012    2.0975458
#> Bird           -1.820474 -0.25362479   0.7182343  -0.59602300    0.6373376
#> 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.5495132 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.0633024
#> 
#> Residual Deviance: 164.878 
#> Log-likelihood: -82.43901 
#> AIC: 204.878

all_median_sparse <- update(all_ml_sparse, type = "AS_median")
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.57656367  -1.6099398   0.87918613    0.9772025
#> Reptile        -3.665433 -0.89728837   2.1898215   2.82443510    2.9355873
#> Bird           -1.991952 -0.28581023   0.8094639  -1.37237185    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.4236128 1.1696370   2.5995648    2.5308829    2.5413660
#> Bird           1.0941192 1.0683124   1.2317627    2.5549130    1.4277852
#> Other          0.8695131 0.7700378   0.9449311    1.2711051    1.0676074
#> 
#> Residual Deviance: 162.9522 
#> Log-likelihood: -81.47608 
#> AIC: 202.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 and the associated paper. You can find information on how to do this by typing citation("brglm2").

References

Agresti, A. 2002. Categorical Data Analysis. 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): 109–16. http://www.jstor.org/stable/2345845.