Adjacent category probit ordinal models expected predictions

Interfaces
1

I am unsure how to get the expected predictions from an adjacent category probit ordinal model with a categorical predictor and category specific effects. @paul.buerkner

Here’s the code of the model:

library(tudyverse)
library(brms)

acceptability <- read_csv("https://osf.io/download/xmkvg/")

acc_bm_2 <- brm(
  Response ~ cs(Condition) + (Condition | Participant) + (1 | Item),
  data = acceptability,
  family = acat("probit"),
  cores = 4,
  seed = 7823
)

Response has 7 categories, Condition is canonical or noncanonical.

Here’s the output of fixef():

                           Estimate Est.Error       Q2.5      Q97.5
Intercept[1]             -1.4340338 0.4908987 -2.4498975 -0.5079088
Intercept[2]             -0.6190555 0.3761126 -1.3387427  0.1205419
Intercept[3]             -1.6452514 0.2921813 -2.2186126 -1.0948489
Intercept[4]             -0.8509881 0.1924520 -1.2247104 -0.4828520
Intercept[5]             -1.1222150 0.1475053 -1.4220528 -0.8414859
Intercept[6]             -0.4568727 0.1134987 -0.6824405 -0.2305248
Conditionnoncanonical[1] -0.1406638 0.5366772 -1.2352887  0.8659495
Conditionnoncanonical[2] -0.1425097 0.4109545 -0.9254142  0.6671332
Conditionnoncanonical[3] -0.7785740 0.3209673 -1.4098297 -0.1599417
Conditionnoncanonical[4] -0.1738039 0.2078171 -0.5804308  0.2391148
Conditionnoncanonical[5] -0.6850219 0.1469132 -0.9824149 -0.4053308
Conditionnoncanonical[6] -0.3287998 0.1021691 -0.5301052 -0.1251982

I think to get the mean expected probability of category 6 when Condition = canonical (the reference level), I can do (?):

as_draws_df(acc_bm_2) |>
  mutate(e6 = pnorm(`b_Intercept[6]`)) |>
  summarise(mean(e6))

But for Condition = noncanonical the following doesn’t seem to be correct.

as_draws_df(acc_bm_2) |>
  mutate(e6 = pnorm(`b_Intercept[6]` - `bcs_Conditionnoncanonical[6]`)) |>
  summarise(mean(e6))
> # A tibble: 1 × 1
>  `mean(e6)`
>       <dbl>
> 1      0.449

This is the output of conditional_effects() for reference.

Screenshot 2025-01-08 at 15.57.07
Screenshot 2025-01-08 at 15.57.071360×842 34 KB

I am clearly doing that wrong, but not sure what the right way is.

2

Yeah, the acat model is a bit weird in that regard. you can read more about how its probabilities are computed in https://osf.io/preprints/psyarxiv/x8swp

Also, you can trust that brms (hopefully) computes probabilities correctly via posterior_epred (or fitted if you directly want the posterior summaries).

3

Thanks Paul! I did indeed read through that paper and I gathered that Eq (13) is the answer, but I might have misunderstood it.

\Pr(Y = k \mid Y \in \{k, k+1\}, \eta) = F(\tau_k - \eta)

If F is \Phi i.e. the CDF of the standard Gaussian distribution and \eta is the linear predictor which is the sum of the regression coefficients, then shouldn’t the following work?

as_draws_df(acc_bm_2) |>
  mutate(e6 = pnorm(`b_Intercept[6]` - `bcs_Conditionnoncanonical[6]`)) |>
  summarise(mean(e6))

It does seem to work when I don’t use category specific effects (i.e. when there is only one coefficient for the predictor).

It’s probably my lack of mathematical understanding, so forgive me if the question is trivial!

4

The problem is that what the equation gives you is \Pr(Y = k \mid Y \in \{k, k+1\}, \eta) but what I believe you want is \Pr(Y = k \mid \eta) that is without the condition. In other words, you should check a different equation in the paper.

5

That’s why you (me) should read the appendix… So, if I got it right, the equation is (A22)

P(Y = k \mid \eta) = \frac{\exp\left(\sum_{j=1}^{k-1} (\eta - \tau_j)\right)}{\sum_{r=1}^{K+1} \exp\left(\sum_{j=1}^{r-1} (\eta - \tau_j)\right)}

Thanks!

1 Like
6

Yes, exactly! :-)

1 Like