I’m trying to speed up the model below. I’ve already tried QR reparamaterisation of the data and noncentered parameterisation of the level 2 effects to no avail. I even tried a map_rect approach that calculates the right answer but made things amazingly slower over all.

Can anyone suggest any other tips or tricks to speed things up? I was particularly wondering if there is any matrix magic that could improve on the nested loops? Alternatively what would be a smart way to implememnt map_rect in this model?

Thanks in advance.

```
data{
int<lower=0> NvarsX; // num independent variables
int<lower=0> NvarsY; // num dependent variables
int<lower=0> N; // Total number of rows
int<lower=0> N_ind; // Number of individuals/participants
int<lower=1> indiv[N]; // data to identify individuals
matrix[N, NvarsX] x; // data for independent vars
vector[NvarsY] y [N]; // data for dependent vars
}
parameters{
cholesky_factor_corr[NvarsY] L_Omega; //For the correlations between outcome variables
vector<lower=0>[NvarsY] L_sigma; //Residual SD of outcome variables
matrix[NvarsY, N_ind] Zbeta0_ind; //Intercepts at individual level
matrix[NvarsX, N_ind] Zbeta_ind [NvarsY]; //Coefficients at the individual level
vector<lower=0> [NvarsY] Tau0; //Random effect for individual intercepts
matrix<lower=0> [NvarsY, NvarsX] Tau; //Random effect for indvidiual coefficients
vector[NvarsY] Beta0; // Level 2 intercepts for each Y
vector[NvarsX] Beta [NvarsY]; // Level 2 coefficients for each X vs each Y
}
transformed parameters{
vector[NvarsY] mu[N];
matrix[NvarsY, NvarsY] L_Sigma;
matrix[NvarsY, N_ind] beta0_ind;
matrix[NvarsX, N_ind] beta_ind [NvarsY];
L_Sigma = diag_pre_multiply(L_sigma, L_Omega);
// Define Individual level effects - non parametric specification
for (s in 1:N_ind){
for (dv in 1:NvarsY){
beta0_ind[dv, s] = Zbeta0_ind[dv, s] * Tau0[dv] + Beta0[dv];
for (iv in 1:NvarsX){
beta_ind[dv, iv, s] = Zbeta_ind[dv, iv, s] * Tau[dv, iv] + Beta[dv, iv];
}
}
}
// Define level 2 effects
for (i in 1:N){
for (dv in 1:NvarsY){
mu[i, dv] = beta0_ind[dv, indiv[i]] +
dot_product( to_vector( beta_ind[dv, 1:NvarsX, indiv[i]] ), x[i, 1:NvarsX] );
}
}
}
model{
// Priors
L_Omega ~ lkj_corr_cholesky(1);
L_sigma ~ normal(0, 1);
Tau0 ~ normal(0,1);
to_vector(Tau) ~ normal(0,1);
for (s in 1:N_ind){
for (dv in 1:NvarsY){
Zbeta0_ind[dv, s] ~ normal(0,1);
Zbeta_ind[dv, 1:NvarsX, s] ~ normal(0,1);
}
}
to_vector(Beta0) ~ cauchy(0, 5);
//for( dv in 1:NvarsY){
// to_vector(Beta[dv, 1:NvarsX]) ~ normal(0, 2);
//}
// Likelihood
for (i in 1:N){
y[i, 1:NvarsY] ~ multi_normal_cholesky(mu[i, 1:NvarsY], L_Sigma);
}
}
generated quantities{
matrix[NvarsY, NvarsY] Omega;
matrix[NvarsY, NvarsY] Sigma;
Omega = multiply_lower_tri_self_transpose(L_Omega);
Sigma = quad_form_diag(L_Omega, L_sigma);
}
```