Skew normal superpopulation for multilevel modeling

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?