I’d like to fit a hierarchical model in stan, where the response is *not* univariate. However, I can only find examples

for when the response is univariate. Here’s a super simple example model for the univariate case, just to be clear about what I mean:

```
set.seed(42)
obs_per_group <- 50
data <- list(num_groups = 10)
data$num_obs <- data$num_groups*obs_per_group
data$Y <- rnorm(data$num_obs)
data$group_inds <- rep(1:data$num_groups,each=obs_per_group)
data$Y <- data$Y + rnorm(data$num_groups)[data$group_inds]
library('rstan')
compiled_model <- stan_model(model_code =
"
data {
int<lower=1> num_obs;
int<lower=1> num_groups;
vector[num_obs] Y;
int<lower=1> group_inds[num_obs];
}
parameters {
real mu;
real sigma_offset;
real sigma_obs;
vector[num_groups] offset_per_group;
}
model {
vector[num_obs] yhat;
mu ~ normal(0,1);
sigma_offset ~ lognormal(0,1);
offset_per_group ~ normal(0, sigma_offset);
sigma_obs ~ lognormal(0,1);
yhat = mu + offset_per_group[group_inds];
Y ~ normal(yhat, sigma_obs);
}
")
```

Now, I’d like to fit a similar model, but where Y is a matrix with K columns, one for each response-per-row.

The tricky thing is that, not only do I want to model the population-level covariance/correlation between the responses-- I also want to capture any per-group differences from this average amount of covariance/correlation. I.e. this group is especially correlated between his Y[,1] and Y[,2] responses, this other group is not very correlated between his Y[,1] and Y[,2] responses, etc.

So in the univariate case, including per-group offsets is as easy as add mu + offset_per_group. But I’m not sure how this would work in the case of correlation or covariance matrices? How do I generate the per-group covariance/correlation matrices?