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

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.

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

.

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

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.

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.

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

method

predict(stemcells_ml, type = "probs")

`?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.

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

.

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.