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);
?