Multilevel Ordered Logistic regression in Stan

Question

I'm trying to create a multilevel ordinal logistic regression model in Stan and the following converges:

   stanmodel <- '
data {
  int<lower=2> K;
  int<lower=0> N;
  int<lower=1,upper=K> y[N];
  int<lower=0> Ntests;
  int<lower=1,upper=Ntests> tests[N];
}

parameters {
  ordered[K-1]  CutpointsMean;
  real<lower=0> CutpointsSigma[K-1];
  ordered[K-1]  Cutpoints[Ntests];
}

model {

  CutpointsSigma ~ exponential(1);
  CutpointsMean  ~ normal(2, 3);

  for (i in 1:Ntests) {
    Cutpoints[i][1] ~ normal(CutpointsMean[1] , CutpointsSigma[1]);
    Cutpoints[i][2] ~ normal(CutpointsMean[2] , CutpointsSigma[2]);
    Cutpoints[i][3] ~ normal(CutpointsMean[3] , CutpointsSigma[3]);
  }


  for (i in 1:N)
    y[i] ~ ordered_logistic(0, Cutpoints[tests[i]]);

}

'CutpointsMean' and 'CutpointsSigma' define the population global ordinal response while Cutpoints[i][1:3] is the ordinal response for group i.

As Cutpoints[i] is an ordered vector of 3 elements, what happens when I write directly into Cutpoints[i][2] ? Is the write operation rejected if the constraints are not satisfied or simply the entry in written in the vector and and the results is sorted?

Is this the correct way of modelling a multilevel ordinal response in Stan?

Answer 1

If you're using rstan or have access to R, I'd suggest loading the brms package. Then create your model something like:

mod <- brm (tgt ~ foo + bar + baz, family="acat") stancode (mod)

And you'll see what brms does to implement such a model. Note that acat is adjacent categories, and you may want one of the other ordinal model types (cumulative, sratio, or cratio).

brms-generated Stan is fairly straightforward and readable. Or maybe you just want a Stan-based model but don't actually want to learn Stan itself, in which case brms or rstanarm would be very helpful. Of course, this may be too far afield from your original question. Just a suggestion.