Hiya folks. I’ve been trying to better understand the use of the MVN distribution for priors in a hierarchical model. Ultimately, I’m hoping to develop models with MVN priors on regression coefficients and covariance priors informed by group similarity. For now, however, I am starting with trying whats in the SUG here. I first tried to modify the example code to get rid of the group-level predictors because I don’t often deal with those and wanted to keep it simple. I simulated some code form that model and tried running it, only to face fine convergence on reasonable estimates but many divergent transitions. Then, I decided to try to recreate exactly whats in the SUG (including the group-level predictors), but similarly ran into many (180ish) divergent transitions. Here is the Stan code (copied from the above link) that I am working with:

```
data {
int<lower=0> N; // num individuals
int<lower=1> K; // num ind predictors
int<lower=1> J; // num groups
int<lower=1> L; // num group predictors
array[N] int<lower=1, upper=J> jj; // group for individual
matrix[N, K] x; // individual predictors
array[J] row_vector[L] u; // group predictors
vector[N] y; // outcomes
}
parameters {
corr_matrix[K] Omega; // prior correlation
vector<lower=0>[K] tau; // prior scale
matrix[L, K] gamma; // group coeffs
array[J] vector[K] beta; // indiv coeffs by group
real<lower=0> sigma; // prediction error scale
}
model {
tau ~ cauchy(0, 2.5);
Omega ~ lkj_corr(2);
to_vector(gamma) ~ normal(0, 5);
{
array[J] row_vector[K] u_gamma;
for (j in 1:J) {
u_gamma[j] = u[j] * gamma;
}
beta ~ multi_normal(u_gamma, quad_form_diag(Omega, tau));
}
for (n in 1:N) {
y[n] ~ normal(x[n] * beta[jj[n]], sigma);
}
}
```

Here is the code that I wrote to try to simulate an appropriate dataset in R:

```
set.seed(6)
# data
N = 600
K = 3
J = 3
L = 3
jj = rep(1:J, each = N/J)
x = matrix(data = rnorm(N*K), nrow = N, ncol = K)
u = matrix(data = rnorm(J*L), nrow = J, ncol = L)
# parameters
Omega <- rlkjcorr(1, K, 2)
tau <- rhalfcauchy(K, 0.25)
Sigma <- diag(tau) %*% Omega %*% diag(tau)
gamma <- matrix(data = rnorm(L*K, 0, 5), nrow = L, ncol = K)
sigma <- 1
beta <- matrix(NA, nrow = J, ncol = K)
for (j in 1:J) {
u_gamma <- u[j, ] %*% gamma
beta[j, ] <- mvrnorm(1, mu = u_gamma, Sigma = Sigma)
}
y <- numeric(N)
for (n in 1:N) {
y[n] <- rnorm(1, mean = sum(x[n, ] * beta[jj[n], ]), sd = sigma)
}
stan_data <- list(
N = N,
K = K,
J = J,
L = L,
jj = jj,
x = x,
u = u,
y = y
)
test <- stan(file = "~/Downloads/MVN_from_guide.stan",
data = stan_data,
chains = 2, cores = 2)
```

I suspect that the issue has to do with my simulation as I am pretty new to working with MVN distributions. Any help in identifying why I might be struggling to fit the model in the SUG would be greatly appreciated!