I constructed a three level multilevel model in stan using the non-centered parameterization approach. In this way my random parameter `beta_0jk[j]`

looks like the following:

```
transformed parameters {
// Varying intercepts
real beta_0jk[Npars];
real beta_0k[Nk];
// Varying intercepts definition
// Level-3 (10 level-3 random intercepts)
for (k in 1:Nk) {
beta_0k[k] = beta_0 + u_0k[k];
}
// Level-2 (100 level-2 random intercepts)
for (j in 1:Npars) {
beta_0jk[j] = beta_0k[countryLookup[j]] + u_0jk[j];
}
}
model {
// Random effects distribution
u_0k ~ normal(0, sigma_u0k);
u_0jk ~ normal(0, sigma_u0jk);
sigma_u0jk ~ student_t(3, 0, 10);
sigma_u0k ~ student_t(3, 0, 10);
}
```

However, I also want to incorporate more random parameters, which I model just like this:

```
transformed parameters {
// Varying intercepts
real beta_0jk[Npars];
real beta_0k[Nk];
real beta_1jk[Npars];
real beta_1k[Nk];
// Varying intercepts definition
// Level-3 (10 level-3 random intercepts)
for (k in 1:Nk) {
beta_0k[k] = beta_0 + u_0k[k];
beta_1k[k] = beta_1 + u_1k[k];
}
// Level-2 (100 level-2 random intercepts)
for (j in 1:Npars) {
beta_0jk[j] = beta_0k[countryLookup[j]] + u_0jk[j];
beta_1jk[j] = beta_1k[countryLookup[j]] + u_1jk[j];
}
}
model {
// Random effects distribution
u_0k ~ normal(0, sigma_u0k);
u_0jk ~ normal(0, sigma_u0jk);
u_1k ~ normal(0, sigma_u1k);
u_1jk ~ normal(0, sigma_u1jk);
sigma_u0jk ~ student_t(3, 0, 10);
sigma_u0k ~ student_t(3, 0, 10);
sigma_u1jk ~ student_t(3, 0, 10);
sigma_u1k ~ student_t(3, 0, 10);
}
```

Here, `beta_1jk`

is my additional random parameter. This works fine and I can estimate my model which gives fine results. However, I doubted whether I did not make a large simplification to my model by not modelling the covariances between `u_0k`

and `u_1k`

and between `u_0jk`

and `u_1jk`

. Does anyone know how I can model this covariance and whether it will make large impact on my model results?

Thanks in advance for the help!