Nested hierarchical model: specifying multiple parameter covariance matrices


#1

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.