# Hierarchical ordered logit model in stan and brms

I want to model the hierarchical ordered logit model in stan. And I only want the threshold to vary with the group variables. I am not sure if there is an example or vignette in the hierarchical ordered logit model
First I use the brms code like this

`````` brm(data = data,
family = cumulative(),
y|thres(gr=yr) ~ x1+x2+x3,
iter = 4000, warmup = 2000, chains = 4, cores = 4,
seed = 23
)
``````

The `yr` is group variable means I have 3 different groups.
`x1,x2,x3`means predicted variables
And I use Stan to model a threshold varying hierarchical model like this

``````functions {
real induced_dirichlet_lpdf(vector c, vector alpha, real phi) {
int K = num_elements(c) + 1;
vector[K - 1] sigma = inv_logit(phi - c);
vector[K] p;
matrix[K, K] J = rep_matrix(0, K, K);

// Induced ordinal probabilities
p[1] = 1 - sigma[1];
for (k in 2:(K - 1))
p[k] = sigma[k - 1] - sigma[k];
p[K] = sigma[K - 1];

// Baseline column of Jacobian
for (k in 1:K) J[k, 1] = 1;

// Diagonal entries of Jacobian
for (k in 2:K) {
real rho = sigma[k - 1] * (1 - sigma[k - 1]);
J[k, k] = - rho;
J[k - 1, k] = rho;
}

return   dirichlet_lpdf(p | alpha)
+ log_determinant(J);
}
}
pper=K> y[N]; // Observed ordinals
row_vector[D] x[N];
int g[N];
}

parameters {
vector[D] beta;       // Latent effect
ordered[K - 1] c; // (Internal) cut points
real a[M];
real<lower=0,upper=10> sigma;
}

model {
// Prior model
//gamma ~ normal(0, 1);
a ~ normal(0,sigma);
c ~ induced_dirichlet(rep_vector(1, K), 0);

// Observational model
for (n in 1:N){
y[n] ~ ordered_logistic(x[n]* beta + a[g[n]], c);
}

}

generated quantities {
int<lower=1, upper=K> y_ppc[N];
for (n in 1:N)
y_ppc[n] = ordered_logistic_rng(x[n,]*beta, c);
}

``````

I use awkward coding to separate data to 3 groups and input to stan.
I have two questions : 1. what is the different with brms code and stan in my hierarchical model
2. if I did wrong in writing this stan code for my model

I don’t speak brms, so it’s very hard to tell without seeing a mathematical definition of what you are trying to achieve.

May I ask why you’re writing your own “induced Dirichlet” distribution and what it’s supposed to be doing?

I’m the opposite of Bob in that I speak brms, but not raw Stan. The left-hand side of your model formula

indicates your thresholds will vary across your group variable `yr`. However, I don’t see anything in your formula to indicate your model is hierarchical.

And I edit my code to a new concept

``````functions {
real induced_dirichlet_lpdf(vector c, vector alpha, real phi) {
int K = num_elements(c) + 1;
vector[K - 1] sigma = inv_logit(phi - c);
vector[K] p;
matrix[K, K] J = rep_matrix(0, K, K);

// Induced ordinal probabilities
p[1] = 1 - sigma[1];
for (k in 2:(K - 1))
p[k] = sigma[k - 1] - sigma[k];
p[K] = sigma[K - 1];

// Baseline column of Jacobian
for (k in 1:K) J[k, 1] = 1;

// Diagonal entries of Jacobian
for (k in 2:K) {
real rho = sigma[k - 1] * (1 - sigma[k - 1]);
J[k, k] = - rho;
J[k - 1, k] = rho;
}

return   dirichlet_lpdf(p | alpha)
+ log_determinant(J);
}
}

data {
int<lower=2> N1;             // Number of observations
int<lower=2> N2;
int<lower=2> N3;
int<lower=1> N4;
int<lower=1> N5;
int<lower=1> N6;
int<lower=1> N7;
int<lower=1> N8;
int<lower=1> N9;
int<lower=1> N10;
int<lower=1> N11;
int<lower=1> N12;
int<lower=1> N13;
int<lower=2> K;             // Number of ordinal categories
int<lower=2> D;
int<lower=2> T;
int<lower=1, upper=K> y1[N1]; // Observed ordinals
int<lower=1, upper=K> y2[N2];
int<lower=1, upper=K> y3[N3];
int<lower=1, upper=K> y4[N4];
int<lower=1, upper=K> y5[N5];
int<lower=1, upper=K> y6[N6];
int<lower=1, upper=K> y7[N7];
int<lower=1, upper=K> y8[N8];
int<lower=1, upper=K> y9[N9];
int<lower=1, upper=K> y10[N10];
int<lower=1, upper=K> y11[N11];
int<lower=1, upper=K> y12[N12];
int<lower=1, upper=K> y13[N13];

matrix[N1,D] x1;
matrix[N2,D] x2;
matrix[N3,D] x3;
matrix[N4,D] x4;
matrix[N5,D] x5;
matrix[N6,D] x6;
matrix[N7,D] x7;
matrix[N8,D] x8;
matrix[N9,D] x9;
matrix[N10,D] x10;
matrix[N11,D] x11;
matrix[N12,D] x12;
matrix[N13,D] x13;
}

parameters {

vector[D] beta;       // Latent effect

real a[K-1];
real<lower=0,upper=10> sigma[K-1];
vector[13] omega[K-1];
}
transformed parameters{
matrix[K-1,13] mu;
mu[1,] = exp(a[1]+ sigma[1] * omega[1]');
for (i in 2: K-1){
mu[i,] = mu[i-1,] +exp(a[i]+ sigma[i] * omega[i]');
}

}
model {

for (i in 1: (K-1)){
omega[i] ~ normal(0,1);
}

// Observational model
y1 ~ ordered_logistic(x1* beta, mu[,1]);
y2 ~ ordered_logistic(x2* beta, mu[,2]);
y3 ~ ordered_logistic(x3* beta, mu[,3]);
y4 ~ ordered_logistic(x4* beta, mu[,4]);
y5 ~ ordered_logistic(x5* beta, mu[,5]);
y6 ~ ordered_logistic(x6* beta, mu[,6]);
y7 ~ ordered_logistic(x7* beta, mu[,7]);
y8 ~ ordered_logistic(x8* beta, mu[,8]);
y9 ~ ordered_logistic(x9* beta, mu[,9]);
y10 ~ ordered_logistic(x10* beta, mu[,10]);
y11 ~ ordered_logistic(x11* beta, mu[,11]);
y12 ~ ordered_logistic(x12* beta, mu[,12]);
y13 ~ ordered_logistic(x13* beta, mu[,13]);

}

``````

which I split my data to 13 different sub-samples.