Brms sampling in multi-level multinomial models

Hi all,

I am running a multi-level multinomial model, but having some issues fitting it efficiently with brms. First off, I am aware fitting Bayesian models can be computationally heavy and it requires time, but my models run more than 30 days without finishing and I am afraid I have overlooked something important.

The model
I want to model the mapping between seven decorrelated continuous variables and 7 categories. Moreover, I want to describe how this mapping differs between demographic variables (like sex, country or language). I am using the following intercept-free model:

model <- brm(
  category ~ 0 + RC1 + RC2 + RC3 + RC4 + RC5 + RC6 + RC7 + (0 + RC1 + RC2 + RC3 + RC4 + RC5 + RC6 + RC7 | sex + country:language + speaker),
  data = dat_list,
  family = categorical(link = "logit"),
  prior = set_prior("normal(0,1)", class = "b"),
  control = list(adapt_delta = 0.999, max_treedepth = 15))
)

The full dataset consists of 52k samples. Most group-level effects contain a few levels, but one factor contains 432 levels (speakers).

Observations
Somehow the model does not scale well:

  • The full model as specified above runs in reasonable time and without warnings on a subset of the data
  • An intercept-only model on all data also runs in reasonable time and without warnings (RC1 + RC2 + RC3 + RC4 + RC5 + RC6 + RC7 + (1 | sex) + (1 | country:language) + (1 | speaker))
  • A population-level model only runs through very fast on the full dataset (RC1 + RC2 + RC3 + RC4 + RC5 + RC6 + RC7)
  • But sampling from the full dataset with the full model does not work well (often it does not even reach the 10% marker)

I played around with a divide and conquer strategy (only taking into account the population effects, came across the idea here), splitting the data up in equal bins and running smaller models on the subgroups and later average over the posteriors. The result was that the estimates were often off and the estimated distribution was too wide. Also, I am not sure if such strategy is advisable given we use a hierarchical model, pooling across groups in all data and not just subsets, so I rejected this approach.

A simple way to speedup the sampling is to simply run on more CPUs at once (moved from 4 to 32 CPUs), use one chain per CPU and obtain less posterior samples per chain.

Another way to increase sampling speed, is to experiment with the number of warmup iterations, to keep this phase as short as possible while having reliable estimates. Although I noted that when using this many chains the chains took a different time to complete, which shouldn’t really happen, right?

Finally, using less variables would probably speedup the model too, but this conflicts with my scientific interest.

Work in progress
Currently, I am looking into improving the grain size and investigating if it is possible to profile the execution of the model, to investigate what is going wrong.

Open questions
It might be wise to specify better priors, but I have a hard time coming up with a good prior. I picked a standard normal prior, because the variables are centered around 0 and roughly fall into the range of a standard normal across all categories. I would be interested if you could suggest a prior with better sampling properties.

Finally, I would be interested if you have some suggestions how to improve sampling speed and if I overlooked some possibilities.

Thank you for taking the time to read this topic!

Hi!

Are you able to share the Stan code generated by brms? You may be able to get some help optimizing the code, which will in turn speed up sampling.

I received some fantastic help in this thread that took days off of sampling time.

2 Likes

@polvanrijn follow the suggestions in the thread that @JLC linked, and also:

  1. Vectorize every for loop in your code
  2. You can apply within-chain parallelization in brms. See this vignette
  3. Tell brms to perform a QR decomposition on the data (this helps to speed things up in Hierarchical Models): bf(formula, family = XX, decomp = "QR")
  4. Try using CdmStanR instead of RStan (for me sometimes I see a ~3x speedup): brm(formula, data, backend = "cmdstanr").
  5. You can tell brms to use log unormalized PDF (*_lupdf, *_lpdm) instead of the log normalized PDF (*_lpdf) and *_lpmf). This is available with the latest cmdstanr and the GitHub version of brms: brm(..., backend = "cmdstanr", normalize = FALSE). See this(issue and this pull request.

Those 5 things may speedup things. Also, if you can, please provide us with the code that brms generated, so people can help to do more speedup.

3 Likes

Dear all,
Thanks for your quick replies. @JLC, thanks for sharing this amazing thread and @storopoli thanks for showing many possible improvements :-) Sorry for not getting back to you earlier, but I wanted to try your suggestions before replying. I will go through the generated Stan code block-wise and indicate what I changed/want to change. The code was generated with:

make_stancode(
  emotion ~ 0 + RC1 + RC2 + RC3 + RC4 + RC5 + RC6 + RC7 + (0 + RC1 + RC2 + RC3 + RC4 + RC5 + RC6 + RC7 | sex + country:language + speaker),
  data = dat,
  family = categorical(link = "logit"),
  prior = set_prior("normal(0,1)", class = "b")
)

Functions

// generated with brms 2.14.0
functions {
  /* turn a vector into a matrix of defined dimension
   * stores elements in row major order
   * Args:
   *   X: a vector
   *   N: first dimension of the desired matrix
   *   K: second dimension of the desired matrix
   * Returns:
   *   a matrix of dimension N x K
   */
  matrix as_matrix(vector X, int N, int K) {
    matrix[N, K] Y;
    for (i in 1:N) {
      Y[i] = to_row_vector(X[((i - 1) * K + 1):(i * K)]);
    }
    return Y;
  }
 /* compute correlated group-level effects
  * Args:
  *   z: matrix of unscaled group-level effects
  *   SD: vector of standard deviation parameters
  *   L: cholesky factor correlation matrix
  * Returns:
  *   matrix of scaled group-level effects
  */
  matrix scale_r_cor(matrix z, vector SD, matrix L) {
    // r is stored in another dimension order than z
    return transpose(diag_pre_multiply(SD, L) * z);
  }
}

It surprised me that we need scale_r_cor whereas this does not occur in @JLC 's thread. @JLC and @storopoli, do you know why correlated group-level effects are computed in this model, but not in the model in the linked thread? The reason is that in this model we compute group-level coefficients and not only group-level intercepts (when replacing the group-level coefficients with an intercept in make_stancode, the correlation matrixes disappear).

Data

