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:
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]
compiled_model <- stan_model(model_code =
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?