Hi,

I’m trying to fit a multivariate related-outcomes multilevel model. One purpose of such a model is to obtain covariances between random effects and covariances between residual errors across my two outcomes (repeated measurements: meas1, meas2) (if interested, see Baldwin et al. (2014)).

The model below works fine as long as my data have the same number of observations per participant for meas1 and meas2. (You can see that I’ve been trying to set it up to allow for unbalanced data though.) All participants have multiple observations for meas1 and meas2.

```
stan_code <- "
data {
int<lower=1> NOBS1; //number of obs for meas 1
int<lower=1> NOBS2; //number of obs for meas 2
int<lower=1> NSUB; //number of subjects
int<lower=1,upper=NSUB> id1[NOBS1]; //id membership for meas 1
int<lower=1,upper=NSUB> id2[NOBS2]; //id membership for meas 2
int<lower=1> meas1_trl[NOBS1]; //trl membership for meas 1
int<lower=1> meas2_trl[NOBS1]; //trl membership for meas 2
vector[NOBS1] meas1; //scores for meas 1
vector[NOBS2] meas2; //scores for meas 2
}
transformed data {
vector[NOBS1] meas1_idv = rep_vector(1,NOBS1); //random effs int meas 1
vector[NOBS2] meas2_idv = rep_vector(1,NOBS2); //random effs int meas 2
vector[2] Y[NOBS1]; //store meas1 and meas2 into same vector for multi_normal_cholesky
//only works b/c NOBS1=NOBS2
for (n in 1:NOBS1) //only works b/c NOBS1=NOBS2
Y[n] = [meas1[n], meas2[n]]';
}
parameters {
vector[2] pop_int; //population-level intercepts
vector<lower=0>[2] sig_err; //population-level residual
cholesky_factor_corr[2] Lrescor; //population-level residuals
cholesky_factor_corr[4] L_omega; //random effects
vector<lower=0>[4] L_sigma; //person-level standard deviations
matrix[4, NSUB] ind_groeffs; //unscaled person-level effects
}
transformed parameters {
matrix[NSUB, 4] rs = (diag_pre_multiply(L_sigma, L_omega) * ind_groeffs)';
//correlations for random effects
vector[2] L_sigma_res = [sig_err[1], sig_err[2]]';
}
model {
//vector[NOBS1] mu1;
//vector[NOBS2] mu2;
vector[2] Mu[NOBS1]; //only possible b/c NOBS1=NOBS2
for (n in 1:NOBS1) {
Mu[n,1] = pop_int[1] +
(rs[,1][id1[n]]) * meas1_idv[n] +
(rs[,2][id1[n]]) * meas1_trl[n];
}
for (n in 1:NOBS2) {
Mu[n,2] = pop_int[2] +
(rs[,3][id2[n]]) * meas2_idv[n] +
(rs[,4][id2[n]]) * meas2_trl[n];
}
//set priors
pop_int ~ normal(0,5);
sig_err ~ cauchy(0,10);
L_sigma ~ cauchy(0,10);
Lrescor ~ lkj_corr_cholesky(2);
L_omega ~ lkj_corr_cholesky(2);
to_vector(ind_groeffs) ~ normal(0,1);
Y ~ multi_normal_cholesky(Mu, diag_pre_multiply(L_sigma_res, Lrescor));
}
generated quantities {
matrix[4,4] omega_ind;
matrix[4,4] sigma_ind;
matrix[2,2] omega_res;
matrix[2,2] sigma_res;
omega_ind = multiply_lower_tri_self_transpose(L_omega);
sigma_ind = quad_form_diag(omega_ind, L_sigma);
omega_res = multiply_lower_tri_self_transpose(Lrescor);
sigma_res = quad_form_diag(omega_res, L_sigma_res);
}
"
```

The problem is that I can’t figure out how to use the multi_normal_cholesky line without making Y contain an equal number of observations for both measurements.

```
Y ~ multi_normal_cholesky(Mu, diag_pre_multiply(L_sigma_res, Lrescor));
```

If I try to estimate them separately the errors separately using something like…

```
meas1 ~ normal(mu1,sig_err[1]);
meas2 ~ normal(mu2,sig_err[2]);
```

… then I’m unable to estimate the covariance between sig_err[1] and sig_err[2], which I need.

Does anyone have any suggestions about the code I need to use in order to get the covariances between residuals once I put in a different number of observations for meas1 and meas2?

Thanks!

Peter