Trouble with Mathematical Model

Interfaces brms
1

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

1 Like
2

Sorry, your question fell through despite being well written and relevant.

The complicated code is just a way to write Z matrix from the matrix notation of a mixed model (Mixed model - Wikipedia). Sometimes, just stating the model in matrix notation is sufficient for the paper. If you want to be more detailed, you can decompose the matrix in individual coefficients. An example of how that might look like is on page 8 of my recent publication

Definitely also check out
Understanding the Mathematical Model and Mathematical description of a brms model where a similar problem is discussed.

Those question arise from time to time. The problem is that this is easy only for some models. But Stan (and even just brms) are very general and doing that generally is hard (and would often produce useless output anyway as it would just rewrite the code in greek which would be hard to understand). I’ve definitely thought about adding that capability to at least a common subset of brms models, but upon inspection it turned out to be likely quite a demanding project (but IMHO suitable for a motivated undergraduate student).

I would also mention that a growing part of academic readership would IMHO be faster to understand your model from the R formula than from the math notation, so - depending on your target audience - doing away with the math notation completely might be feasible, though a bit riskier as it is not (yet) commonly done.

1 Like
3

Thank you very much for your informative answer! I think I will either go with the matrix notation or do away with the mathematical model entirely since the target audience probably won’t be all that familiar with Bayesian statistics in general, at least from what I can gather. Some rather basic explanations are therefore required anyway :)

Thanks again!