Understanding the Mathematical Model

I’d like to have a clear understanding of the mathematical model behind a formula and priors. Take for example the following call:

fit1 <- brm(brmsformula(time | cens(censored) ~ (1 | 0 | group), shape ~ (1 | 0 | group)),
data = mydata, family = weibull(),
prior = c(set_prior("normal(0,5)", class = "Intercept"),
          set_prior("normal(0,4)", dpar="shape", class = "Intercept"),
          set_prior("normal(0,0.1)", class = "sd", group="group")),
warmup = 1000, iter = 2000, chains = 4,
control = list(adapt_delta = 0.95))

What is the mathematical model here?

time_g = a + a_g
time_g \sim Weibull(\beta_g, \eta) and right-censored
a \sim N(0,5)
a_g \sim N(0,0.1)
\beta_g = b
b\sim N(0,4)

Is this correct?

Also, after fitting the model and calling summary(), I see cor and sd. How do those relate to the mathematical model?

The mathematical model looks somewhat different than what you suggested. The model can be written down as follows:

time_n \sim Weibull(\beta_n, \eta_n) \text{ and right censored} \\ \eta_n = a + a_{g_n} \\ \beta_n = b + b_{g_n} \\ a \sim N(0, 5) \\ b \sim N(0, 4) \\ (a_g, b_g) \sim MN(0, \Sigma) \\ \Sigma = \left( \begin{matrix} \sigma_a^2 & \sigma_a \sigma_b \rho \\ \sigma_a \sigma_b \rho & \sigma_b^2 \end{matrix} \right) \\ \sigma_a \sim N(0, 0.1) \\ \sigma_b \sim N(0, 0.1) \\ \rho \sim LKJ(1) \text{ (default prior for correlation matrices)}

I also suggest you take a look at stancode(fit1).

2 Likes

Thanks Paul! I edited my question two times or so, not sure if you noticed it. Are both \beta and \eta group-level parameters? I only define a formula for the shape in addition to time.

Please look at stancode(fit1) to see the details. I also suggest to take a look at the parameterization of the weibull family in vignette("brms_families").

Thanks! There are still a few things I don’t understand about the Stan model description I get with stancode. For example, there’s this line:

target += normal_lpdf(to_vector(z_1) | 0, 1);

where z_1 appears in other lines. Under parameters:

matrix[M_1, N_1] z_1; // unscaled group-level effects

and transformed parameters:

matrix[N_1, M_1] r_1 = (diag_pre_multiply(sd_1, L_1) * z_1)';

I don’t know what that z_1 parameter corresponds to and why it’s \sim N(0,1) seemingly by default.

In your model, r_1 contains the a_g and b_g parameters. The parameterization via z_1 is called “non-centered” paramerization, which is used because it often leads to better convergence than setting priors on the varying effects r_1 directly.

Paul, the correlation matrix you wrote seems to relate to a_g and b_g of the same group g. Is this what’s called “intra-class correlation”? Does it make sense to model the correlation between parameters of different groups (inter-class)? If so, how can it be coded with brms?

No, the ICC is something rather different. I believe that if you google it, you will find some good resources that explain ICCs.

This question is too general to be reasonably answered. If you don’t have a clear reasoning for that, I would probably not bother with it.

Paul, how do I specify a prior on the Cholesky correlation factor in this example? I tried adding a prior to the list of priors,

set_prior("lkj(2)", class = "cor")

But I get an error:

Error: Duplicated prior specifications are not allowed.

Apparently you try to specify a prior on the correlations twice ;-)

I figured out the problem. I was reloading a saved model and then trying to set a prior on that parameter, but the reloaded model already has the default prior set with lkj(1). But I only get that error when I have the cor prior in the list of priors. If I try to update the model without the cor prior, I don’t the get the error.