Error message with user defined priors in brms

Modeling
1

I am new to brms and have been getting a warning message when I specify priors.

The model is a re-analysis of a randomised trial published in New England J Medicine. It is a simple logistic regression model (dead ~ treat) comparing ECMO (a form of artificial lung support for patients with pneumonia) vs control on mortality outcome in patients in ICU. Based on neutral beta priors (converted to normal priors) I specified the priors as:

prior_neutral <- c(
  set_prior("normal(0, 0.33132171)",class = "b", coef = ""),
  set_prior("normal(0, 0.33132171)", class = "b", coef = "treatECMO"),
  set_prior("student_t(3, 0, 2.5)",class = "Intercept", coef = "") # use default for this prior
)

The model is:

                     bernoulli(link = "logit"), 
                     data = EOLIA,
                     prior = prior_neutral,
                     chains = 2, # nb of chains
                     iter = 5000, # nb of iterations, including burnin
                     warmup = 1000, # burnin
                     thin = 1)

When I run the model, I get the following warning:
Warning message:
The global prior ‘normal(0, 0.33132171)’ of class ‘b’ will not be used in the model as all related coefficients have individual priors already. If you did not set those priors yourself, then maybe brms has assigned default priors. See ?set_prior and ?get_prior for more details.

However, when I run the following code:

I get:
prior_summary(EOLIA_neutral)
prior class coef group resp dpar nlpar lb ub source
normal(0, 0.33132171) b user
normal(0, 0.33132171) b treatECMO user
student_t(3, 0, 2.5) Intercept user

This seems to imply the model is using my defined priors.

So, (finally) my question is: what does the warning message mean and are my priors being used by the model (which seems to be the case.

Any advice gratefully received.

Many thanks,
David

2

Realised I have not included the entire model. It should be


EOLIA_neutral <- brm(dead ~ treat, 
                     bernoulli(link = "logit"), 
                     data = EOLIA,
                     prior = prior_neutral,
                     chains = 2, # nb of chains
                     iter = 5000, # nb of iterations, including burnin
                     warmup = 1000, # burnin
                     thin = 1)

Thanks,
David

3

Try this:

prior_neutral <- c(
  set_prior("student_t(3, 0, 2.5)", class = "Intercept"),
  set_prior("normal(0, 0.33132171)", class = "b", coef = "treatECMO")
)
4

Hi Solomon,

Thanks very much. No error message with your prior. The prior summary for your prior is:

                 prior     class      coef group resp dpar nlpar lb ub  source
                (flat)         b                                       default
 normal(0, 0.33132171)         b treatECMO                                user
  student_t(3, 0, 2.5) Intercept                                          user

The prior summary for my prior is:

                 prior     class      coef group resp dpar nlpar lb ub source
 normal(0, 0.33132171)         b                                         user
 normal(0, 0.33132171)         b treatECMO                               user
  student_t(3, 0, 2.5) Intercept                                         user

So it looks like your prior for the beta coefficient is flat, while mine is N((0, 0.33132171). But the coefficient for treatECMO is the same from the two models. Which I don’t understand.

What I want to do is specify a prior for the control group and the treatment group separately. I have assumed that by setting the prior on the beta coefficient I am setting my prior for the control group and by setting my prior for the beta coefficient for treatECMO, I am setting my prior for the intervention group.
In the model I have posted, the two priors are the same (ie, neutral) but I have models where I have specified different priors for beta (control) and beta/treatECMO (treatment group).

Would appreciate a steer. I am misinterpreting my prior on the beta (control group)?

Thanks,
David

5

Hey @dsidebotham, let’s work through this slowly.

Your model formula is dead ~ treat, which means you only have one predictor variable, treat. Based on the output you’ve shown, treat is a 2-level factor or character variable for which the second level is ECMO. Based on your initial post, the reference category is for the control condition.

When you set priors with brms, the priors for the non-intercept \beta coefficients are of class = b, which you can see from your prior summary (second column). brms allows users to set a generic prior for all \beta coefficients and/or priors for specific predictors. If you want to set a generic prior for all \beta coefficients, leave the coef argument empty–and by “empty” I mean don’t type that argument in at all. You can see the coef information in the third column of your prior summary. If you’d like to set a prior for a specific \beta coefficient, such as with your treat predictor, you can also use the coef argument. But as you can see from the prior summary output, brms has the additional complication that when your variable is a character or factor, it will tack on the name of the non-reference category as a suffix. So in your case, your treat variable got the suffix ECMO. This is why I recommended the syntax set_prior("normal(0, 0.33132171)", class = "b", coef = "treatECMO").

Also note that in the case of your model with only one non-intercept \beta coefficient, the predictor-specific prior is identical to the global \beta prior. So when you’re writing your prior syntax, choose only one. That is, either do this

prior_neutral <- c(
  set_prior("student_t(3, 0, 2.5)", class = "Intercept"),
  set_prior("normal(0, 0.33132171)", class = "b")
)

or do this

prior_neutral <- c(
  set_prior("student_t(3, 0, 2.5)", class = "Intercept"),
  set_prior("normal(0, 0.33132171)", class = "b", coef = "treatECMO")
)

but do NOT do this

prior_neutral <- c(
  set_prior("student_t(3, 0, 2.5)", class = "Intercept"),
  set_prior("normal(0, 0.33132171)", class = "b", coef = ""),
  set_prior("normal(0, 0.33132171)", class = "b", coef = "treatECMO")
)

You should also know that when you use a coefficient-specific prior with the coef argument, that overrides any settings of the global prior. This might not make sense in the context of a single-predictor model like yours, but it’s very handy with a model with many predictors. For example, say I have the model

\begin{align} y_i & \sim \mathcal N (\mu_i, \sigma) \\ \mu_i & = \beta_0 + {\color{blue}{\beta_1 x_{1i}}} + \beta_2 x_{2i} + \beta_3 x_{3i} \end{align}

where {\color{blue}{x_1}} is my focal predictor, which I’ve tried to indicate with blue font, and x_2 and x_3 are some covariates I need to improve the model, but don’t really care about from an interpretive perspective. Let’s also, for the sake of simplicity, assume all variables are standardized continuous variables. In such a case, I might want to use a theory-specific prior for {\color{blue}{x_1}}, but just some generic weakly-regularizing prior for x_2 and x_3. Here’s what that prior code might look like:

prior_neutral <- c(
  set_prior("student_t(3, 0, 2.5)", class = "Intercept"),
  set_prior("normal(0, 1)", class = "b"),
  set_prior("normal(0, 0.33132171)", class = "b", coef = "x1")
)

With this syntax, my focal variable {\color{blue}{x_1}} got its own theory-based prior, \mathcal N (0, 0.33132171), but my two less-interesting covariates x_2 and x_3 both got the generic weakly-regularizing prior \mathcal N (0, 1).

1 Like
6

Hi Solomon,

Thanks for taking the time for such a detailed reply. Greatly appreciated and very helpful.
I think I see my mistake. For the treatECMO prior I need the effect size - not the prior for the event rate in the ECMO group.

I have four models using different priors (flat, neutral, enthusiastic, very enthusiastic). I calculated the normal priors from beta priors using a logit transformation. I am now just setting the treatECMO prior using priors for the effect size and I get no error messages. Fantastic!

Interestingly, the credible intervals for the four models agree almost completely with the intervals obtained from a beta-binomial model.
As you may have guessed, I am using these examples to gain confidence in brms.

Once again, thanks very much for your detailed explanation.

David

1 Like