Proportional odds assumptions

Interfaces brms
1

I have used brms package to fit an ordered multilevel logistic model (shown below). I would wish to determine whether my model is still valid by testing proportional odds assumption. However, despite reading the brms manual repeatedly, I don’t seem to find out how to test this critical assumption.

polr.model<- brm(orderedScore ~ comobidities
+ cadre
+ severity.levels
+ cadre.gender
+ child.sex
+ age.rescaled
+ (cadre | hospital/clinician),
family = “cumulative”
, prior = c(set_prior(“normal(0,3)”, class=“Intercept”)
, set_prior(“normal(0,3)”, class=“b”), set_prior(“cauchy(0,2)”, class=“sd”))
, control = list(adapt_delta = 0.99)
, data=dataToUSe
, cores=parallel::detectCores()
)

Kindly provide guidance.
Thank,
Morris

  • Operating System: linux
  • brms Version: 2.4.2
2

A general way to test an assumption is to fit a second model that comes without it and compare the two models against each other. For the proportional odds assumption, this comes down to fitting category specific effects for the predictors, which you may achieve by wrapping predictors in cs(), for instance cs(age.rescaled).

However, the cumulative model with category specific effects (i.e. the non- or partial- proportional odds model) is ill defined as it may imply negative probabilities. Thus, you need to use a different ordinal family, for instance "sratio" (sequential model) or "acat" (adjacent category model). You can read more about all of this in my tutorial paper about ordinal models: https://psyarxiv.com/x8swp/

1 Like
3

Many thanks for the quick reply Paul.
I have read the materials you referred me to and i very much appreciate too. A quick question, on the assessment of the proportional odds assumptions, is there a published paper/book one can cite about this proposed method of comparing category specific effect model with cumulative model?

4

You may very well cite the paper I referenced in the above post.

5

I have read (ordinal regression and IRT with BRMS) tutorials and checked stan/brms forum for a solid example for testing the PO assumption (with cumulative link). I am hoping to get some feedback towards forming a working example for the application of the available information about PO assumption testing with brms. In particular I will appreciate feedback as to how to compare models properly , and how a comparison of two models of different families (cumulative and adjacent cats.) makes sense.

Here is my baseline version of the example (subject to corrections):The data-set includes multiple person/item (paired) measurements (ordered categories) before and after an intervention.
The variables of the model:

  • pl: response variable-ordinal-3 point likert item
  • time: categorical (before,after)
  • Persoon: categorical
  • compItem: categorical

Cumulative link with PO assumption is the best model that suits the ordering relation between categories.

fitc <- brm(pl ~ time+(1|Persoon)+(1|compItem), data = mydata, prior= prior_ma, family=cumulative(link= “logit”, threshold=“flexible”))

I have created another model, this time with familiy =acat and with CS effects for comparison, hence to check the differences between the the cumulative and acat models:

fita2 <- brm(pl ~ cs(time)+(1|Persoon)+(1|compItem), data = mydata, prior= prior_ma, family=acat(link= “logit”, threshold=“flexible”))

The comparison of the two models against each other does not suggest that addition of CS effects makes an important difference (?), so I can assume PO assumption holds.

Output 1 (Cumulative)

Family: cumulative
Links: mu = logit; disc = identity
Formula: pl ~ time + (1 | Persoon) + (1 | compItem)
Data: assdata (Number of observations: 1224)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000

Group-Level Effects:
~compItem (Number of levels: 52)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 1.18 0.14 0.93 1.49 1057 1.00

~Persoon (Number of levels: 16)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 1.57 0.35 1.02 2.37 929 1.00

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept[1] -1.58 0.43 -2.48 -0.74 540 1.01
Intercept[2] 1.03 0.43 0.14 1.90 549 1.01
timet2 0.60 0.12 0.37 0.84 5682 1.00

Output 2 (Acat with CS effects)
Family: acat
Links: mu = logit; disc = identity
Formula: pl ~ cs(time) + (1 | Persoon) + (1 | compItem)
Data: assdata (Number of observations: 1224)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000

Group-Level Effects:
~compItem (Number of levels: 52)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 0.95 0.12 0.74 1.20 1282 1.00

~Persoon (Number of levels: 16)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept) 1.24 0.29 0.80 1.89 1276 1.01

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept[1] -1.10 0.36 -1.82 -0.40 682 1.01
Intercept[2] 0.69 0.36 -0.02 1.41 731 1.00
timet2[1] 0.54 0.16 0.23 0.86 8089 1.00
timet2[2] 0.52 0.16 0.22 0.82 7964 1.00

6

For a fair comparison, you also need to fit the acat() models without CS-effects as explained in the ordinal tutorial paper. Otherwise you may conflate the distribution with the CS assumption.

2 Likes
7

Thank you for your timely response and sorry for my late reply. Now I understand the approach and also run the following for prop-odds assumption:

fita1  <- brm(pl ~ time+(1|Persoon)+(1|compItem), data = mydata, 
  prior= prior_ma, family=acat(link= “logit”, threshold=“flexible”))

I have one related question :
I was able to plot population marginal_effects(fitc, categorical = TRUE) .Any recommendations for plotting crossed level effects for a category levels (e.g., for the levels of factor Persoon). ?

I am exploring brms API and I learn alot from it, brms+stan+bayesian is addictive!

8

Can you please explain what exactly you mean by cross-level effects in the context of your model?

9

Note that you can use triple backticks (```) to mark code blocks in yout posts. I edited the last post for you to show this. (I have little to say about your particular model, sorry).

Best of luck!

1 Like
10

Hi Paul, for example, I want to plot posterior mean estimate of the probability of ordinal responses for each compItem category e.g., a plot like marginal_effects populated for each compItem category.

11

You need to use the conditions argument and set re_formula = NULL. See the examples section of ?marginal_effects.

2 Likes