I am working to become proficient in Bayesian methods. I’ve been applying them to my work using the brms package. My interest is this: I have fit a multilevel model to some data which results in a superpopulation for intercepts and slopes. Hence normal distributions describing group to group variability in the intercept and slope. My data exhibits left skewed behavior in the observed group intercepts and slopes. How can I specify skewed normal distributions for this superpopulation model on the intercepts and slopes?

I’ve pulled some of the stan code from the brm fit that I believe will address this need. It is below.

model {

// likelihood including constants

if (!prior_only) {

// initialize linear predictor term

vector[N] mu = Intercept + rep_vector(0.0, N);

for (n in 1:N) {

// add more terms to the linear predictor

mu[n] += r_1_1[J_1[n]] * Z_1_1[n] + r_1_2[J_1[n]] * Z_1_2[n];

}

target += normal_id_glm_lpdf(Y | Xc, mu, b, sigma);

}

// priors including constants

target += student_t_lpdf(Intercept | 3, 20.4, 2.5);

target += student_t_lpdf(sigma | 3, 0, 2.5)

- 1 * student_t_lccdf(0 | 3, 0, 2.5);

target += student_t_lpdf(sd_1 | 3, 0, 2.5)

- 2 * student_t_lccdf(0 | 3, 0, 2.5);

target += std_normal_lpdf(to_vector(z_1));

target += lkj_corr_cholesky_lpdf(L_1 | 1);

}

I’m not perfectly sure but I believe I need to alter this code for the target distribution on “Intercept” and “z_1” and am wondering how I can do this to use ** skew_normal_lpdf** in place of the distributions that are defaults. Can this be done? How do I do this?