# Stan: Multilevel Ordinal Logistic Regression

I'm trying to create a multilevel ordinal logistic regression model in Stan and the following code would seem to work, in the sense that Stan seems to convergence to sensible answers:

stanmodel <- '
data {
int<lower=2> K;  // ordinal response with 4 values, 3 cutpoints
int<lower=0> N;  // number of measurements
int<lower=1,upper=K> y[N]; // response

int<lower=0> Ntests;         // number of groups
int<lower=1,upper=Ntests> tests[N];  // groups
}

parameters {
// population cutpoints and associated variance.
ordered[K-1]  CutpointsMean;
real<lower=0> CutpointsSigma[K-1];

ordered[K-1]  Cutpoints[Ntests];   // ordinal response cutpoints for groups
}

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]]);

}
'


I have removed the part relating to the covariates for clarity.

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

The idea is that for example the first 'cutpoint' of each group is generated from a normal distribution centered on the first 'cutpoint' of the overall population.

The second 'cutpoint' of each group is generated from a normal distribution centered on the second 'cutpoint' of the overall population and so on.

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 is written in the vector and the results is sorted?

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

• You can't write to the elements of something declared in the parameters block; you will get a parser error. What you have looks like legal syntax, but I don't think a normal prior for the cutpoints makes much sense. At a minimum, you need a truncated normal. – Ben Goodrich Jul 31 '16 at 18:29
• HI Ben, thanks for your reply. I have also created another model where the first cutpoint is normal and the other two cutpoints are half normal. The cutpoints are then recovered by adding the offsets back together: ordererd_cuts[1] <- cuts[1]; ordererd_cuts[2] <- cuts[1] + cuts[2]; ordererd_cuts[3] <- cuts[1] + cuts[2] + cuts[3]; where cuts[2] and cuts[3] are positive half normal – Anton Aug 1 '16 at 11:49
• Try stan_polr from rstanarm! – Zach May 2 '18 at 1:06