data {
  int<lower=1> N;  // total number of observations
  int<lower=2> ncat;  // number of categories
  int Y[N];  // response variable
  int<lower=1> K_muANG;  // number of population-level effects
  matrix[N, K_muANG] X_muANG;  // population-level design matrix
  int<lower=1> K_muDIS;  // number of population-level effects
  matrix[N, K_muDIS] X_muDIS;  // population-level design matrix
  int<lower=1> K_muFER;  // number of population-level effects
  matrix[N, K_muFER] X_muFER;  // population-level design matrix
  int<lower=1> K_muHAP;  // number of population-level effects
  matrix[N, K_muHAP] X_muHAP;  // population-level design matrix
  int<lower=1> K_muSAD;  // number of population-level effects
  matrix[N, K_muSAD] X_muSAD;  // population-level design matrix
  int<lower=1> K_muSUR;  // number of population-level effects
  matrix[N, K_muSUR] X_muSUR;  // population-level design matrix
  // data for group-level effects of ID 1
  int<lower=1> N_1;  // number of grouping levels
  int<lower=1> M_1;  // number of coefficients per level
  int<lower=1> J_1[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_1_muANG_1;
  vector[N] Z_1_muANG_2;
  vector[N] Z_1_muANG_3;
  vector[N] Z_1_muANG_4;
  vector[N] Z_1_muANG_5;
  vector[N] Z_1_muANG_6;
  vector[N] Z_1_muANG_7;
  int<lower=1> NC_1;  // number of group-level correlations
  // data for group-level effects of ID 2
  int<lower=1> N_2;  // number of grouping levels
  int<lower=1> M_2;  // number of coefficients per level
  int<lower=1> J_2[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_2_muANG_1;
  vector[N] Z_2_muANG_2;
  vector[N] Z_2_muANG_3;
  vector[N] Z_2_muANG_4;
  vector[N] Z_2_muANG_5;
  vector[N] Z_2_muANG_6;
  vector[N] Z_2_muANG_7;
  int<lower=1> NC_2;  // number of group-level correlations
  // data for group-level effects of ID 3
  int<lower=1> N_3;  // number of grouping levels
  int<lower=1> M_3;  // number of coefficients per level
  int<lower=1> J_3[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_3_muANG_1;
  vector[N] Z_3_muANG_2;
  vector[N] Z_3_muANG_3;
  vector[N] Z_3_muANG_4;
  vector[N] Z_3_muANG_5;
  vector[N] Z_3_muANG_6;
  vector[N] Z_3_muANG_7;
  int<lower=1> NC_3;  // number of group-level correlations
  // data for group-level effects of ID 4
  int<lower=1> N_4;  // number of grouping levels
  int<lower=1> M_4;  // number of coefficients per level
  int<lower=1> J_4[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_4_muDIS_1;
  vector[N] Z_4_muDIS_2;
  vector[N] Z_4_muDIS_3;
  vector[N] Z_4_muDIS_4;
  vector[N] Z_4_muDIS_5;
  vector[N] Z_4_muDIS_6;
  vector[N] Z_4_muDIS_7;
  int<lower=1> NC_4;  // number of group-level correlations
  // data for group-level effects of ID 5
  int<lower=1> N_5;  // number of grouping levels
  int<lower=1> M_5;  // number of coefficients per level
  int<lower=1> J_5[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_5_muDIS_1;
  vector[N] Z_5_muDIS_2;
  vector[N] Z_5_muDIS_3;
  vector[N] Z_5_muDIS_4;
  vector[N] Z_5_muDIS_5;
  vector[N] Z_5_muDIS_6;
  vector[N] Z_5_muDIS_7;
  int<lower=1> NC_5;  // number of group-level correlations
  // data for group-level effects of ID 6
  int<lower=1> N_6;  // number of grouping levels
  int<lower=1> M_6;  // number of coefficients per level
  int<lower=1> J_6[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_6_muDIS_1;
  vector[N] Z_6_muDIS_2;
  vector[N] Z_6_muDIS_3;
  vector[N] Z_6_muDIS_4;
  vector[N] Z_6_muDIS_5;
  vector[N] Z_6_muDIS_6;
  vector[N] Z_6_muDIS_7;
  int<lower=1> NC_6;  // number of group-level correlations
  // data for group-level effects of ID 7
  int<lower=1> N_7;  // number of grouping levels
  int<lower=1> M_7;  // number of coefficients per level
  int<lower=1> J_7[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_7_muFER_1;
  vector[N] Z_7_muFER_2;
  vector[N] Z_7_muFER_3;
  vector[N] Z_7_muFER_4;
  vector[N] Z_7_muFER_5;
  vector[N] Z_7_muFER_6;
  vector[N] Z_7_muFER_7;
  int<lower=1> NC_7;  // number of group-level correlations
  // data for group-level effects of ID 8
  int<lower=1> N_8;  // number of grouping levels
  int<lower=1> M_8;  // number of coefficients per level
  int<lower=1> J_8[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_8_muFER_1;
  vector[N] Z_8_muFER_2;
  vector[N] Z_8_muFER_3;
  vector[N] Z_8_muFER_4;
  vector[N] Z_8_muFER_5;
  vector[N] Z_8_muFER_6;
  vector[N] Z_8_muFER_7;
  int<lower=1> NC_8;  // number of group-level correlations
  // data for group-level effects of ID 9
  int<lower=1> N_9;  // number of grouping levels
  int<lower=1> M_9;  // number of coefficients per level
  int<lower=1> J_9[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_9_muFER_1;
  vector[N] Z_9_muFER_2;
  vector[N] Z_9_muFER_3;
  vector[N] Z_9_muFER_4;
  vector[N] Z_9_muFER_5;
  vector[N] Z_9_muFER_6;
  vector[N] Z_9_muFER_7;
  int<lower=1> NC_9;  // number of group-level correlations
  // data for group-level effects of ID 10
  int<lower=1> N_10;  // number of grouping levels
  int<lower=1> M_10;  // number of coefficients per level
  int<lower=1> J_10[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_10_muHAP_1;
  vector[N] Z_10_muHAP_2;
  vector[N] Z_10_muHAP_3;
  vector[N] Z_10_muHAP_4;
  vector[N] Z_10_muHAP_5;
  vector[N] Z_10_muHAP_6;
  vector[N] Z_10_muHAP_7;
  int<lower=1> NC_10;  // number of group-level correlations
  // data for group-level effects of ID 11
  int<lower=1> N_11;  // number of grouping levels
  int<lower=1> M_11;  // number of coefficients per level
  int<lower=1> J_11[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_11_muHAP_1;
  vector[N] Z_11_muHAP_2;
  vector[N] Z_11_muHAP_3;
  vector[N] Z_11_muHAP_4;
  vector[N] Z_11_muHAP_5;
  vector[N] Z_11_muHAP_6;
  vector[N] Z_11_muHAP_7;
  int<lower=1> NC_11;  // number of group-level correlations
  // data for group-level effects of ID 12
  int<lower=1> N_12;  // number of grouping levels
  int<lower=1> M_12;  // number of coefficients per level
  int<lower=1> J_12[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_12_muHAP_1;
  vector[N] Z_12_muHAP_2;
  vector[N] Z_12_muHAP_3;
  vector[N] Z_12_muHAP_4;
  vector[N] Z_12_muHAP_5;
  vector[N] Z_12_muHAP_6;
  vector[N] Z_12_muHAP_7;
  int<lower=1> NC_12;  // number of group-level correlations
  // data for group-level effects of ID 13
  int<lower=1> N_13;  // number of grouping levels
  int<lower=1> M_13;  // number of coefficients per level
  int<lower=1> J_13[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_13_muSAD_1;
  vector[N] Z_13_muSAD_2;
  vector[N] Z_13_muSAD_3;
  vector[N] Z_13_muSAD_4;
  vector[N] Z_13_muSAD_5;
  vector[N] Z_13_muSAD_6;
  vector[N] Z_13_muSAD_7;
  int<lower=1> NC_13;  // number of group-level correlations
  // data for group-level effects of ID 14
  int<lower=1> N_14;  // number of grouping levels
  int<lower=1> M_14;  // number of coefficients per level
  int<lower=1> J_14[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_14_muSAD_1;
  vector[N] Z_14_muSAD_2;
  vector[N] Z_14_muSAD_3;
  vector[N] Z_14_muSAD_4;
  vector[N] Z_14_muSAD_5;
  vector[N] Z_14_muSAD_6;
  vector[N] Z_14_muSAD_7;
  int<lower=1> NC_14;  // number of group-level correlations
  // data for group-level effects of ID 15
  int<lower=1> N_15;  // number of grouping levels
  int<lower=1> M_15;  // number of coefficients per level
  int<lower=1> J_15[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_15_muSAD_1;
  vector[N] Z_15_muSAD_2;
  vector[N] Z_15_muSAD_3;
  vector[N] Z_15_muSAD_4;
  vector[N] Z_15_muSAD_5;
  vector[N] Z_15_muSAD_6;
  vector[N] Z_15_muSAD_7;
  int<lower=1> NC_15;  // number of group-level correlations
  // data for group-level effects of ID 16
  int<lower=1> N_16;  // number of grouping levels
  int<lower=1> M_16;  // number of coefficients per level
  int<lower=1> J_16[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_16_muSUR_1;
  vector[N] Z_16_muSUR_2;
  vector[N] Z_16_muSUR_3;
  vector[N] Z_16_muSUR_4;
  vector[N] Z_16_muSUR_5;
  vector[N] Z_16_muSUR_6;
  vector[N] Z_16_muSUR_7;
  int<lower=1> NC_16;  // number of group-level correlations
  // data for group-level effects of ID 17
  int<lower=1> N_17;  // number of grouping levels
  int<lower=1> M_17;  // number of coefficients per level
  int<lower=1> J_17[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_17_muSUR_1;
  vector[N] Z_17_muSUR_2;
  vector[N] Z_17_muSUR_3;
  vector[N] Z_17_muSUR_4;
  vector[N] Z_17_muSUR_5;
  vector[N] Z_17_muSUR_6;
  vector[N] Z_17_muSUR_7;
  int<lower=1> NC_17;  // number of group-level correlations
  // data for group-level effects of ID 18
  int<lower=1> N_18;  // number of grouping levels
  int<lower=1> M_18;  // number of coefficients per level
  int<lower=1> J_18[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_18_muSUR_1;
  vector[N] Z_18_muSUR_2;
  vector[N] Z_18_muSUR_3;
  vector[N] Z_18_muSUR_4;
  vector[N] Z_18_muSUR_5;
  vector[N] Z_18_muSUR_6;
  vector[N] Z_18_muSUR_7;
  int<lower=1> NC_18;  // number of group-level correlations
  int prior_only;  // should the likelihood be ignored?

Simplifying X

Replace

  int<lower=1> K_muANG;  // number of population-level effects
  matrix[N, K_muANG] X_muANG;  // population-level design matrix
  int<lower=1> K_muDIS;  // number of population-level effects
  matrix[N, K_muDIS] X_muDIS;  // population-level design matrix
  int<lower=1> K_muFER;  // number of population-level effects
  matrix[N, K_muFER] X_muFER;  // population-level design matrix
  int<lower=1> K_muHAP;  // number of population-level effects
  matrix[N, K_muHAP] X_muHAP;  // population-level design matrix
  int<lower=1> K_muSAD;  // number of population-level effects
  matrix[N, K_muSAD] X_muSAD;  // population-level design matrix
  int<lower=1> K_muSUR;  // number of population-level effects
  matrix[N, K_muSUR] X_muSUR;  // population-level design matrix

with

  int<lower=1> K;  // number of population-level effects
  matrix[N, K] X;  // population-level design matrix

Since K will always be 7 (seven factors) and the design matrix (X) will be the same for all emotions.

Simplifying group-level effects

As described in the thread, we only need one Z per group-level effect and not for every emotion separately

  // data for group-level effects of ID 1
  int<lower=1> N_1;  // number of grouping levels
  int<lower=1> M_1;  // number of coefficients per level
  int<lower=1> J_1[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_1_muANG_1;
  vector[N] Z_1_muANG_2;
  vector[N] Z_1_muANG_3;
  vector[N] Z_1_muANG_4;
  vector[N] Z_1_muANG_5;
  vector[N] Z_1_muANG_6;
  vector[N] Z_1_muANG_7;

ID 1 refers to the interaction between countries and languages (for simplicity, I’ll call it ‘culture’ here), ID 2 sex and ID 3 single speakers.

Since the group-level effects are computed for the same coefficients as the population effects, we can omit M_1, M_2 and M_3 since it will always be 7 (number of factors). So this part can be rewritten to:

  // data for group-level effects of culture
  int<lower=1> N_culture;  // number of cultures
  int<lower=1> J_culture[N];  // culture grouping indicator per observation
  // culture group-level predictor values
  vector[N] Z_culture_RC1;
  vector[N] Z_culture_RC2;
  vector[N] Z_culture_RC3;
  vector[N] Z_culture_RC4;
  vector[N] Z_culture_RC5;
  vector[N] Z_culture_RC6;
  vector[N] Z_culture_RC7;
  
  // data for group-level effects of sex
  int<lower=1> N_sex;  // number of sexes
  int<lower=1> J_sex[N];  // sex grouping indicator per observation
  // sex group-level predictor values
  vector[N] Z_sex_RC1;
  vector[N] Z_sex_RC2;
  vector[N] Z_sex_RC3;
  vector[N] Z_sex_RC4;
  vector[N] Z_sex_RC5;
  vector[N] Z_sex_RC6;
  vector[N] Z_sex_RC7;

  // data for group-level effects of speaker
  int<lower=1> N_speaker;  // number of speakers
  int<lower=1> J_speaker[N];  // speaker grouping indicator per observation
  // speaker group-level predictor values
  vector[N] Z_speaker_RC1;
  vector[N] Z_speaker_RC2;
  vector[N] Z_speaker_RC3;
  vector[N] Z_speaker_RC4;
  vector[N] Z_speaker_RC5;
  vector[N] Z_speaker_RC6;
  vector[N] Z_speaker_RC7;

@JLC and @storopoli, why do we need to define Z anyway? Isn’t this information already contained in X?

Simplifying NC

NC_1 to NC_18 is always 21 because this is the number of items in the upper diagonal of a correlation matrix of the size 7x7, so one parameter NC would suffice.

Part 2:

Transformed data

transformed data {
}

Is empty by default, but when applying QR decomposition it will be

transformed data {
  // matrices for QR decomposition
  matrix[N, K_muANG] XQ_muANG;
  matrix[K_muANG, K_muANG] XR_muANG;
  matrix[K_muANG, K_muANG] XR_muANG_inv;
  
  // matrices for QR decomposition
  matrix[N, K_muDIS] XQ_muDIS;
  matrix[K_muDIS, K_muDIS] XR_muDIS;
  matrix[K_muDIS, K_muDIS] XR_muDIS_inv;
  
  // matrices for QR decomposition
  matrix[N, K_muFER] XQ_muFER;
  matrix[K_muFER, K_muFER] XR_muFER;
  matrix[K_muFER, K_muFER] XR_muFER_inv;
  
  // matrices for QR decomposition
  matrix[N, K_muHAP] XQ_muHAP;
  matrix[K_muHAP, K_muHAP] XR_muHAP;
  matrix[K_muHAP, K_muHAP] XR_muHAP_inv;
  
  // matrices for QR decomposition
  matrix[N, K_muSAD] XQ_muSAD;
  matrix[K_muSAD, K_muSAD] XR_muSAD;
  matrix[K_muSAD, K_muSAD] XR_muSAD_inv;
  
  // matrices for QR decomposition
  matrix[N, K_muSUR] XQ_muSUR;
  matrix[K_muSUR, K_muSUR] XR_muSUR;
  matrix[K_muSUR, K_muSUR] XR_muSUR_inv;
  
  // compute and scale QR decomposition
  XQ_muANG = qr_thin_Q(X_muANG) * sqrt(N - 1);
  XR_muANG = qr_thin_R(X_muANG) / sqrt(N - 1);
  XR_muANG_inv = inverse(XR_muANG);
  
  // compute and scale QR decomposition
  XQ_muDIS = qr_thin_Q(X_muDIS) * sqrt(N - 1);
  XR_muDIS = qr_thin_R(X_muDIS) / sqrt(N - 1);
  XR_muDIS_inv = inverse(XR_muDIS);
  
  // compute and scale QR decomposition
  XQ_muFER = qr_thin_Q(X_muFER) * sqrt(N - 1);
  XR_muFER = qr_thin_R(X_muFER) / sqrt(N - 1);
  XR_muFER_inv = inverse(XR_muFER);
  
  // compute and scale QR decomposition
  XQ_muHAP = qr_thin_Q(X_muHAP) * sqrt(N - 1);
  XR_muHAP = qr_thin_R(X_muHAP) / sqrt(N - 1);
  XR_muHAP_inv = inverse(XR_muHAP);
  
  // compute and scale QR decomposition
  XQ_muSAD = qr_thin_Q(X_muSAD) * sqrt(N - 1);
  XR_muSAD = qr_thin_R(X_muSAD) / sqrt(N - 1);
  XR_muSAD_inv = inverse(XR_muSAD);

  // compute and scale QR decomposition
  XQ_muSUR = qr_thin_Q(X_muSUR) * sqrt(N - 1);
  XR_muSUR = qr_thin_R(X_muSUR) / sqrt(N - 1);
  XR_muSUR_inv = inverse(XR_muSUR);
}

Is this the correct way to simplify this?

transformed data {
  // matrices for QR decomposition
  matrix[N, K] XQ;
  matrix[K, K] XR;
  matrix[K, K] XR_inv;
  
  // compute and scale QR decomposition
  real sqrt_N_1 = sqrt(N - 1)
  XQ = qr_thin_Q(X) * sqrt_N_1;
  XR = qr_thin_R(X) / sqrt_N_1;
  XR_inv = inverse(XR);
}

Parameters

parameters {
  vector[K_muANG] b_muANG;  // population-level effects
  vector[K_muDIS] b_muDIS;  // population-level effects
  vector[K_muFER] b_muFER;  // population-level effects
  vector[K_muHAP] b_muHAP;  // population-level effects
  vector[K_muSAD] b_muSAD;  // population-level effects
  vector[K_muSUR] b_muSUR;  // population-level effects

  vector<lower=0>[M_1] sd_1;  // group-level standard deviations
  matrix[M_1, N_1] z_1;  // standardized group-level effects
  cholesky_factor_corr[M_1] L_1;  // cholesky factor of correlation matrix
  vector<lower=0>[M_2] sd_2;  // group-level standard deviations
  matrix[M_2, N_2] z_2;  // standardized group-level effects
  cholesky_factor_corr[M_2] L_2;  // cholesky factor of correlation matrix
  vector<lower=0>[M_3] sd_3;  // group-level standard deviations
  matrix[M_3, N_3] z_3;  // standardized group-level effects
  cholesky_factor_corr[M_3] L_3;  // cholesky factor of correlation matrix
  vector<lower=0>[M_4] sd_4;  // group-level standard deviations
  matrix[M_4, N_4] z_4;  // standardized group-level effects
  cholesky_factor_corr[M_4] L_4;  // cholesky factor of correlation matrix
  vector<lower=0>[M_5] sd_5;  // group-level standard deviations
  matrix[M_5, N_5] z_5;  // standardized group-level effects
  cholesky_factor_corr[M_5] L_5;  // cholesky factor of correlation matrix
  vector<lower=0>[M_6] sd_6;  // group-level standard deviations
  matrix[M_6, N_6] z_6;  // standardized group-level effects
  cholesky_factor_corr[M_6] L_6;  // cholesky factor of correlation matrix
  vector<lower=0>[M_7] sd_7;  // group-level standard deviations
  matrix[M_7, N_7] z_7;  // standardized group-level effects
  cholesky_factor_corr[M_7] L_7;  // cholesky factor of correlation matrix
  vector<lower=0>[M_8] sd_8;  // group-level standard deviations
  matrix[M_8, N_8] z_8;  // standardized group-level effects
  cholesky_factor_corr[M_8] L_8;  // cholesky factor of correlation matrix
  vector<lower=0>[M_9] sd_9;  // group-level standard deviations
  matrix[M_9, N_9] z_9;  // standardized group-level effects
  cholesky_factor_corr[M_9] L_9;  // cholesky factor of correlation matrix
  vector<lower=0>[M_10] sd_10;  // group-level standard deviations
  matrix[M_10, N_10] z_10;  // standardized group-level effects
  cholesky_factor_corr[M_10] L_10;  // cholesky factor of correlation matrix
  vector<lower=0>[M_11] sd_11;  // group-level standard deviations
  matrix[M_11, N_11] z_11;  // standardized group-level effects
  cholesky_factor_corr[M_11] L_11;  // cholesky factor of correlation matrix
  vector<lower=0>[M_12] sd_12;  // group-level standard deviations
  matrix[M_12, N_12] z_12;  // standardized group-level effects
  cholesky_factor_corr[M_12] L_12;  // cholesky factor of correlation matrix
  vector<lower=0>[M_13] sd_13;  // group-level standard deviations
  matrix[M_13, N_13] z_13;  // standardized group-level effects
  cholesky_factor_corr[M_13] L_13;  // cholesky factor of correlation matrix
  vector<lower=0>[M_14] sd_14;  // group-level standard deviations
  matrix[M_14, N_14] z_14;  // standardized group-level effects
  cholesky_factor_corr[M_14] L_14;  // cholesky factor of correlation matrix
  vector<lower=0>[M_15] sd_15;  // group-level standard deviations
  matrix[M_15, N_15] z_15;  // standardized group-level effects
  cholesky_factor_corr[M_15] L_15;  // cholesky factor of correlation matrix
  vector<lower=0>[M_16] sd_16;  // group-level standard deviations
  matrix[M_16, N_16] z_16;  // standardized group-level effects
  cholesky_factor_corr[M_16] L_16;  // cholesky factor of correlation matrix
  vector<lower=0>[M_17] sd_17;  // group-level standard deviations
  matrix[M_17, N_17] z_17;  // standardized group-level effects
  cholesky_factor_corr[M_17] L_17;  // cholesky factor of correlation matrix
  vector<lower=0>[M_18] sd_18;  // group-level standard deviations
  matrix[M_18, N_18] z_18;  // standardized group-level effects
  cholesky_factor_corr[M_18] L_18;  // cholesky factor of correlation matrix
}

We can rename some variables here, (e.g. due to the QR decomposition, the regression coefficients will be renamed, i.e., b_muANG to bQ_muANG, b_muDIS to bQ_muDIS, etc.), but I think the code cannot be further simplified. But, I think I – like addressed earlier – I am not understanding yet why we need to compute group-level correlations here.

Transformed parameters

transformed parameters {
  matrix[N_1, M_1] r_1;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_1] r_1_muANG_1;
  vector[N_1] r_1_muANG_2;
  vector[N_1] r_1_muANG_3;
  vector[N_1] r_1_muANG_4;
  vector[N_1] r_1_muANG_5;
  vector[N_1] r_1_muANG_6;
  vector[N_1] r_1_muANG_7;
  matrix[N_2, M_2] r_2;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_2] r_2_muANG_1;
  vector[N_2] r_2_muANG_2;
  vector[N_2] r_2_muANG_3;
  vector[N_2] r_2_muANG_4;
  vector[N_2] r_2_muANG_5;
  vector[N_2] r_2_muANG_6;
  vector[N_2] r_2_muANG_7;
  matrix[N_3, M_3] r_3;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_3] r_3_muANG_1;
  vector[N_3] r_3_muANG_2;
  vector[N_3] r_3_muANG_3;
  vector[N_3] r_3_muANG_4;
  vector[N_3] r_3_muANG_5;
  vector[N_3] r_3_muANG_6;
  vector[N_3] r_3_muANG_7;
  matrix[N_4, M_4] r_4;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_4] r_4_muDIS_1;
  vector[N_4] r_4_muDIS_2;
  vector[N_4] r_4_muDIS_3;
  vector[N_4] r_4_muDIS_4;
  vector[N_4] r_4_muDIS_5;
  vector[N_4] r_4_muDIS_6;
  vector[N_4] r_4_muDIS_7;
  matrix[N_5, M_5] r_5;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_5] r_5_muDIS_1;
  vector[N_5] r_5_muDIS_2;
  vector[N_5] r_5_muDIS_3;
  vector[N_5] r_5_muDIS_4;
  vector[N_5] r_5_muDIS_5;
  vector[N_5] r_5_muDIS_6;
  vector[N_5] r_5_muDIS_7;
  matrix[N_6, M_6] r_6;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_6] r_6_muDIS_1;
  vector[N_6] r_6_muDIS_2;
  vector[N_6] r_6_muDIS_3;
  vector[N_6] r_6_muDIS_4;
  vector[N_6] r_6_muDIS_5;
  vector[N_6] r_6_muDIS_6;
  vector[N_6] r_6_muDIS_7;
  matrix[N_7, M_7] r_7;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_7] r_7_muFER_1;
  vector[N_7] r_7_muFER_2;
  vector[N_7] r_7_muFER_3;
  vector[N_7] r_7_muFER_4;
  vector[N_7] r_7_muFER_5;
  vector[N_7] r_7_muFER_6;
  vector[N_7] r_7_muFER_7;
  matrix[N_8, M_8] r_8;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_8] r_8_muFER_1;
  vector[N_8] r_8_muFER_2;
  vector[N_8] r_8_muFER_3;
  vector[N_8] r_8_muFER_4;
  vector[N_8] r_8_muFER_5;
  vector[N_8] r_8_muFER_6;
  vector[N_8] r_8_muFER_7;
  matrix[N_9, M_9] r_9;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_9] r_9_muFER_1;
  vector[N_9] r_9_muFER_2;
  vector[N_9] r_9_muFER_3;
  vector[N_9] r_9_muFER_4;
  vector[N_9] r_9_muFER_5;
  vector[N_9] r_9_muFER_6;
  vector[N_9] r_9_muFER_7;
  matrix[N_10, M_10] r_10;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_10] r_10_muHAP_1;
  vector[N_10] r_10_muHAP_2;
  vector[N_10] r_10_muHAP_3;
  vector[N_10] r_10_muHAP_4;
  vector[N_10] r_10_muHAP_5;
  vector[N_10] r_10_muHAP_6;
  vector[N_10] r_10_muHAP_7;
  matrix[N_11, M_11] r_11;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_11] r_11_muHAP_1;
  vector[N_11] r_11_muHAP_2;
  vector[N_11] r_11_muHAP_3;
  vector[N_11] r_11_muHAP_4;
  vector[N_11] r_11_muHAP_5;
  vector[N_11] r_11_muHAP_6;
  vector[N_11] r_11_muHAP_7;
  matrix[N_12, M_12] r_12;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_12] r_12_muHAP_1;
  vector[N_12] r_12_muHAP_2;
  vector[N_12] r_12_muHAP_3;
  vector[N_12] r_12_muHAP_4;
  vector[N_12] r_12_muHAP_5;
  vector[N_12] r_12_muHAP_6;
  vector[N_12] r_12_muHAP_7;
  matrix[N_13, M_13] r_13;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_13] r_13_muSAD_1;
  vector[N_13] r_13_muSAD_2;
  vector[N_13] r_13_muSAD_3;
  vector[N_13] r_13_muSAD_4;
  vector[N_13] r_13_muSAD_5;
  vector[N_13] r_13_muSAD_6;
  vector[N_13] r_13_muSAD_7;
  matrix[N_14, M_14] r_14;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_14] r_14_muSAD_1;
  vector[N_14] r_14_muSAD_2;
  vector[N_14] r_14_muSAD_3;
  vector[N_14] r_14_muSAD_4;
  vector[N_14] r_14_muSAD_5;
  vector[N_14] r_14_muSAD_6;
  vector[N_14] r_14_muSAD_7;
  matrix[N_15, M_15] r_15;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_15] r_15_muSAD_1;
  vector[N_15] r_15_muSAD_2;
  vector[N_15] r_15_muSAD_3;
  vector[N_15] r_15_muSAD_4;
  vector[N_15] r_15_muSAD_5;
  vector[N_15] r_15_muSAD_6;
  vector[N_15] r_15_muSAD_7;
  matrix[N_16, M_16] r_16;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_16] r_16_muSUR_1;
  vector[N_16] r_16_muSUR_2;
  vector[N_16] r_16_muSUR_3;
  vector[N_16] r_16_muSUR_4;
  vector[N_16] r_16_muSUR_5;
  vector[N_16] r_16_muSUR_6;
  vector[N_16] r_16_muSUR_7;
  matrix[N_17, M_17] r_17;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_17] r_17_muSUR_1;
  vector[N_17] r_17_muSUR_2;
  vector[N_17] r_17_muSUR_3;
  vector[N_17] r_17_muSUR_4;
  vector[N_17] r_17_muSUR_5;
  vector[N_17] r_17_muSUR_6;
  vector[N_17] r_17_muSUR_7;
  matrix[N_18, M_18] r_18;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_18] r_18_muSUR_1;
  vector[N_18] r_18_muSUR_2;
  vector[N_18] r_18_muSUR_3;
  vector[N_18] r_18_muSUR_4;
  vector[N_18] r_18_muSUR_5;
  vector[N_18] r_18_muSUR_6;
  vector[N_18] r_18_muSUR_7;
  // compute actual group-level effects
  r_1 = scale_r_cor(z_1, sd_1, L_1);
  r_1_muANG_1 = r_1[, 1];
  r_1_muANG_2 = r_1[, 2];
  r_1_muANG_3 = r_1[, 3];
  r_1_muANG_4 = r_1[, 4];
  r_1_muANG_5 = r_1[, 5];
  r_1_muANG_6 = r_1[, 6];
  r_1_muANG_7 = r_1[, 7];
  // compute actual group-level effects
  r_2 = scale_r_cor(z_2, sd_2, L_2);
  r_2_muANG_1 = r_2[, 1];
  r_2_muANG_2 = r_2[, 2];
  r_2_muANG_3 = r_2[, 3];
  r_2_muANG_4 = r_2[, 4];
  r_2_muANG_5 = r_2[, 5];
  r_2_muANG_6 = r_2[, 6];
  r_2_muANG_7 = r_2[, 7];
  // compute actual group-level effects
  r_3 = scale_r_cor(z_3, sd_3, L_3);
  r_3_muANG_1 = r_3[, 1];
  r_3_muANG_2 = r_3[, 2];
  r_3_muANG_3 = r_3[, 3];
  r_3_muANG_4 = r_3[, 4];
  r_3_muANG_5 = r_3[, 5];
  r_3_muANG_6 = r_3[, 6];
  r_3_muANG_7 = r_3[, 7];
  // compute actual group-level effects
  r_4 = scale_r_cor(z_4, sd_4, L_4);
  r_4_muDIS_1 = r_4[, 1];
  r_4_muDIS_2 = r_4[, 2];
  r_4_muDIS_3 = r_4[, 3];
  r_4_muDIS_4 = r_4[, 4];
  r_4_muDIS_5 = r_4[, 5];
  r_4_muDIS_6 = r_4[, 6];
  r_4_muDIS_7 = r_4[, 7];
  // compute actual group-level effects
  r_5 = scale_r_cor(z_5, sd_5, L_5);
  r_5_muDIS_1 = r_5[, 1];
  r_5_muDIS_2 = r_5[, 2];
  r_5_muDIS_3 = r_5[, 3];
  r_5_muDIS_4 = r_5[, 4];
  r_5_muDIS_5 = r_5[, 5];
  r_5_muDIS_6 = r_5[, 6];
  r_5_muDIS_7 = r_5[, 7];
  // compute actual group-level effects
  r_6 = scale_r_cor(z_6, sd_6, L_6);
  r_6_muDIS_1 = r_6[, 1];
  r_6_muDIS_2 = r_6[, 2];
  r_6_muDIS_3 = r_6[, 3];
  r_6_muDIS_4 = r_6[, 4];
  r_6_muDIS_5 = r_6[, 5];
  r_6_muDIS_6 = r_6[, 6];
  r_6_muDIS_7 = r_6[, 7];
  // compute actual group-level effects
  r_7 = scale_r_cor(z_7, sd_7, L_7);
  r_7_muFER_1 = r_7[, 1];
  r_7_muFER_2 = r_7[, 2];
  r_7_muFER_3 = r_7[, 3];
  r_7_muFER_4 = r_7[, 4];
  r_7_muFER_5 = r_7[, 5];
  r_7_muFER_6 = r_7[, 6];
  r_7_muFER_7 = r_7[, 7];
  // compute actual group-level effects
  r_8 = scale_r_cor(z_8, sd_8, L_8);
  r_8_muFER_1 = r_8[, 1];
  r_8_muFER_2 = r_8[, 2];
  r_8_muFER_3 = r_8[, 3];
  r_8_muFER_4 = r_8[, 4];
  r_8_muFER_5 = r_8[, 5];
  r_8_muFER_6 = r_8[, 6];
  r_8_muFER_7 = r_8[, 7];
  // compute actual group-level effects
  r_9 = scale_r_cor(z_9, sd_9, L_9);
  r_9_muFER_1 = r_9[, 1];
  r_9_muFER_2 = r_9[, 2];
  r_9_muFER_3 = r_9[, 3];
  r_9_muFER_4 = r_9[, 4];
  r_9_muFER_5 = r_9[, 5];
  r_9_muFER_6 = r_9[, 6];
  r_9_muFER_7 = r_9[, 7];
  // compute actual group-level effects
  r_10 = scale_r_cor(z_10, sd_10, L_10);
  r_10_muHAP_1 = r_10[, 1];
  r_10_muHAP_2 = r_10[, 2];
  r_10_muHAP_3 = r_10[, 3];
  r_10_muHAP_4 = r_10[, 4];
  r_10_muHAP_5 = r_10[, 5];
  r_10_muHAP_6 = r_10[, 6];
  r_10_muHAP_7 = r_10[, 7];
  // compute actual group-level effects
  r_11 = scale_r_cor(z_11, sd_11, L_11);
  r_11_muHAP_1 = r_11[, 1];
  r_11_muHAP_2 = r_11[, 2];
  r_11_muHAP_3 = r_11[, 3];
  r_11_muHAP_4 = r_11[, 4];
  r_11_muHAP_5 = r_11[, 5];
  r_11_muHAP_6 = r_11[, 6];
  r_11_muHAP_7 = r_11[, 7];
  // compute actual group-level effects
  r_12 = scale_r_cor(z_12, sd_12, L_12);
  r_12_muHAP_1 = r_12[, 1];
  r_12_muHAP_2 = r_12[, 2];
  r_12_muHAP_3 = r_12[, 3];
  r_12_muHAP_4 = r_12[, 4];
  r_12_muHAP_5 = r_12[, 5];
  r_12_muHAP_6 = r_12[, 6];
  r_12_muHAP_7 = r_12[, 7];
  // compute actual group-level effects
  r_13 = scale_r_cor(z_13, sd_13, L_13);
  r_13_muSAD_1 = r_13[, 1];
  r_13_muSAD_2 = r_13[, 2];
  r_13_muSAD_3 = r_13[, 3];
  r_13_muSAD_4 = r_13[, 4];
  r_13_muSAD_5 = r_13[, 5];
  r_13_muSAD_6 = r_13[, 6];
  r_13_muSAD_7 = r_13[, 7];
  // compute actual group-level effects
  r_14 = scale_r_cor(z_14, sd_14, L_14);
  r_14_muSAD_1 = r_14[, 1];
  r_14_muSAD_2 = r_14[, 2];
  r_14_muSAD_3 = r_14[, 3];
  r_14_muSAD_4 = r_14[, 4];
  r_14_muSAD_5 = r_14[, 5];
  r_14_muSAD_6 = r_14[, 6];
  r_14_muSAD_7 = r_14[, 7];
  // compute actual group-level effects
  r_15 = scale_r_cor(z_15, sd_15, L_15);
  r_15_muSAD_1 = r_15[, 1];
  r_15_muSAD_2 = r_15[, 2];
  r_15_muSAD_3 = r_15[, 3];
  r_15_muSAD_4 = r_15[, 4];
  r_15_muSAD_5 = r_15[, 5];
  r_15_muSAD_6 = r_15[, 6];
  r_15_muSAD_7 = r_15[, 7];
  // compute actual group-level effects
  r_16 = scale_r_cor(z_16, sd_16, L_16);
  r_16_muSUR_1 = r_16[, 1];
  r_16_muSUR_2 = r_16[, 2];
  r_16_muSUR_3 = r_16[, 3];
  r_16_muSUR_4 = r_16[, 4];
  r_16_muSUR_5 = r_16[, 5];
  r_16_muSUR_6 = r_16[, 6];
  r_16_muSUR_7 = r_16[, 7];
  // compute actual group-level effects
  r_17 = scale_r_cor(z_17, sd_17, L_17);
  r_17_muSUR_1 = r_17[, 1];
  r_17_muSUR_2 = r_17[, 2];
  r_17_muSUR_3 = r_17[, 3];
  r_17_muSUR_4 = r_17[, 4];
  r_17_muSUR_5 = r_17[, 5];
  r_17_muSUR_6 = r_17[, 6];
  r_17_muSUR_7 = r_17[, 7];
  // compute actual group-level effects
  r_18 = scale_r_cor(z_18, sd_18, L_18);
  r_18_muSUR_1 = r_18[, 1];
  r_18_muSUR_2 = r_18[, 2];
  r_18_muSUR_3 = r_18[, 3];
  r_18_muSUR_4 = r_18[, 4];
  r_18_muSUR_5 = r_18[, 5];
  r_18_muSUR_6 = r_18[, 6];
  r_18_muSUR_7 = r_18[, 7];
}

As in parameters, there is not much to do except renaming (e.g., r_1_muANG_1 to r_culture_muANG_RC1) also maybe variable declaration and assignment can be combined, so instead of:

matrix[N_1, M_1] r_culture_ANG;
r_culture_ANG = scale_r_cor(z_culture_ANG, sd_culture_ANG, L_culture_ANG);

One could write:

matrix[N_1, M_1] r_culture_ANG = scale_r_cor(z_culture_ANG, sd_culture_ANG, L_culture_ANG);

Model

model {
  // likelihood including all constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N] muANG = X_muANG * b_muANG;
    // initialize linear predictor term
    vector[N] muDIS = X_muDIS * b_muDIS;
    // initialize linear predictor term
    vector[N] muFER = X_muFER * b_muFER;
    // initialize linear predictor term
    vector[N] muHAP = X_muHAP * b_muHAP;
    // initialize linear predictor term
    vector[N] muSAD = X_muSAD * b_muSAD;
    // initialize linear predictor term
    vector[N] muSUR = X_muSUR * b_muSUR;
    // linear predictor matrix
    vector[ncat] mu[N];
    for (n in 1:N) {
      // add more terms to the linear predictor
      muANG[n] += r_1_muANG_1[J_1[n]] * Z_1_muANG_1[n] + r_1_muANG_2[J_1[n]] * Z_1_muANG_2[n] + r_1_muANG_3[J_1[n]] * Z_1_muANG_3[n] + r_1_muANG_4[J_1[n]] * Z_1_muANG_4[n] + r_1_muANG_5[J_1[n]] * Z_1_muANG_5[n] + r_1_muANG_6[J_1[n]] * Z_1_muANG_6[n] + r_1_muANG_7[J_1[n]] * Z_1_muANG_7[n] + r_2_muANG_1[J_2[n]] * Z_2_muANG_1[n] + r_2_muANG_2[J_2[n]] * Z_2_muANG_2[n] + r_2_muANG_3[J_2[n]] * Z_2_muANG_3[n] + r_2_muANG_4[J_2[n]] * Z_2_muANG_4[n] + r_2_muANG_5[J_2[n]] * Z_2_muANG_5[n] + r_2_muANG_6[J_2[n]] * Z_2_muANG_6[n] + r_2_muANG_7[J_2[n]] * Z_2_muANG_7[n] + r_3_muANG_1[J_3[n]] * Z_3_muANG_1[n] + r_3_muANG_2[J_3[n]] * Z_3_muANG_2[n] + r_3_muANG_3[J_3[n]] * Z_3_muANG_3[n] + r_3_muANG_4[J_3[n]] * Z_3_muANG_4[n] + r_3_muANG_5[J_3[n]] * Z_3_muANG_5[n] + r_3_muANG_6[J_3[n]] * Z_3_muANG_6[n] + r_3_muANG_7[J_3[n]] * Z_3_muANG_7[n];
    }
    for (n in 1:N) {
      // add more terms to the linear predictor
      muDIS[n] += r_4_muDIS_1[J_4[n]] * Z_4_muDIS_1[n] + r_4_muDIS_2[J_4[n]] * Z_4_muDIS_2[n] + r_4_muDIS_3[J_4[n]] * Z_4_muDIS_3[n] + r_4_muDIS_4[J_4[n]] * Z_4_muDIS_4[n] + r_4_muDIS_5[J_4[n]] * Z_4_muDIS_5[n] + r_4_muDIS_6[J_4[n]] * Z_4_muDIS_6[n] + r_4_muDIS_7[J_4[n]] * Z_4_muDIS_7[n] + r_5_muDIS_1[J_5[n]] * Z_5_muDIS_1[n] + r_5_muDIS_2[J_5[n]] * Z_5_muDIS_2[n] + r_5_muDIS_3[J_5[n]] * Z_5_muDIS_3[n] + r_5_muDIS_4[J_5[n]] * Z_5_muDIS_4[n] + r_5_muDIS_5[J_5[n]] * Z_5_muDIS_5[n] + r_5_muDIS_6[J_5[n]] * Z_5_muDIS_6[n] + r_5_muDIS_7[J_5[n]] * Z_5_muDIS_7[n] + r_6_muDIS_1[J_6[n]] * Z_6_muDIS_1[n] + r_6_muDIS_2[J_6[n]] * Z_6_muDIS_2[n] + r_6_muDIS_3[J_6[n]] * Z_6_muDIS_3[n] + r_6_muDIS_4[J_6[n]] * Z_6_muDIS_4[n] + r_6_muDIS_5[J_6[n]] * Z_6_muDIS_5[n] + r_6_muDIS_6[J_6[n]] * Z_6_muDIS_6[n] + r_6_muDIS_7[J_6[n]] * Z_6_muDIS_7[n];
    }
    for (n in 1:N) {
      // add more terms to the linear predictor
      muFER[n] += r_7_muFER_1[J_7[n]] * Z_7_muFER_1[n] + r_7_muFER_2[J_7[n]] * Z_7_muFER_2[n] + r_7_muFER_3[J_7[n]] * Z_7_muFER_3[n] + r_7_muFER_4[J_7[n]] * Z_7_muFER_4[n] + r_7_muFER_5[J_7[n]] * Z_7_muFER_5[n] + r_7_muFER_6[J_7[n]] * Z_7_muFER_6[n] + r_7_muFER_7[J_7[n]] * Z_7_muFER_7[n] + r_8_muFER_1[J_8[n]] * Z_8_muFER_1[n] + r_8_muFER_2[J_8[n]] * Z_8_muFER_2[n] + r_8_muFER_3[J_8[n]] * Z_8_muFER_3[n] + r_8_muFER_4[J_8[n]] * Z_8_muFER_4[n] + r_8_muFER_5[J_8[n]] * Z_8_muFER_5[n] + r_8_muFER_6[J_8[n]] * Z_8_muFER_6[n] + r_8_muFER_7[J_8[n]] * Z_8_muFER_7[n] + r_9_muFER_1[J_9[n]] * Z_9_muFER_1[n] + r_9_muFER_2[J_9[n]] * Z_9_muFER_2[n] + r_9_muFER_3[J_9[n]] * Z_9_muFER_3[n] + r_9_muFER_4[J_9[n]] * Z_9_muFER_4[n] + r_9_muFER_5[J_9[n]] * Z_9_muFER_5[n] + r_9_muFER_6[J_9[n]] * Z_9_muFER_6[n] + r_9_muFER_7[J_9[n]] * Z_9_muFER_7[n];
    }
    for (n in 1:N) {
      // add more terms to the linear predictor
      muHAP[n] += r_10_muHAP_1[J_10[n]] * Z_10_muHAP_1[n] + r_10_muHAP_2[J_10[n]] * Z_10_muHAP_2[n] + r_10_muHAP_3[J_10[n]] * Z_10_muHAP_3[n] + r_10_muHAP_4[J_10[n]] * Z_10_muHAP_4[n] + r_10_muHAP_5[J_10[n]] * Z_10_muHAP_5[n] + r_10_muHAP_6[J_10[n]] * Z_10_muHAP_6[n] + r_10_muHAP_7[J_10[n]] * Z_10_muHAP_7[n] + r_11_muHAP_1[J_11[n]] * Z_11_muHAP_1[n] + r_11_muHAP_2[J_11[n]] * Z_11_muHAP_2[n] + r_11_muHAP_3[J_11[n]] * Z_11_muHAP_3[n] + r_11_muHAP_4[J_11[n]] * Z_11_muHAP_4[n] + r_11_muHAP_5[J_11[n]] * Z_11_muHAP_5[n] + r_11_muHAP_6[J_11[n]] * Z_11_muHAP_6[n] + r_11_muHAP_7[J_11[n]] * Z_11_muHAP_7[n] + r_12_muHAP_1[J_12[n]] * Z_12_muHAP_1[n] + r_12_muHAP_2[J_12[n]] * Z_12_muHAP_2[n] + r_12_muHAP_3[J_12[n]] * Z_12_muHAP_3[n] + r_12_muHAP_4[J_12[n]] * Z_12_muHAP_4[n] + r_12_muHAP_5[J_12[n]] * Z_12_muHAP_5[n] + r_12_muHAP_6[J_12[n]] * Z_12_muHAP_6[n] + r_12_muHAP_7[J_12[n]] * Z_12_muHAP_7[n];
    }
    for (n in 1:N) {
      // add more terms to the linear predictor
      muSAD[n] += r_13_muSAD_1[J_13[n]] * Z_13_muSAD_1[n] + r_13_muSAD_2[J_13[n]] * Z_13_muSAD_2[n] + r_13_muSAD_3[J_13[n]] * Z_13_muSAD_3[n] + r_13_muSAD_4[J_13[n]] * Z_13_muSAD_4[n] + r_13_muSAD_5[J_13[n]] * Z_13_muSAD_5[n] + r_13_muSAD_6[J_13[n]] * Z_13_muSAD_6[n] + r_13_muSAD_7[J_13[n]] * Z_13_muSAD_7[n] + r_14_muSAD_1[J_14[n]] * Z_14_muSAD_1[n] + r_14_muSAD_2[J_14[n]] * Z_14_muSAD_2[n] + r_14_muSAD_3[J_14[n]] * Z_14_muSAD_3[n] + r_14_muSAD_4[J_14[n]] * Z_14_muSAD_4[n] + r_14_muSAD_5[J_14[n]] * Z_14_muSAD_5[n] + r_14_muSAD_6[J_14[n]] * Z_14_muSAD_6[n] + r_14_muSAD_7[J_14[n]] * Z_14_muSAD_7[n] + r_15_muSAD_1[J_15[n]] * Z_15_muSAD_1[n] + r_15_muSAD_2[J_15[n]] * Z_15_muSAD_2[n] + r_15_muSAD_3[J_15[n]] * Z_15_muSAD_3[n] + r_15_muSAD_4[J_15[n]] * Z_15_muSAD_4[n] + r_15_muSAD_5[J_15[n]] * Z_15_muSAD_5[n] + r_15_muSAD_6[J_15[n]] * Z_15_muSAD_6[n] + r_15_muSAD_7[J_15[n]] * Z_15_muSAD_7[n];
    }
    for (n in 1:N) {
      // add more terms to the linear predictor
      muSUR[n] += r_16_muSUR_1[J_16[n]] * Z_16_muSUR_1[n] + r_16_muSUR_2[J_16[n]] * Z_16_muSUR_2[n] + r_16_muSUR_3[J_16[n]] * Z_16_muSUR_3[n] + r_16_muSUR_4[J_16[n]] * Z_16_muSUR_4[n] + r_16_muSUR_5[J_16[n]] * Z_16_muSUR_5[n] + r_16_muSUR_6[J_16[n]] * Z_16_muSUR_6[n] + r_16_muSUR_7[J_16[n]] * Z_16_muSUR_7[n] + r_17_muSUR_1[J_17[n]] * Z_17_muSUR_1[n] + r_17_muSUR_2[J_17[n]] * Z_17_muSUR_2[n] + r_17_muSUR_3[J_17[n]] * Z_17_muSUR_3[n] + r_17_muSUR_4[J_17[n]] * Z_17_muSUR_4[n] + r_17_muSUR_5[J_17[n]] * Z_17_muSUR_5[n] + r_17_muSUR_6[J_17[n]] * Z_17_muSUR_6[n] + r_17_muSUR_7[J_17[n]] * Z_17_muSUR_7[n] + r_18_muSUR_1[J_18[n]] * Z_18_muSUR_1[n] + r_18_muSUR_2[J_18[n]] * Z_18_muSUR_2[n] + r_18_muSUR_3[J_18[n]] * Z_18_muSUR_3[n] + r_18_muSUR_4[J_18[n]] * Z_18_muSUR_4[n] + r_18_muSUR_5[J_18[n]] * Z_18_muSUR_5[n] + r_18_muSUR_6[J_18[n]] * Z_18_muSUR_6[n] + r_18_muSUR_7[J_18[n]] * Z_18_muSUR_7[n];
    }
    for (n in 1:N) {
      mu[n] = transpose([0, muANG[n], muDIS[n], muFER[n], muHAP[n], muSAD[n], muSUR[n]]);
    }
    for (n in 1:N) {
      target += categorical_logit_lpmf(Y[n] | mu[n]);
    }
  }
  // priors including all constants
  target += normal_lpdf(b_muANG | 0,1);
  target += normal_lpdf(b_muDIS | 0,1);
  target += normal_lpdf(b_muFER | 0,1);
  target += normal_lpdf(b_muHAP | 0,1);
  target += normal_lpdf(b_muSAD | 0,1);
  target += normal_lpdf(b_muSUR | 0,1);
  target += student_t_lpdf(sd_1 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_1));
  target += lkj_corr_cholesky_lpdf(L_1 | 1);
  target += student_t_lpdf(sd_2 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_2));
  target += lkj_corr_cholesky_lpdf(L_2 | 1);
  target += student_t_lpdf(sd_3 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_3));
  target += lkj_corr_cholesky_lpdf(L_3 | 1);
  target += student_t_lpdf(sd_4 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_4));
  target += lkj_corr_cholesky_lpdf(L_4 | 1);
  target += student_t_lpdf(sd_5 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_5));
  target += lkj_corr_cholesky_lpdf(L_5 | 1);
  target += student_t_lpdf(sd_6 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_6));
  target += lkj_corr_cholesky_lpdf(L_6 | 1);
  target += student_t_lpdf(sd_7 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_7));
  target += lkj_corr_cholesky_lpdf(L_7 | 1);
  target += student_t_lpdf(sd_8 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_8));
  target += lkj_corr_cholesky_lpdf(L_8 | 1);
  target += student_t_lpdf(sd_9 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_9));
  target += lkj_corr_cholesky_lpdf(L_9 | 1);
  target += student_t_lpdf(sd_10 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_10));
  target += lkj_corr_cholesky_lpdf(L_10 | 1);
  target += student_t_lpdf(sd_11 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_11));
  target += lkj_corr_cholesky_lpdf(L_11 | 1);
  target += student_t_lpdf(sd_12 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_12));
  target += lkj_corr_cholesky_lpdf(L_12 | 1);
  target += student_t_lpdf(sd_13 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_13));
  target += lkj_corr_cholesky_lpdf(L_13 | 1);
  target += student_t_lpdf(sd_14 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_14));
  target += lkj_corr_cholesky_lpdf(L_14 | 1);
  target += student_t_lpdf(sd_15 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_15));
  target += lkj_corr_cholesky_lpdf(L_15 | 1);
  target += student_t_lpdf(sd_16 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_16));
  target += lkj_corr_cholesky_lpdf(L_16 | 1);
  target += student_t_lpdf(sd_17 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_17));
  target += lkj_corr_cholesky_lpdf(L_17 | 1);
  target += student_t_lpdf(sd_18 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);
  target += std_normal_lpdf(to_vector(z_18));
  target += lkj_corr_cholesky_lpdf(L_18 | 1);
}

