Treedepth satisfactory for all transitions.
Checking sampler transitions for divergences.
Checking E-BFMI - sampler transitions HMC potential energy.
Split R-hat values satisfactory all parameters.
Processing complete, no problems detected.

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