// generated with brms 2.16.1
functions {
 /* compute correlated group-level effects with 'by' variables
  * in the presence of a within-group covariance matrix
  * Args: 
  *   z: matrix of unscaled group-level effects
  *   SD: matrix of standard deviation parameters
  *   L: an array of cholesky factor correlation matrices
  *   Jby: index which grouping level belongs to which by level
  *   Lcov: cholesky factor of within-group correlation matrix
  * Returns: 
  *   matrix of scaled group-level effects
  */ 
  matrix scale_r_cor_by_cov(matrix z, matrix SD, matrix[] L, 
                            int[] Jby, matrix Lcov) {
    vector[num_elements(z)] z_flat = to_vector(z);
    vector[num_elements(z)] r = rep_vector(0, num_elements(z));
    matrix[rows(L[1]), cols(L[1])] LC[size(L)];
    int rows_z = rows(z);
    int rows_L = rows(L[1]);
    for (i in 1:size(LC)) {
      LC[i] = diag_pre_multiply(SD[, i], L[i]);
    }
    // kronecker product of cholesky factors times a vector
    for (icov in 1:rows(Lcov)) {
      for (jcov in 1:icov) {
        if (Lcov[icov, jcov] > 1e-10) { 
          // avoid calculating products between unrelated individuals
          for (i in 1:rows_L) {
            for (j in 1:i) {
              // incremented element of the output vector
              int k = (rows_L * (icov - 1)) + i;
              // applied element of the input vector
              int l = (rows_L * (jcov - 1)) + j;
              // column number of z to which z_flat[l] belongs
              int m = (l - 1) / rows_z + 1;
              r[k] = r[k] + Lcov[icov, jcov] * LC[Jby[m]][i, j] * z_flat[l];
            }
          }
        }
      }
    }
    // r is returned in another dimension order than z
    return to_matrix(r, cols(z), rows(z), 0);
  }
 /* 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);
  }
}
data {
  int<lower=1> N;  // total number of observations
  int<lower=1> N_bwt;  // number of observations
  vector[N_bwt] Y_bwt;  // response variable
  int<lower=1> K_bwt;  // number of population-level effects
  matrix[N_bwt, K_bwt] X_bwt;  // population-level design matrix
  int<lower=1> K_sigma_bwt;  // number of population-level effects
  matrix[N_bwt, K_sigma_bwt] X_sigma_bwt;  // population-level design matrix
  int<lower=1> N_tarsus;  // number of observations
  vector[N_tarsus] Y_tarsus;  // response variable
  int<lower=1> K_tarsus;  // number of population-level effects
  matrix[N_tarsus, K_tarsus] X_tarsus;  // population-level design matrix
  int<lower=1> K_sigma_tarsus;  // number of population-level effects
  matrix[N_tarsus, K_sigma_tarsus] X_sigma_tarsus;  // population-level design matrix
  int<lower=1> nresp;  // number of responses
  int nrescor;  // number of residual correlations
  // 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_bwt[N_bwt];  // grouping indicator per observation
  int<lower=1> J_1_tarsus[N_tarsus];  // grouping indicator per observation
  int<lower=1> Nby_1;  // number of by-factor levels
  int<lower=1> Jby_1[N_1];  // by-factor indicator per observation
  matrix[N_1, N_1] Lcov_1;  // cholesky factor of known covariance matrix
  // group-level predictor values
  vector[N_bwt] Z_1_bwt_1;
  vector[N_tarsus] Z_1_tarsus_2;
  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_bwt[N_bwt];  // grouping indicator per observation
  int<lower=1> J_2_tarsus[N_tarsus];  // grouping indicator per observation
  // group-level predictor values
  vector[N_bwt] Z_2_bwt_1;
  vector[N_tarsus] Z_2_tarsus_2;
  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_bwt[N_bwt];  // grouping indicator per observation
  int<lower=1> J_3_tarsus[N_tarsus];  // grouping indicator per observation
  // group-level predictor values
  vector[N_bwt] Z_3_bwt_1;
  vector[N_tarsus] Z_3_tarsus_2;
  int<lower=1> NC_3;  // number of group-level correlations
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
  int Kc_bwt = K_bwt - 1;
  matrix[N_bwt, Kc_bwt] Xc_bwt;  // centered version of X_bwt without an intercept
  vector[Kc_bwt] means_X_bwt;  // column means of X_bwt before centering
  int Kc_tarsus = K_tarsus - 1;
  matrix[N_tarsus, Kc_tarsus] Xc_tarsus;  // centered version of X_tarsus without an intercept
  vector[Kc_tarsus] means_X_tarsus;  // column means of X_tarsus before centering
  vector[nresp] Y[N];  // response array
  for (i in 2:K_bwt) {
    means_X_bwt[i - 1] = mean(X_bwt[, i]);
    Xc_bwt[, i - 1] = X_bwt[, i] - means_X_bwt[i - 1];
  }
  for (i in 2:K_tarsus) {
    means_X_tarsus[i - 1] = mean(X_tarsus[, i]);
    Xc_tarsus[, i - 1] = X_tarsus[, i] - means_X_tarsus[i - 1];
  }
  for (n in 1:N) {
    Y[n] = transpose([Y_bwt[n], Y_tarsus[n]]);
  }
}
parameters {
  vector[Kc_bwt] b_bwt;  // population-level effects
  real Intercept_bwt;  // temporary intercept for centered predictors
  vector[K_sigma_bwt] b_sigma_bwt;  // population-level effects
  vector[Kc_tarsus] b_tarsus;  // population-level effects
  real Intercept_tarsus;  // temporary intercept for centered predictors
  vector[K_sigma_tarsus] b_sigma_tarsus;  // population-level effects
  cholesky_factor_corr[nresp] Lrescor;  // parameters for multivariate linear models
  matrix<lower=0>[M_1, Nby_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[Nby_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
}
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_bwt_1;
  vector[N_1] r_1_tarsus_2;
  matrix[N_2, M_2] r_2;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_2] r_2_bwt_1;
  vector[N_2] r_2_tarsus_2;
  matrix[N_3, M_3] r_3;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_3] r_3_bwt_1;
  vector[N_3] r_3_tarsus_2;
  // compute actual group-level effects
  r_1 = scale_r_cor_by_cov(z_1, sd_1, L_1, Jby_1, Lcov_1);
  r_1_bwt_1 = r_1[, 1];
  r_1_tarsus_2 = r_1[, 2];
  // compute actual group-level effects
  r_2 = scale_r_cor(z_2, sd_2, L_2);
  r_2_bwt_1 = r_2[, 1];
  r_2_tarsus_2 = r_2[, 2];
  // compute actual group-level effects
  r_3 = scale_r_cor(z_3, sd_3, L_3);
  r_3_bwt_1 = r_3[, 1];
  r_3_tarsus_2 = r_3[, 2];
}
model {
  // likelihood including constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N_bwt] mu_bwt = Intercept_bwt + Xc_bwt * b_bwt;
    // initialize linear predictor term
    vector[N_bwt] sigma_bwt = X_sigma_bwt * b_sigma_bwt;
    // initialize linear predictor term
    vector[N_tarsus] mu_tarsus = Intercept_tarsus + Xc_tarsus * b_tarsus;
    // initialize linear predictor term
    vector[N_tarsus] sigma_tarsus = X_sigma_tarsus * b_sigma_tarsus;
    // multivariate predictor array
    vector[nresp] Mu[N];
    vector[nresp] sigma[N];
    // cholesky factor of residual covariance matrix
    matrix[nresp, nresp] LSigma[N];
    for (n in 1:N_bwt) {
      // add more terms to the linear predictor
      mu_bwt[n] += r_1_bwt_1[J_1_bwt[n]] * Z_1_bwt_1[n] + r_2_bwt_1[J_2_bwt[n]] * Z_2_bwt_1[n] + r_3_bwt_1[J_3_bwt[n]] * Z_3_bwt_1[n];
    }
    for (n in 1:N_tarsus) {
      // add more terms to the linear predictor
      mu_tarsus[n] += r_1_tarsus_2[J_1_tarsus[n]] * Z_1_tarsus_2[n] + r_2_tarsus_2[J_2_tarsus[n]] * Z_2_tarsus_2[n] + r_3_tarsus_2[J_3_tarsus[n]] * Z_3_tarsus_2[n];
    }
    for (n in 1:N_bwt) {
      // apply the inverse link function
      sigma_bwt[n] = exp(sigma_bwt[n]);
    }
    for (n in 1:N_tarsus) {
      // apply the inverse link function
      sigma_tarsus[n] = exp(sigma_tarsus[n]);
    }
    // combine univariate parameters
    for (n in 1:N) {
      Mu[n] = transpose([mu_bwt[n], mu_tarsus[n]]);
      sigma[n] = transpose([sigma_bwt[n], sigma_tarsus[n]]);
      LSigma[n] = diag_pre_multiply(sigma[n], Lrescor);
    }
    for (n in 1:N) {
      target += multi_normal_cholesky_lpdf(Y[n] | Mu[n], LSigma[n]);
    }
  }
  // priors including constants
  target += student_t_lpdf(Intercept_bwt | 3, 7.3, 2.8);
  target += student_t_lpdf(Intercept_tarsus | 3, 20.4, 5.2);
  target += lkj_corr_cholesky_lpdf(Lrescor | 1);
  target += student_t_lpdf(sd_1[1] | 3, 0, 2.8)
    - 1 * student_t_lccdf(0 | 3, 0, 2.8);
  target += student_t_lpdf(sd_1[2] | 3, 0, 5.2)
    - 1 * student_t_lccdf(0 | 3, 0, 5.2);
  target += std_normal_lpdf(to_vector(z_1));
  target += lkj_corr_cholesky_lpdf(L_1[1] | 1);
  target += lkj_corr_cholesky_lpdf(L_1[2] | 1);
  target += student_t_lpdf(sd_2[1] | 3, 0, 2.8)
    - 1 * student_t_lccdf(0 | 3, 0, 2.8);
  target += student_t_lpdf(sd_2[2] | 3, 0, 5.2)
    - 1 * student_t_lccdf(0 | 3, 0, 5.2);
  target += std_normal_lpdf(to_vector(z_2));
  target += lkj_corr_cholesky_lpdf(L_2 | 1);
  target += student_t_lpdf(sd_3[1] | 3, 0, 2.8)
    - 1 * student_t_lccdf(0 | 3, 0, 2.8);
  target += student_t_lpdf(sd_3[2] | 3, 0, 5.2)
    - 1 * student_t_lccdf(0 | 3, 0, 5.2);
  target += std_normal_lpdf(to_vector(z_3));
  target += lkj_corr_cholesky_lpdf(L_3 | 1);
}
generated quantities {
  // actual population-level intercept
  real b_bwt_Intercept = Intercept_bwt - dot_product(means_X_bwt, b_bwt);
  // actual population-level intercept
  real b_tarsus_Intercept = Intercept_tarsus - dot_product(means_X_tarsus, b_tarsus);
  // residual correlations
  corr_matrix[nresp] Rescor = multiply_lower_tri_self_transpose(Lrescor);
  vector<lower=-1,upper=1>[nrescor] rescor;
  // compute group-level correlations
  corr_matrix[M_1] Cor_1_1 = multiply_lower_tri_self_transpose(L_1[1]);
  vector<lower=-1,upper=1>[NC_1] cor_1_1;
  // compute group-level correlations
  corr_matrix[M_1] Cor_1_2 = multiply_lower_tri_self_transpose(L_1[2]);
  vector<lower=-1,upper=1>[NC_1] cor_1_2;
  // 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;
  // extract upper diagonal of correlation matrix
  for (k in 1:nresp) {
    for (j in 1:(k - 1)) {
      rescor[choose(k - 1, 2) + j] = Rescor[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_1) {
    for (j in 1:(k - 1)) {
      cor_1_1[choose(k - 1, 2) + j] = Cor_1_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_1) {
    for (j in 1:(k - 1)) {
      cor_1_2[choose(k - 1, 2) + j] = Cor_1_2[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];
    }
  }
}