I am trying to fit an ordinal model with twist. At this point, I am using only simulated data, so things should be working well. The context is a randomized control trial, with randomization occuring in multiple sites. The twist is that there are three different interventions which are quite similar but not exact. Each site has chooses one intervention. So, I have simulated 9 sites, 3 with intervention A, 3 with intervention B, and 3 with intervention C. I am actually interested in a pooled treatment effect (across all sites and interventions). I model individual site-level effects and intervention-type effects. And since this is ordinal, I need to estimate cut-points - I actually have specified site-specific cut-points rather than using a site-specific random effect, because we donâ€™t want to force a proportional odds assumption across sites, only within a site. I think that should be enough information to understand the model.

My question: I am getting â€śdivergenceâ€ť messages and I want to know if there is a better way for me to parameterize this model to smooth things out a bit. The posterior plots look pretty reasonable, but the trace plots for the treatment effects have some clear spikes. Any suggestions would be much appreciated.

Here is the stan code:

```
data {
int<lower=0> N; // number of observations
int<lower=2> L; // number of ordinal categories
int<lower=1> K; // number of sites
int<lower=1,upper=L> y[N]; // vector of categorical outcomes
int<lower=1,upper=K> kk[N]; // site for individual
int<lower=0,upper=1> ctrl[N]; // treatment or control
int<lower=1,upper=3> cc[K]; // specific control for site
}
parameters {
ordered[L-1] tau[K]; // cut-points for cumulative odds model
vector[K] delta_k; // site specific treatment effect
real<lower=0> eta_0; // sd of delta_k (around delta)
vector[3] delta; // control-specific effect
real<lower=0> eta; // sd of delta (around Delta)
real Delta; // overall control effect
}
transformed parameters{
vector[N] yhat;
for (i in 1:N)
yhat[i] = ctrl[i] * delta_k[kk[i]];
}
model {
// priors
eta_0 ~ cauchy(0, 25);
eta ~ cauchy(0, 25);
for (k in 1:K)
delta_k[k] ~ normal(delta[cc[k]], eta_0);
for (c in 1:3)
delta[c] ~ normal(Delta, eta);
Delta ~ normal(0, 1);
// I also used an induced_dirichlet approach here with
// very similar results:
for (k in 1:K)
for (l in 1:(L-1))
tau[k, l] ~ uniform(-5, 5);
// outcome model
for (i in 1:N)
y[i] ~ ordered_logistic(yhat[i], tau[kk[i]]);
}
```

If it would help for me to submit the simulation code, Iâ€™d be happy to do that - it is a little bit involved, but not too bad.

Here are the trace plots for the treatment effects:

And the trace plots for the variance parameters:

And finally - the density plots for the cut-points (tau in the model) - one panel for each site, and there are 5 levels in the ordinal outcome, hence 4 cut-points: