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?
Best Answer
If you're using
rstan
or have access to R, I'd suggest loading thebrms
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 thatacat
is adjacent categories, and you may want one of the other ordinal model types (cumulative
,sratio
, orcratio
).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 casebrms
orrstanarm
would be very helpful. Of course, this may be too far afield from your original question. Just a suggestion.