QR decomposition

This will change the linear predictor term

    // initialize linear predictor term
    vector[N] muANG = XQ * bQ_muANG;
    vector[N] muDIS = XQ * bQ_muDIS;
    vector[N] muFER = XQ * bQ_muFER;
    vector[N] muHAP = XQ * bQ_muHAP;
    vector[N] muSAD = XQ * bQ_muSAD;
    vector[N] muSUR = XQ * bQ_muSUR;

Vectorize for-loops

We can simplify

    // initialize linear predictor term
    vector[N] muANG = XQ * bQ_muANG;
    for (n in 1:N) {
      // add more terms to the linear predictor
      muANG[n] += r_1_muANG_1[J_1[n]] * Z_1_muANG_1[n] + r_1_muANG_2[J_1[n]] * Z_1_muANG_2[n] + r_1_muANG_3[J_1[n]] * Z_1_muANG_3[n] + r_1_muANG_4[J_1[n]] * Z_1_muANG_4[n] + r_1_muANG_5[J_1[n]] * Z_1_muANG_5[n] + r_1_muANG_6[J_1[n]] * Z_1_muANG_6[n] + r_1_muANG_7[J_1[n]] * Z_1_muANG_7[n] + r_2_muANG_1[J_2[n]] * Z_2_muANG_1[n] + r_2_muANG_2[J_2[n]] * Z_2_muANG_2[n] + r_2_muANG_3[J_2[n]] * Z_2_muANG_3[n] + r_2_muANG_4[J_2[n]] * Z_2_muANG_4[n] + r_2_muANG_5[J_2[n]] * Z_2_muANG_5[n] + r_2_muANG_6[J_2[n]] * Z_2_muANG_6[n] + r_2_muANG_7[J_2[n]] * Z_2_muANG_7[n] + r_3_muANG_1[J_3[n]] * Z_3_muANG_1[n] + r_3_muANG_2[J_3[n]] * Z_3_muANG_2[n] + r_3_muANG_3[J_3[n]] * Z_3_muANG_3[n] + r_3_muANG_4[J_3[n]] * Z_3_muANG_4[n] + r_3_muANG_5[J_3[n]] * Z_3_muANG_5[n] + r_3_muANG_6[J_3[n]] * Z_3_muANG_6[n] + r_3_muANG_7[J_3[n]] * Z_3_muANG_7[n];
    }

