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
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
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)
?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
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.