Setting lower bound for intercept parameter in brms?

Interfaces brms
1

With set_prior, brms allows to set lower bounds for regression weights b but not for the Intercept parameter.

Is there a different way to set a lower bound for a parameter?

One use case is regression models with a likelihood function that is constrained to positive values (e.g. Gamma) in combination with an identity link (default e.g. for Gamma).

PS:
I am aware of following workaround:

  1. Fit an initial model without bounds
  2. Modify model code to implement bound
  3. Generate new stanfit-object with updated model code
  4. Replace stanfit object in original brmsfit object with the new stanfit object
  5. Update the now modified brmsfit model to obtain samples from the updated model (with compile = F).

This is what I am currently doing. It works, but it is kind of hacky.

1 Like
2

Here is a simpler way to use a custom Stan model to implement lower bounds (using rstan as the backend):

  1. Create brmsfit object without stanfit object by using the empty option
    bfit = brms(my_formula, my_data, my_family, empty = T)

  2. Update the Stan model, e.g.
    bfit$model = gsub("real Intercept;","real<lower=0> Intercept;",bfit$model)

  3. Fit the model
    my_standata = make_standata(my_formula, my_data, my_family)
    bfit$fit = stan(model_code = bfit$model, data = my_standata

ht @paul.buerkner

1 Like
3

FWIW I think it would probably also work with cmdstanr as the backend if you change step 3 to something like this:

mod = cmdstan_model(stan_file = write_stan_tempfile(bfit$model))
fit = mod$sample(data = my_standata)
bfit$fit = rstan::read_stan_csv(fit$output_files())
1 Like
4

Agreed, I didn’t want to imply to imply that it doesn’t work with cmdstanr. I just didn’t have the time to write it down then and wanted to avoid that somebody tries using the code without being aware that this example works with ratan.

1 Like
5

Another solution is to use the formula = y ~ 0 + intercept + … syntax where the intercept parameter is of class = b.

1 Like
6

I thought of this solution (should have mentioned it) and rejected it for the following reasons:

  • I want to give the intercept a lower bound, but not the other regression weights
  • There are often good reasons that the prior for the intercept is wider than the prior for regression weights

These two points are of course only correct if one sets only one prior for all regression weights. I never checked if brms allows to set different priors for different regression weights like this:

my_prior = 
  set_prior(normal(0,10), class = "b", coef = "my_intercept", lb = 0) + 
  set_prior(normal(0,2), class = "b", coef = "x1")

If this is possible, the solution to with y ~ 0 + ... is clearly more elegant.

Edit: OK, I checked, and @Solomon’s solution does not work (except one thinks it’s OK giving intercept and regression weights the same prior and boundaries), because

Argument ‘coef’ may not be specified when using boundaries.

What does work is to use different priors, and not to specify boundaries i.e. :

my_prior = 
  set_prior(normal(0,10), class = "b", coef = "my_intercept") + 
  set_prior(normal(0,2), class = "b", coef = "x1")

However, this means that Stan can use negative values for my_intercept during initialisation, which is problematic for the data set I was am currently working with. More generally, it is my understanding that it is good practice that the one should add the correct boundaries when specifying model parameters.

In sum, the solution @Solomon proposes gets one half the way. If one wants correct specification of prior distribution and parameter boundaries, as of now one has to work with modified stan_code as described above.

2 Likes
7

Interesting. I was not aware that Argument ‘coef’ may not be specified when using boundaries. Good to know.

8

I am thinking there has to be a reason this limitation exists. Is it simply because it has not been implemented yet, or is it because it would make things much more complicated for the sampler?

9

It’s not a problem for the sampler.

10

It is simply not implemented yet. The bound features of set_prior and not really far developed and will get a complete overhaul at some point to support boundaries more generally rather than just for a few specific parameter classes.

2 Likes
12

would love to see the lb = option implemented for other classes including sd

1 Like