Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.
Checking sampler transitions for divergences.
23 of 4000 (0.57%) transitions ended with a divergence.
These divergent transitions indicate that HMC is not fully able to explore the posterior distribution.
Try increasing adapt delta closer to 1.
If this doesn't remove all divergences, try to reparameterize the model.
Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.
Effective sample size satisfactory.
Split R-hat values satisfactory all parameters.
Processing complete.

// For constrained partial PO model with Dirichlet prior for transformed intercepts
vector pw_log_lik(vector alpha, vector beta, vector tau, vector pposcore,
vector gamma, array[] row_vector X, array[] row_vector Z, array[,] int y,
int k = max(max(y[,1]), max(y[,2])); // assumes all possible categories are observed
if(Nc > 0) eta = eta + gamma[cluster[n]];
a = y[n, 1];    // Y interval censored in [a,b] (uncensored if a=b)
if(a == 1 && b == k) ll = 0.;  // [min,max] interval censoring uninformative
else if (a == 2)   ca  =   alpha[1]     + eta;
else               ca  =   alpha[a - 1] + eta;
if(a > 1 && a < k) ca1 =   alpha[a]     + eta;
if (a == 1)        ca  = -( alpha[1]     + eta + pposcore[2]   * zeta);
else if (a == 2)   ca  =    alpha[1]     + eta + pposcore[2]   * zeta;
else               ca  =    alpha[a - 1] + eta + pposcore[a]   * zeta;
if(a > 1 && a < k) ca1 =    alpha[a]     + eta + pposcore[a+1] * zeta;
if (a == 1 || a == k) ll = log_inv_logit(ca);
else if(b == k) {   // right censored at a
ll = log_inv_logit(alpha[a-1] + eta + r);
else if(a == 1) {   // left censored at b
ll = log_inv_logit(-(alpha[b] + eta + r));
int<lower = 1> N;   // number of observations
int<lower = 1> p;   // number of predictors
int<lower = 0> q;   // number of non-PO predictors in Z
int<lower = 2> k;   // number of outcome categories
int<lower = 0> cn;  // number of contrasts needing priors
int<lower = 0> cn2; // number of non-PO contrasts needing priors
int<lower = 0, upper = k> lpposcore;  // extent of pposcore (1=PO)
matrix[N, p] X;     // matrix of CENTERED predictors
matrix[N, q] Z;     // matrix of CENTERED PPO predictors
array[N, 2] int<lower = 1, upper = k> y; // 2-column outcome on 1 ... k
vector[lpposcore] pposcore; // scores for constrained partial PO
int<lower = 0> Nc;  // number of clusters (0=no clustering)
array[Nc == 0 ? 0 : N] int<lower = 1, upper = Nc> cluster;
int<lower = 1, upper = 2> psigma;  // 1=t(4, rsdmean[1], rsdsd[1]), 2=exponential
array[Nc == 0 ? 0 : 1] real<lower = 0> rsdmean;
array[Nc == 0 || psigma == 2 ? 0 : 1] real<lower = 0> rsdsd;
simplex[k] pi;  // category probabilities for a person w/ average predictors
vector[Nc] gamma_raw;   // unscaled random effects
array[Nc == 0 ? 0 : 1] real<lower = 0> sigmag; // SD of random effects
vector[k - 1] alpha;                               // intercepts
vector[Nc] gamma = Nc == 0 ? gamma_raw : sigmag[1] * gamma_raw;
vector[N] log_lik;                                 // log-likelihood pieces
for (j in 2:k) alpha[j - 1] = logit(sum(pi[j:k])); // predictors are CENTERED
log_lik = pw_log_lik(alpha, beta, tau, pposcore, gamma, Xr, Zr, y, cluster);
gamma_raw ~ std_normal();  // implies gamma ~ normal(0, sigmag)
if(psigma == 1) sigmag ~ student_t(4, rsdmean[1], rsdsd[1]);
else sigmag ~ exponential(1. / rsdmean[1]);
target += dirichlet_lpdf(pi | rep_vector(conc, k));
if(cn > 0)  target += normal_lpdf(C * beta | cmus, csds);
if(cn2 > 0) target += normal_lpdf(C2 * tau | cmus2, csds2);