data{ int P; // total number of periods int N; // total number of doctor int M; // total number of data points int K; // number of demographic variables (with intercept) int R; //number of individual-level coefficients (with intercept) int y[M]; // outcome vector (TRx) int t[M]; // period vector int d[M]; // details vector int id[M]; // doctor id vector real l[M]; //log(lagged TRx+1) } transformed data{ matrix[M, R] X = append_col(rep_vector(1,M),append_col(to_vector(d),to_vector(l))); } parameters{ real phi; //dispersion parameter for NBD matrix[R, N] eta; // normal errors for individual-level coefficients cholesky_factor_corr[R] L_Omega; // cholesky corr. for individual-level coefficients row_vector[R] mu_beta; // mean for individual-level coefficients vector[R] tau_unif; // scaled variance of individual-level coefficients //parameters to determine non-randomness real gam0; real gam1; real gam2; } transformed parameters{ matrix[N,R] zeta; // individual coefficients vector[R] tau; // variance of individual coefficients vector[N] nu; for (r in 1:R) tau[r] = 2.5 * tan(tau_unif[r]); zeta = rep_matrix(mu_beta,N) + (diag_pre_multiply(tau,L_Omega) * eta)'; for (i in 1:N) nu[i] = gam0+gam1*(zeta[i,1]/(1-zeta[i,3]))+gam2*(zeta[i,2]/(1-zeta[i,3])); } model{ vector[M] theta=rep_vector(0,M); phi ~ gamma(2,1); L_Omega ~ lkj_corr_cholesky(2); to_vector(eta)~std_normal(); to_vector(mu_beta) ~ normal(0,5); // NBD likelihood for (r in 1:R) theta += zeta[id,r].*X[,r]; // target += neg_binomial_2_log_lpmf(y|theta,rep_vector(phi,M)); target += poisson_log_lpmf(d|nu[id]); } generated quantities { matrix[R,R] Omega; matrix[R,R] Sigma; Omega = multiply_lower_tri_self_transpose(L_Omega); Sigma = quad_form_diag(Omega, tau); }