Help with priors

Interfaces brms
1

Dear everyone,

I am very new to Bayesian stats and even newer to Stan/brms.
So i guess my question is basic; I apologize for that already…!

I am first trying to re-create an analysis that was performed by a colleague in Python using pyMC3. So the question is not whether the priors and model make sense, but about how I can translate it to Stan/brms.

I have the following data:

my_data <- structure(list(Material = c("Flint", "Flint", "Flint", "Flint", "Flint", "Flint", "Flint", "Quartzite", "Quartzite", "Quartzite", "Quartzite", "Quartzite", "Quartzite", "Quartzite", "Quartzite"), Brush_Dirt = c("No_Brush_No_Dirt", "No_Brush_Yes_Dirt", "No_Brush_Yes_Dirt", "Yes_Brush_No_Dirt", "Yes_Brush_No_Dirt", "Yes_Brush_Yes_Dirt", "Yes_Brush_Yes_Dirt", "No_Brush_No_Dirt", "No_Brush_No_Dirt", "No_Brush_Yes_Dirt", "No_Brush_Yes_Dirt", "Yes_Brush_No_Dirt", "Yes_Brush_No_Dirt", "Yes_Brush_Yes_Dirt", "Yes_Brush_Yes_Dirt"), Z = c(0.244351382982314, -0.0015361217248987, -0.877775019539761, -0.498374272882404, -0.659639015247898, 0.488286794228672, -0.794444097789949, 0.909867363652624, 2.05669202585517, 0.019788897546955, 1.31082564239473, -0.0230158910253944, -0.62363353368801, 0.498543577208916, -2.04993773197107)), row.names = c(NA, 15L), class = "data.frame")

The model is described in the accompanying graph (as PDF).

For now, I have this unfinished model (“work in progress”):

fit <- brm(data = my_data, family = gaussian,
           Z ~ 1 + Material * Brush_Dirt,
           prior = c(prior(normal(0, sigma0), class = Intercept),                    # beta0
                     prior(normal(0, sigma1), class = b, lb = 0, coef = Material),   # beta1
                     prior(normal(0, sigma2), class = b, lb = 0, coef = Brush_Dirt), # beta2
                     prior(uniform(0, ErrorMax), class = sigma))                     # sigma

I have three questions:

  1. How do I define the priors for the interaction (M)?

  2. How do I define the priors for sigma0, sigma1, sigma2 and sigmaM [~HN(0, sX)]? Do I need several formulas with the bf() function (one for each sigmaX) and specify the dpar argument (= sigmaX)?
    But then, what would be the class for the prior() statement?

I have also realized that priors do not necessarily need to be specified in the model. I guess default values will be used.
3) How can I find out what priors will be used for a given model? Is it what the function get_prior() is for? Do I run it before or after brm()?

I have tried to search online for that, but I found it very difficult to know what and where to search since I am not familiar with many of the technical terms.

I would be very grateful if you could help me with that, and also point me to some useful resources.

Thank you in advance!

Bayesian-model_brushing.pdf (33.3 KB)

2
  1. Use get_prior() to find the coefficients names, then assign priors to those.
model <- bf(Z ~ 1 + Material * Brush_Dirt)
get_prior(model, my_data)
# For example
prior_int1 <- prior(normal(0, 1), class = "b", coef = "MaterialQuartzite:Brush_DirtNo_Brush_Yes_Dirt")

  1. You can’t with brms. You need to fix the hyperparameters so e.g. prior(normal(0, 1), class = Intercept) instead of using sigma0 as the sd prior. (edit: there may be a way around this now within brms, but you can do it by modifying the Stan code output by make_stancode() to achieve this.)

  2. You can use get_prior() as shown above.

3

Thank you for the quick and very helpful answer!

  1. At my level, it would make it even more difficult if I have to switch back and forth between brms and Stan. I mean, I would have to learn another language to make it work.
    So I hope there is indeed a way to specify priors for priors within brms, and that someone from this forum will be able to help me with that :)
4

Short question related to customs on this forum: should I “like” your post (it does help me) and/or should I mark that it solves the problem (even though it solves it partially)?

5
  1. I thought I had understood your answer, but apparently not. I have checked the get_prior() output but there are things I don’t get:
  • Why only 1 material (quartzite) shows up…? Is “MaterialFlint” missing in the column coef for rows 1-4?
  • Where are the rows for the main effect (“Material” and “Brush_Dirt”), i.e. beta1 and beta2?
  • I would like to specify the same normal priors for all interaction terms. Can I do that in one go? I would like to avoid specifying 8 priors…
6

I don’t think anyone will be able to help you put priors on priors with brms, because you can’t do that with brms. However, it is relatively straightforward to 1) create your brms model without said priors, 2) edit the Stan code of said model 3) put edited Stan code back to the brm model, 4) sample, 5) profit. (Note I am not sure if that model makes sense.) See ?brm. Learning Stan is tremendously helpful but I understand it takes time that you might not have. In that case you may consider hiring a specialist if no one here is willing to volunteer their time.

You can mark an answer whenever you feel like it appropriately answers your question.

I thought I had understood your answer, but apparently not. I have checked the get_prior() output but there are things I don’t get:

  • Why only 1 material (quartzite) shows up…? Is “MaterialFlint” missing in the column coef for rows 1-4?

Flint is used as the baseline because by default categorical variables are dummy coded (and Flint is alphabetically before Quartzite).

  • Where are the rows for the main effect (“Material” and “Brush_Dirt”), i.e. beta1 and beta2?

Row 2: Brush_DirtNo_Brush_Yes_Dirt; row 5: MaterialQuartzite.

  • I would like to specify the same normal priors for all interaction terms. Can I do that in one go? I would like to avoid specifying 8 priors…

(There are 3 interaction terms in your model.) You can’t specify a prior for all interaction terms with a one-liner. However, you can take advantage of the fact that get_prior() returns a data frame, so something like this:

priors <- get_prior(model, my_data)
priors <- dplyr::mutate(
  priors, 
  prior = ifelse(str_detect(coef, ":"), "normal(0, 1)", prior)
)
priors[,1:3]

                  prior     class                                           coef
1                               b                                               
2                               b                    Brush_DirtNo_Brush_Yes_Dirt
3                               b                    Brush_DirtYes_Brush_No_Dirt
4                               b                   Brush_DirtYes_Brush_Yes_Dirt
5                               b                              MaterialQuartzite
6          normal(0, 1)         b  MaterialQuartzite:Brush_DirtNo_Brush_Yes_Dirt
7          normal(0, 1)         b  MaterialQuartzite:Brush_DirtYes_Brush_No_Dirt
8          normal(0, 1)         b MaterialQuartzite:Brush_DirtYes_Brush_Yes_Dirt
9  student_t(3, 0, 2.5) Intercept                                               
10 student_t(3, 0, 2.5)     sigma
3 Likes
7

I was off for 2 weeks…
Thank you for the detailed answer! I will look into it and I guess I’ll have to learn Stan indeed.

8

It looks like you could achieve question 2. in your original post by using the set_prior() and stanvar() functions: https://github.com/paul-buerkner/brms/issues/459#issuecomment-399041581

1 Like
9

It looks very promising! Thank you for the pointer!