to

// ANG linear predictor term
vector[N] muANG = XQ * bQ_muANG
                + r_culture_muANG[J_culture] .* Z_culture
                + r_sex_muANG[J_sex] .* Z_sex
                + r_speaker_muANG[J_speaker] .* Z_speaker;

Vectorizing mu for-loop

Also, we can save one for loop for computing mu, instead of writing

    for (n in 1:N) {
      mu[n] = transpose([0, muANG[n], muDIS[n], muFER[n], muHAP[n], muSAD[n], muSUR[n]]);
    }
    for (n in 1:N) {
      target += categorical_logit_lpmf(Y[n] | mu[n]);
    }

We can write:

    for (n in 1:N) {
      mu[n] = transpose([0, muANG[n], muDIS[n], muFER[n], muHAP[n], muSAD[n], muSUR[n]]);
      target += categorical_logit_lpmf(Y[n] | mu[n]);
    }

Priors

Instead of:

 // priors including all constants
  target += normal_lpdf(b_muANG | 0,1);
  target += normal_lpdf(b_muDIS | 0,1);
  target += normal_lpdf(b_muFER | 0,1);
  target += normal_lpdf(b_muHAP | 0,1);
  target += normal_lpdf(b_muSAD | 0,1);
  target += normal_lpdf(b_muSUR | 0,1);

