Hello, I am a beginner of STAN and am reading User’s Guide and Reference Manual now.

I would like to ask any opinion or thoughts to extend the multivariate outcomes model at Chapter 9 to multivariate ‘multilevel’ model.

For example, I had a data set which had two scores (math and science) at level-1 nested students at level-2, middle school at level-3, and high school at level-4.

Because I will not add error term at level-1 because the level-1 error is omitted from the model because level one is only used to partition the multivariate outcomes’ structure into the model (Hox, 2010; Goldstein, 2010).

Therefore, it will be actually 3-level multivariate model as shown below.

Y = (Gamma1 + u11 + u21 + r1)*dummy1 + (Gamma2 + u12 + u22 + r2)*dummy2

(Actually, the middle schools are not purely nested to high school thus it is cross-classified, but as I search in the forum, it would not require specific codes to set the characteristic. If I am wrong, please let me know.)

However, I am struggle with extending the multivariate model to multivariate multilevel model because I cannot find the example to make 3-level multilevel model with multiple responses…

Could you give any your thoughts or suggestion about it? Or, is there any study or example I can follow to understand how the code can be set up to multivariate multilevel model?

Here is the code what I saw at the Manual (Chapter 9.15 Multivariate outcomes)

data {

int<lower=1> K;

int<lower=1> J;

int<lower=0> N;

vector[J] x[N];

vector[K] y[N];

}

parameters {

matrix[K, J] beta;

cholesky_factor_corr[K] L_Omega;

vector<lower=0>[K] L_sigma;

}

model {

vector[K] mu[N];

matrix[K, K] L_Sigma;

for (n in 1:N)

mu[n] = beta * x[n];

L_Sigma = diag_pre_multiply(L_sigma, L_Omega);

to_vector(beta) ~ normal(0, 5);

L_Omega ~ lkj_corr_cholesky(4);

L_sigma ~ cauchy(0, 2.5);

y ~ multi_normal_cholesky(mu, L_Sigma);

}

Thank you for your time in advance.