Looking into the part of the code that generates the coding of the two parameters

Interfaces brms
1

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?

2

@paul.buerkner @fweber144 Please help me! :(

3

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
4

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!

5

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
6

Thanks for your answer @andrjohns.

  1. The mo() function is dependent on brmsformula ? And, the link which I have attached for the mo() is correct?
  2. And, we can use all the brms formulas for that mo() function?
7

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
8

Thanks a lot! I appreciate your support.