Adjacent category logit models using **brglm2**

In brglm2: Bias Reduction in Generalized Linear Models

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

bracl

The brglm2 R package provides bracl() which is a wrapper of brglmFit() for fitting adjacent category models for ordinal responses using either maximum likelihood or maximum penalized likelihood or any of the various bias reduction methods described in brglmFit(). There is a formal equivalence between adjacent category logit models for ordinal responses and multinomial logistic regression models (see, e.g. the Multinomial logistic regression using brglm2 vignette and the brmultinom() function). bracl() utilizes that equivalence and fits the corresponding 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.

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

Opinion on stem cell research and religious fundamentalism

The stemcell data set ships with brglm2. @agresti:15[, Section 4.1] provides a detailed description of the variables recorded in this data set (see also ?stemcell).

library("brglm2")
data("stemcell", package = "brglm2")
stem <- within(stemcell, religion <- as.numeric(religion))

Maximum likelihood estimation

The following chunk of code fits an adjacent category logit model with proportional odds and reproduces @agresti:10[, Table 4.2]. Note that the intercept parameters are different because @agresti:10[, Table 4.2] uses different contrasts for the intercept parameters.

stem_formula <- research ~ religion + gender
stemcells_ml <- bracl(stem_formula, weights = frequency, data = stem,
                      parallel = TRUE, type = "ML")
summary(stemcells_ml)

stemcells_ml is an object inheriting from

class(stemcells_ml)

brglm2 implements print, coef, fitted, predict, summary, vcov and logLik methods for

We can check if a model with non-proportional odds fits the data equally well by fitting it and carrying out a likelihood ration test.

stemcells_ml_full <- bracl(stem_formula, weights = frequency, data = stemcell,
                           parallel = FALSE, type = "ML")
summary(stemcells_ml_full)

The value of the log likelihood ratio statistic here is

(lrt <- deviance(stemcells_ml) - deviance(stemcells_ml_full))

and has an asymptotic chi-squared distribution with

(df1 <- df.residual(stemcells_ml) - df.residual(stemcells_ml_full))

The p-value from testing the hypothesis that stemcells_ml_full is an as good fit as stemcells_ml is

pchisq(lrt, df1, lower.tail = FALSE)

hence, the simpler model is found to be as adequate as the full model is.

Mean and median bias reduction

We can use bracl() to fit the adjacent category model using estimators with smaller mean or median bias. For mean bias reduction we do

summary(update(stemcells_ml, type = "AS_mean"))

and for median

summary(update(stemcells_ml, type = "AS_median"))

The estimates from mean and median bias reduction are similar to the maximum likelihood ones, indicating that estimation bias is not a major issue here.

Prediction

We can predict the category probabilities using the predict() method

predict(stemcells_ml, type = "probs")

Relevant resources

?brglmFit and ?brglm_control provide 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

brglm2 documentation built on Oct. 12, 2023, 1:07 a.m.