**EDIT: Just noticed I’m not using the covariances anywhere in the model… still, I haven’t figured out how to add that just yet, any suggestions on how to do that are welcome** (currently trying to understand the multivariate probit regression example in the manual, feel this might be the right way)

I am building a model that will be used to analyze 13 Likert items that supposedly measure the same latent construct.

For now, I just have 2 simulated items. The observations are produced by taking samples from a multivariate normal and then cutting at arbitrary points. I’ve set a high correlation (0.9) between both items, and the simulated data reflects this:

```
1 2 3 4
1 4 3 0 0
2 0 3 0 0
3 0 15 4 3
4 0 1 3 14
```

The model runs without issue* and it recovers most parameters (latent mean, position of cutpoints, variance) just fine, but the covariance remains at zero.

Here’s the model:

```
data {
int<lower=0> N; // # obs
int<lower=2> K; // # item options
int<lower=0> D; // # items
int y[N, D]; // answers
}
parameters {
ordered[K-3] ci[D]; // unknown cutpoints
vector[D] mu; // latent means
cov_matrix[D] Sigma; // latent covariance
}
transformed parameters{
vector[K-1] c[D];
for(d in 1:D){
c[d, 1] = 1.5;
c[d, K-1] = K-0.5;
for(k in 1:(K-3)){
c[d, 1+k] = ci[d, k];
}
}
}
model {
matrix[N, K] p[D]; // item probabilities
//priors
mu ~ normal(rep_vector(0.5*(K+1), D), K*0.5);
for(d in 1:D){
for (k in 1:(K-3)){
ci[d,k] ~ normal(1.5 + k, 1);
}
}
// calculate probabilities
for(d in 1:D){
for(n in 1:N){
p[d,n,1] = Phi((c[d, 1] - mu[d])/diagonal(Sigma)[d]);
p[d,n,K] = 1 - Phi((c[d, K-1] - mu[d])/diagonal(Sigma)[d]);
for(k in 3:K){
p[d,n,k-1] = Phi((c[d, k-1] - mu[d])/diagonal(Sigma)[d]) - Phi((c[d, k-2] - mu[d])/diagonal(Sigma)[d]);
}
y[n,d] ~ categorical(p[d,n]');
}
}
}
```

* I get divergent transitions when the sample size is small and the correlation low, not sure if that’s something to worry about just yet