Item Response Model definition in brms for cumulative logit link

In his article on item response theory in brms (, Paul shows how to use his package for a so-called two-parameter graded response model. These are cumulative logit models with two item-specific parameters. An easiness parameter, say b_i, and a discrimination parameter, say a_i, for each item i.

A simple version of my current code to fit such a model would be

family <- brmsfamily("cumulative", "logit")
formula <- bf(score ~ 1 + (1 | item) + (1 | rater),  disc  ~ 1 + (1 | item))
model <- brm(formula = formula, data = df, family = family)

However, mathematically, I can see two possible models that could be defined by this syntax. Let’s assume score has three ordinal categories Y_{i,j} \in \{1,2,3\} and that raters are enumerated by j = 1, ..., J and their ability is \theta_j. The probability to score a 1 (lowest category) could be defined in two ways now.

Version 1

Pr(Y_{i,j} = 1) = \text{expit}(\beta_1 - ((a_i) * (b_i + \theta_j)))

Version 2

Pr(Y_{i,j} = 1) = \text{expit}(a_i (\beta_1 - (b_i + \theta_j)))

In the first, the threshold parameter \beta_1 is not multiplied by the discrimination parameter, in the second one, it is. Of course, it’s the same for the definition of the probability of the second category, but this is enough to make the point.

To me, it is not clear which of the two models brms is fitting here, even after carefully studying the associated paper. It seems they would refer to different underlying model definitions?

I have found the solution in another paper by Paul (Ordinal Regression Models in Psychology: A Tutorial, equation 36), it is Version 2. It can also be seen in the Stan code

  real cumulative_logit_lpmf(int y, real mu, real disc, vector thres) {
    int nthres = num_elements(thres);
    if (y == 1) {
      return log_inv_logit(disc * (thres[1] - mu));
    } else if (y == nthres + 1) {
      return log1m_inv_logit(disc * (thres[nthres] - mu));
    } else {
      return log_diff_exp(log_inv_logit(disc * (thres[y] - mu)),
                          log_inv_logit(disc * (thres[y - 1] - mu)));