We write (b_muANG to bQ_muANG because of QR decomposition):

  // priors including all constants
  bQ_muANG ~ std_normal();
  bQ_muDIS ~ std_normal();
  bQ_muFER ~ std_normal();
  bQ_muHAP ~ std_normal();
  bQ_muSAD ~ std_normal();
  bQ_muSUR ~ std_normal();

@storopoli, if I understand the ticket and PR correctly, each of these blocks

  target += student_t_lpdf(sd_1 | 3, 0, 2.5)
    - 7 * student_t_lccdf(0 | 3, 0, 2.5);

could be rewritten. However, I am not sure if this is the correct way to rewrite this. Is it?

  target += student_t_lupdf(sd_1 | 3, 0, 2.5);

Also should I convert

  target += std_normal_lpdf(to_vector(z_1));
  target += lkj_corr_cholesky_lpdf(L_1 | 1);

to

  target += std_normal_lupdf(to_vector(z_1));
  target += lkj_corr_cholesky_lupdf(L_1 | 1);

?

Part 3:

Generated quantities

generated quantities {
  // compute group-level correlations
  corr_matrix[M_1] Cor_1 = multiply_lower_tri_self_transpose(L_1);
  vector<lower=-1,upper=1>[NC_1] cor_1;
  // compute group-level correlations
  corr_matrix[M_2] Cor_2 = multiply_lower_tri_self_transpose(L_2);
  vector<lower=-1,upper=1>[NC_2] cor_2;
  // compute group-level correlations
  corr_matrix[M_3] Cor_3 = multiply_lower_tri_self_transpose(L_3);
  vector<lower=-1,upper=1>[NC_3] cor_3;
  // compute group-level correlations
  corr_matrix[M_4] Cor_4 = multiply_lower_tri_self_transpose(L_4);
  vector<lower=-1,upper=1>[NC_4] cor_4;
  // compute group-level correlations
  corr_matrix[M_5] Cor_5 = multiply_lower_tri_self_transpose(L_5);
  vector<lower=-1,upper=1>[NC_5] cor_5;
  // compute group-level correlations
  corr_matrix[M_6] Cor_6 = multiply_lower_tri_self_transpose(L_6);
  vector<lower=-1,upper=1>[NC_6] cor_6;
  // compute group-level correlations
  corr_matrix[M_7] Cor_7 = multiply_lower_tri_self_transpose(L_7);
  vector<lower=-1,upper=1>[NC_7] cor_7;
  // compute group-level correlations
  corr_matrix[M_8] Cor_8 = multiply_lower_tri_self_transpose(L_8);
  vector<lower=-1,upper=1>[NC_8] cor_8;
  // compute group-level correlations
  corr_matrix[M_9] Cor_9 = multiply_lower_tri_self_transpose(L_9);
  vector<lower=-1,upper=1>[NC_9] cor_9;
  // compute group-level correlations
  corr_matrix[M_10] Cor_10 = multiply_lower_tri_self_transpose(L_10);
  vector<lower=-1,upper=1>[NC_10] cor_10;
  // compute group-level correlations
  corr_matrix[M_11] Cor_11 = multiply_lower_tri_self_transpose(L_11);
  vector<lower=-1,upper=1>[NC_11] cor_11;
  // compute group-level correlations
  corr_matrix[M_12] Cor_12 = multiply_lower_tri_self_transpose(L_12);
  vector<lower=-1,upper=1>[NC_12] cor_12;
  // compute group-level correlations
  corr_matrix[M_13] Cor_13 = multiply_lower_tri_self_transpose(L_13);
  vector<lower=-1,upper=1>[NC_13] cor_13;
  // compute group-level correlations
  corr_matrix[M_14] Cor_14 = multiply_lower_tri_self_transpose(L_14);
  vector<lower=-1,upper=1>[NC_14] cor_14;
  // compute group-level correlations
  corr_matrix[M_15] Cor_15 = multiply_lower_tri_self_transpose(L_15);
  vector<lower=-1,upper=1>[NC_15] cor_15;
  // compute group-level correlations
  corr_matrix[M_16] Cor_16 = multiply_lower_tri_self_transpose(L_16);
  vector<lower=-1,upper=1>[NC_16] cor_16;
  // compute group-level correlations
  corr_matrix[M_17] Cor_17 = multiply_lower_tri_self_transpose(L_17);
  vector<lower=-1,upper=1>[NC_17] cor_17;
  // compute group-level correlations
  corr_matrix[M_18] Cor_18 = multiply_lower_tri_self_transpose(L_18);
  vector<lower=-1,upper=1>[NC_18] cor_18;
  // extract upper diagonal of correlation matrix
  for (k in 1:M_1) {
    for (j in 1:(k - 1)) {
      cor_1[choose(k - 1, 2) + j] = Cor_1[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_2) {
    for (j in 1:(k - 1)) {
      cor_2[choose(k - 1, 2) + j] = Cor_2[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_3) {
    for (j in 1:(k - 1)) {
      cor_3[choose(k - 1, 2) + j] = Cor_3[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_4) {
    for (j in 1:(k - 1)) {
      cor_4[choose(k - 1, 2) + j] = Cor_4[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_5) {
    for (j in 1:(k - 1)) {
      cor_5[choose(k - 1, 2) + j] = Cor_5[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_6) {
    for (j in 1:(k - 1)) {
      cor_6[choose(k - 1, 2) + j] = Cor_6[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_7) {
    for (j in 1:(k - 1)) {
      cor_7[choose(k - 1, 2) + j] = Cor_7[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_8) {
    for (j in 1:(k - 1)) {
      cor_8[choose(k - 1, 2) + j] = Cor_8[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_9) {
    for (j in 1:(k - 1)) {
      cor_9[choose(k - 1, 2) + j] = Cor_9[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_10) {
    for (j in 1:(k - 1)) {
      cor_10[choose(k - 1, 2) + j] = Cor_10[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_11) {
    for (j in 1:(k - 1)) {
      cor_11[choose(k - 1, 2) + j] = Cor_11[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_12) {
    for (j in 1:(k - 1)) {
      cor_12[choose(k - 1, 2) + j] = Cor_12[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_13) {
    for (j in 1:(k - 1)) {
      cor_13[choose(k - 1, 2) + j] = Cor_13[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_14) {
    for (j in 1:(k - 1)) {
      cor_14[choose(k - 1, 2) + j] = Cor_14[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_15) {
    for (j in 1:(k - 1)) {
      cor_15[choose(k - 1, 2) + j] = Cor_15[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_16) {
    for (j in 1:(k - 1)) {
      cor_16[choose(k - 1, 2) + j] = Cor_16[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_17) {
    for (j in 1:(k - 1)) {
      cor_17[choose(k - 1, 2) + j] = Cor_17[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_18) {
    for (j in 1:(k - 1)) {
      cor_18[choose(k - 1, 2) + j] = Cor_18[j, k];
    }
  }
}

The generated quantities block contains the computation of the group-level correlations. Because of the QR decomposition, we would also add this code to the block:

  // obtain the actual coefficients
  vector[K_muANG] b_muANG = XR_inv * bQ_muANG;
  // obtain the actual coefficients
  vector[K_muDIS] b_muDIS = XR_inv * bQ_muDIS;
  // obtain the actual coefficients
  vector[K_muFER] b_muFER = XR_inv * bQ_muFER;
  // obtain the actual coefficients
  vector[K_muHAP] b_muHAP = XR_inv * bQ_muHAP;
  // obtain the actual coefficients
  vector[K_muSAD] b_muSAD = XR_inv * bQ_muSAD;
  // obtain the actual coefficients
  vector[K_muSUR] b_muSUR = XR_inv * bQ_muSUR;

I haven’t tried using CdmStanR as a backend yet, nor have I applied within-chain parallelization yet, but I will do this once the model has been rewritten.

I would also be happy to write a vignette about it once we converged on a solution.

Thanks again for your help and suggestions and I am looking forward to your reply.

In brms X is for the population level stuff only. The Zs are the equivalent multipliers for the hierarchical terms.

So like if we had a term (my_covariate | group), the my_covariate terms would be packed into a Z coefficient.

You can leave the 7 * student_t_lccdf(0 | 3, 0, 2.5) off since it’s a constant term.

lpdf and lupdf are equivalent in that context.

As far as the QR, the way to check this is probably by comparing results from an old version of your model with a new version of your model. The generic way to check models is SBC (here).

It’s hard to check this stuff just by manually reading the code – the easiest way to debug these things is to compare a model you trust vs. the new model.