# Nested hierarchical model: specifying multiple parameter covariance matrices

I am trying to estimate parameters for an ODE model with a nested hierarchical structure, and I have a question regarding the specification of multiple parameter covariance matrices using Cholesky decomposition.

Because my actual model is very large, please consider the following pseudo-code that specifies parameter covariances at two levels: treatment and subjects (which are nested within treatments).

``````data {
.
.
.
int<lower = 1> n_parameters; // number of parameters (in my_model)
int<lower = 1> n_subjects;
int<lower = 1> n_treatments;
int treatment[n_subjects];
}
parameters {
real par[n_parameters];

vector<lower = 0>[n_parameters] sd_t;
vector<lower = 0>[n_parameters] sd_s;

cholesky_factor_corr[n_parameters] L_t;
cholesky_factor_corr[n_parameters] L_s;

matrix[n_parameters,n_treatments] z_t;
matrix[n_parameters,n_subjects] z_s;
.
.
.
}
transformed parameters{
.
.
.
matrix[n_parameters,n_treatments] t; //treatment level random effects
matrix[n_parameters,n_subjects] s; //subject level random effects
vector[n_parameters] theta[n_subjects];

t = diag_pre_multiply(sd_t, L_t) * z_t;
s = diag_pre_multiply(sd_s, L_s) * z_s;

for(i in 1:n_subjects) {
for(j in 1:n_parameters) {
theta[i][j] = par[j] + t[j,treatment[i]] + s[j,i];
}
}
}

model {
.
.
.
to_vector(z_t) ~ normal(0,1);
to_vector(z_s) ~ normal(0,1);
L_t ~ lkj_corr_cholesky(5);
L_s ~ lkj_corr_cholesky(5);

target += map_rect(my_model, phi, theta, x_r, x_i);
}
``````

I would appreciate if you could comment on whether my pseudo-code is conceptually correct.