Hierarchical ordered logit model in stan and brms

Modeling
,
1

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,x3means 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

2

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?

3

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.

4

thanks for your reply
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.