data{ int Nsubj ; // number of participants int Ntotal ; // number of questions int nYlevels ; // number of responses on scale (0-10 scale, so 11 responses) int s[Ntotal] ; // index of subject at a given question int k ; // num. of predictors int yPreA[Ntotal] ; int yPreFb[Ntotal] ; int yJud[Ntotal] ; int yPost[Ntotal] ; int yPreB[Ntotal] ; int yRec[Ntotal] ; } parameters { real mu_tau ; real sigma_tau ; real omega_corr_pos ; vector[k] tau[Nsubj] ; cholesky_factor_corr[k] Omega_pos[Nsubj] ; vector[Nsubj] muPreASub ; real muPreAGroup ; real sigma2PreAGroup ; vector[Nsubj] muPreFbSub ; real muPreFbGroup ; real sigma2PreFbGroup ; vector[Nsubj] muJudPosSub ; real muJudPosGroup ; real sigma2JudPosGroup ; vector[Nsubj] muPostSub ; real muPostGroup ; real sigma2PostGroup ; vector[Nsubj] muPreBSub ; real muPreBGroup ; real sigma2PreBGroup ; vector[Nsubj] muJudRecSub ; real muJudRecGroup ; real sigma2JudRecGroup ; } transformed parameters{ matrix[k, k] quad_pos[Nsubj]; for (i in 1:Nsubj) { quad_pos[i,,] = diag_pre_multiply(tau[i,], Omega_pos[i,,]); } } model { mu_tau ~ cauchy(0, 2.5); sigma_tau ~ uniform(0.001, nYlevels); omega_corr_pos ~ uniform(0.001, nYlevels); for (i in 1:Nsubj) { tau[i,] ~ gamma((mu_tau/sigma_tau)^2, mu_tau/(sigma_tau^2)); Omega_pos[i,,] ~ lkj_corr_cholesky(omega_corr_pos); } sigma2PreAGroup ~ uniform(0.001, nYlevels); muPreAGroup ~ uniform(1 , nYlevels); muPreASub ~ normal(muPreAGroup, sqrt(sigma2PreAGroup)) ; sigma2PreFbGroup ~ uniform(0.001, nYlevels); muPreFbGroup ~ uniform(1 , nYlevels); muPreFbSub ~ normal(muPreFbGroup, sqrt(sigma2PreFbGroup)) ; sigma2JudPosGroup ~ uniform(0.001, nYlevels); muJudPosGroup ~ uniform(1 , nYlevels); muJudPosSub ~ normal(muJudPosGroup, sqrt(sigma2JudPosGroup)) ; sigma2PostGroup ~ uniform(0.001, nYlevels); muPostGroup ~ uniform(1 , nYlevels); muPostSub ~ normal(muPostGroup, sqrt(sigma2PostGroup)) ; sigma2PreBGroup ~ uniform(0.001, nYlevels); muPreBGroup ~ uniform(1 , nYlevels); muPreBSub ~ normal(muPreBGroup, sqrt(sigma2PreBGroup)) ; sigma2JudRecGroup ~ uniform(0.001, nYlevels); muJudRecGroup ~ uniform(1 , nYlevels); muJudRecSub ~ normal(muJudRecGroup, sqrt(sigma2JudRecGroup)) ; for (i in 1:Ntotal) { [yPreA[i], yPreFb[i], yJud[i], yPost[i], yPreB[i], yRec[i]] ~ multi_normal_cholesky([muPreASub[s[i]], muPreFbSub[s[i]], muJudPosSub[s[i]], muPostSub[s[i]], muPreBSub[s[i]], muJudRecSub[s[i]]], quad_pos[s[i],,]); } } generated quantities { vector[Ntotal] log_lik; real yPost_new[Ntotal]; real yPreB_new[Ntotal]; real yRec_new[Ntotal]; vector[k] y_new[Ntotal]; vector[k] D_pos[Nsubj]; corr_matrix[k] Omega_pos_true[Nsubj]; y_new[,1] = yPreA; y_new[,2] = yPreFb; y_new[,3] = yJud; y_new[,4] = yPost_new; y_new[,5] = yPreB_new; y_new[,6] = yRec_new; // convert cov mtx to corr mtx: for (i in 1:Nsubj) { D_pos[i] = sqrt(diagonal(multiply_lower_tri_self_transpose(quad_pos[i,,]))); Omega_pos_true[i] = inverse(diag_matrix(D_pos[i])) * multiply_lower_tri_self_transpose(quad_pos[i,,]) * inverse(diag_matrix(D_pos[i]))'; } // log-likelihood for PSIS-LOO: for (n in 1:Ntotal) { log_lik[n] = multi_normal_cholesky_lpdf([yPreA[n], yPreFb[n], yJud[n], yPost[n], yPreB[n], yRec[n]] | [muPreASub[s[n]], muPreFbSub[s[n]], muJudPosSub[s[n]], muPostSub[s[n]], muPreBSub[s[n]], muJudRecSub[s[n]]], quad_pos[s[n],,]); } // generate simulated data: for (i in 1: Ntotal) { y_new[i,1:k] = multi_normal_cholesky_rng([muPreASub[s[i]], muPreFbSub[s[i]], muJudPosSub[s[i]], muPostSub[s[i]], muPreBSub[s[i]], muJudRecSub[s[i]]], quad_pos[s[i],,]); } }