Hi there,

I’m trying to formulate the mathematical model of a multivariate multilevel model with multiple skew-normally distributed predictors for an upcoming paper and I have trouble turning the Stan code into a proper mathematical model, especially the parts concerning r_1, J_1 and Z_1. I’m also confused about the parameters of the skew-normal distribution, do I report the transformed parameters (i.e. delta), the original ones (i.e. alpha) or both? If it’s the latter, which I suppose it is, is the transformation part of the mathematical model itself?

```
model <- brm(mvbind(outcome1, outcome2) ~ pred1+ pred2+ pred3+ pred4+ pred5+ pred6+ pred7+ pred8+ (pred1+ pred2+ pred3+ pred4+ pred5+ pred6+ pred7+ pred8| subject),
data = df, family = skew_normal(),
iter = 2000, warmup = 1000, chains = 4, cores = 4)
```

The model I have thus far (which isn’t much) is:

outcome1 ~ SN(µ, ω, α)

outcome2 ~ SN(µ, ω, α)

µ is where the problems start, since that’s when r_1 and Z_1 come into play. This is what it looks like in the Stan output:

mu_outcome1[n] += r_1_outcome1_1[J_1_outcome1[n]] * Z_1_outcome1_1[n] + r_1_outcome1_2[J_1_outcome1[n]] * Z_1_outcome1_2[n] + r_1_outcome1_3[J_1_outcome1[n]] * Z_1_outcome1_3[n] + r_1_outcome1_4[J_1_outcome1[n]] * Z_1_outcome1_4[n] + r_1_outcome1_5[J_1_outcome1[n]] * Z_1_outcome1_5[n] + r_1_outcome1_6[J_1_outcome1[n]] * Z_1_outcome1_6[n] + r_1_outcome1_7[J_1_outcome1[n]] * Z_1_outcome1_7[n] + r_1_outcome1_8[J_1_outcome1[n]] * Z_1_outcome1_8[n] + r_1_outcome1_9[J_1_outcome1[n]] * Z_1_outcome1_9[n];

I could also provide the rest of the Stan output, which I omitted from this post due to its length :)

Thank you in advance!

P.S.: I wonder why there isn’t a function that turns Stan code directly into a mathematical model. Maybe I’m overlooking something, but to me, it would seem generally feasible ;) If someone did this, they’d be my personal hero