data { int N; int N_groups; int groups[N]; vector[N] y; } transformed data { matrix [N_groups,N_groups] A; matrix [N_groups,N_groups] Q; matrix [N_groups,N_groups] R; matrix [N_groups,N_groups] R_ast_inverse; A[2:N_groups, 1:N_groups-1] = diag_matrix(rep_vector(1, N_groups - 1)); A[1, 1:N_groups] = rep_row_vector(-1, N_groups); A[, N_groups] = rep_vector(1, N_groups); Q = qr_thin_Q(A) * sqrt(N_groups - 1); R = qr_thin_R(A) / sqrt(N_groups - 1); R_ast_inverse = inverse(R); } parameters { vector[N_groups] group_z; real sigma; real group_sd; real group_mean; //real intercept; } transformed parameters { vector[N_groups] theta = Q * group_z; vector[N_groups] R_group_z = R_ast_inverse * group_z; } model { R_group_z[1:N_groups] ~ normal(group_mean, group_sd); group_sd ~ lognormal(0, 0.3); sigma ~ cauchy(0, 1); //intercept ~ normal(3, 1); //y ~ normal(theta[groups] + intercept, sigma); y ~ normal(theta[groups], sigma); } generated quantities { vector[N_groups] group_r_hat = A[, 1:N_groups] * R_group_z[1:N_groups]; }