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.
Thank you in advance for your input.