So, BRMS, in which the author has coded a simpler alternative that requires the estimation of just two parameters, irrespective of the number of ordinal X categories. As I understand the paper (*Modeling Monotonic effects of ordinal predictors*), one parameter controls the size of the distance between the first and last categories, and the second parameter is a simplex, which sums to unity, and represents the uneven distances between the ordinal categories.

So here I want to find the part of the code that generates the coding of the two parameters. I found this brms/brmsformula.R at master · paul-buerkner/brms · GitHub but I want to confirm whether I am going in the right direction or not?

Please, can somebody help?

Welcome to The Stan Forums, Rohit!

Do you mean that you want to find the brms (R) code that generates the Stan code corresponding to those two parameters? If yes, then I’m sorry because I can’t help (I would have to search through the brms code in the same way as you).

1 Like

Yes, you are right. I tried but I am not sure if it is right or not. Actually, I have not too much experience with R. I am mostly involved in Python & Julia. You and @paul.buerkner were the top contributor to the package so I thought you will have a better idea than I.

Thanks for the welcome!

The approach that you’ve mentioned (modelling monotonic effects of ordinal predictors) is implemented in `brms`

using the `mo()`

functions for `brms`

formulas, with more implementation details outlined in this vignette: Estimating Monotonic Effects with brms

The `brms`

source code responsible for generating the Stan code to handle monotonic effects is held within the `data_sp`

(preparing/coding input data) and `stan_sp`

(Stan code generation) functions.

To see the resulting Stan code from using the `mo()`

framework, you can call `make_stancode()`

with a target `brm`

model. Using one of the examples from the vignette above, this would look like:

```
# Create example data
income_options <- c("below_20", "20_to_40", "40_to_100", "greater_100")
income <- factor(sample(income_options, 100, TRUE),
levels = income_options, ordered = TRUE)
mean_ls <- c(30, 60, 70, 75)
ls <- mean_ls[income] + rnorm(100, sd = 7)
dat <- data.frame(income, ls)
# Request generated Stan code:
make_stancode(ls ~ mo(income), data = dat)
```

Which gives:

```
// generated with brms 2.17.0
functions {
/* compute monotonic effects
* Args:
* scale: a simplex parameter
* i: index to sum over the simplex
* Returns:
* a scalar between 0 and 1
*/
real mo(vector scale, int i) {
if (i == 0) {
return 0;
} else {
return rows(scale) * sum(scale[1:i]);
}
}
}
data {
int<lower=1> N; // total number of observations
vector[N] Y; // response variable
int<lower=1> Ksp; // number of special effects terms
int<lower=1> Imo; // number of monotonic variables
int<lower=1> Jmo[Imo]; // length of simplexes
int Xmo_1[N]; // monotonic variable
vector[Jmo[1]] con_simo_1; // prior concentration of monotonic simplex
int prior_only; // should the likelihood be ignored?
}
transformed data {
}
parameters {
real Intercept; // temporary intercept for centered predictors
simplex[Jmo[1]] simo_1; // monotonic simplex
vector[Ksp] bsp; // special effects coefficients
real<lower=0> sigma; // dispersion parameter
}
transformed parameters {
real lprior = 0; // prior contributions to the log posterior
lprior += student_t_lpdf(Intercept | 3, 63.6, 15.7);
lprior += dirichlet_lpdf(simo_1 | con_simo_1);
lprior += student_t_lpdf(sigma | 3, 0, 15.7)
- 1 * student_t_lccdf(0 | 3, 0, 15.7);
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N] mu = Intercept + rep_vector(0.0, N);
for (n in 1:N) {
// add more terms to the linear predictor
mu[n] += (bsp[1]) * mo(simo_1, Xmo_1[n]);
}
target += normal_lpdf(Y | mu, sigma);
}
// priors including constants
target += lprior;
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept;
}
```

2 Likes

Thanks for your answer @andrjohns.

- The mo() function is dependent on brmsformula ? And, the link which I have attached for the
`mo()`

is correct?
- And, we can use all the brms formulas for that
`mo()`

function?

`mo()`

is simply a helper function to be used within the `formula`

argument of a `brmsformula`

.

See the vignette I linked above for more information around how it can be used

1 Like

Thanks a lot! I appreciate your support.