Hi everyone, I am using a hierarchical model for a meta analysis, with a model very similar to the original Rubin’s Randome Effect model. Some observations that we have in our dataset come from the same RCT, either different treatment arms or different timings of measurement and so I would like to put a correlation structure, while maintaining independece across different RCTs.

I have an issue in the way I declare the block covariance matrix, because Stan does not allow dynamic arrays and groups are of different sizes. I tried some work-arounds by declaring matrices of excessive sizes and then only using the dimension I actually need, but the model does not perform well and I suspect it is because of my cave man programming and approach.

Does anyone have a good fix for it? I am sure the problem is not new, having within group block covariance matrices.

```
data {
int<lower=0> N; // num observations
int<lower=0> J; // num RCTs (groups)
int<lower=1, upper=J> group[N]; // group indicator (which RCT each observation belongs to)
vector[N] y; // outcome
vector[N] se; // vector of standard error estimates
int<lower=0> group_size[J]; // number of observations per group (size of each RCT)
}
parameters {
real t; // common mean for the theta's
real<lower=0> sigma_theta; // standard deviation of group means
// Latent variables for within-group correlation for groups with multiple observations
vector[N] z; // uncorrelated latent variables for observations
cholesky_factor_corr[max(group_size)] L_Omega[J]; // Cholesky factor of correlation matrix for each group
}
transformed parameters {
vector[N] theta; // Observation-level theta with group correlations
// Apply index trick to theta[j] directly for groups with multiple vs single observations
for (n in 1:N) {
theta[n] = t + z[n] * sigma_theta;
}
}
model {
// Priors
t ~ normal(0, 5); // Prior for the common mean
sigma_theta ~ normal(0, 5); // Prior for the group-level standard deviation
for (j in 1:J) {
L_Omega[j] ~ lkj_corr_cholesky(1); // Prior for the Cholesky factor of each group
}
// Correlation for groups with multiple observations
for (j in 1:J) {
if (group_size[j] > 1) {
vector[group_size[j]] z_group;
int start_idx = sum(head(group_size, j - 1)) + 1;
z_group = z[start_idx:(start_idx + group_size[j] - 1)];
// Extract submatrix from L_Omega to match group size
matrix[group_size[j], group_size[j]] L_sub = L_Omega[j][1:group_size[j], 1:group_size[j]];
z_group ~ multi_normal_cholesky(rep_vector(0, group_size[j]), L_sub); // Correlated deviations within groups
} else {
z[j] ~ normal(0, 1); // Single observation groups get standard normal deviations
}
}
// Likelihood for observations
for (n in 1:N) {
y[n] ~ normal(theta[n], se[n]);
}
}
```

Any suggestions? (secretly hoping @betanalpha has an answer)