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