Brms - constrain group-level effects to positive values?

[jonnations]

Hi,

I am estimating a parameter for n-groups using partial pooling. The parameter in my case represents a distance from zero (centroid) in an n-dimensional volume. My sampling for each group is uneven, and the partial pooling approach using group level effects works very well! The problem is that values below zero do not make sense for this parameter, as these groups can not have a negative distance from a centroid. I would like to somehow limit my posteriors to positive values. Below is an example of what I mean, which I took from the Tidyverse website and modified slightly.

pacman::p_load(tidyverse, tidybayes, brms)

set.seed(11)
n = 10
ABC =
  tibble(
    condition = rep(c("A", "B", "C", "D", "E", "F", "G", "H"), n),
    response = rnorm(n * 8, 
                     c(0.3, 1, 2, 1.5, 2.5, 3.5, 2, 1.75), 
                     c(0.17, 0.3, 0.5, 0.5, 0.8, 1, 0.5, 0.7)))

MODEL

m = brm(
  response ~ 0 + (1|condition),
  data = ABC,
  refresh = 0,
  prior = c(
    prior(student_t(3, 0, 1), class = sd),
    prior(student_t(3, 0, 1), class = sigma)),
  control = list(adapt_delta = .99))

PLOT

m %>%
  spread_draws(r_condition[condition,]) %>%
  median_qi(r_condition, .width = c(.95, .66)) %>% 
  ggplot(aes(y = condition, x = r_condition, xmin = .lower, xmax = .upper)) +
  geom_pointinterval()

As you can see, the posterior distribution for group A goes below zero. I understand that upper and lower bounds can be fixed on population-level effects in brms. But this doesn’t seem directly possible with group-level effects, or at least I can’t figure out how.

I came up with a few ideas, but I’m not sure any of them are appropriate. I could use a lognormal family, but I’d rather not deal with the transformations if possible. The second idea is using the exponential gaussian family - exgaussian(). This did work, and the posteriors are all positive, but in my empirical data the results are very left-skewed and the median estimates increased (images below). Also, I’m not sure that this is the appropriate use of this distributional family. A brute force approach would be to take the absolute values of the posteriors (as a negative distance is still a distance?) but I don’t feel really comfortable with that.

Can anyone suggest a good approach to dealing with this problem?

Thanks!

• Operating System: OSX Catalina
• brms Version: 2.13.0

Plots from my empirical examples:

family = gaussian()

family = exgaussian()
Values are positive! But notice left skew and increased median estimates.

[Solomon]

One approach would be to use a lognormal prior. McElreath used this approach a few times in the second edition of his text. For example, here’s a link to my brms translation of his example from Section 4.4.1. Here’s a link to a point in one of his recent lectures where he discusses this example.

[jonnations]

Thanks Solomon! I had tried using different distribution families, including lognormal, but somehow it didn’t occur to me to use a prior that restricts the parameter to positive values.

I’m glad to see you’re working through Rethinking ed. 2! Thanks for the help.

[jonnations]

Ok, maybe I wasn’t as confused as I thought (or maybe I’m more confused). My hope was to restrict the group-level estimates to positive values.
@Solomon , I believe the use of the lognormal prior in the text is for population-level effects. I’m not sure how to do the same thing for group-level effects.
I tried the lognormal prior on the group-level class = sd and I still got negative values for the individual estimates.

#Using simulated data from first post:
m2 <- brm(response ~ 0 + (1|condition),
          data = ABC,
          iter = 5000, 
          refresh = 0,
          prior = c(prior(lognormal(0,1), class = sd),
                    prior(student_t(3, 0, 1), class = sigma)
                    )

I believe that prior just affects the sd value, not the actual group-level estimates, and the sd value is restricted to be positive anyways (so using a normal prior on the group-level sd is really a half-normal).

I did try the lognormal() distribution family, and this actually gave me negative values:

m3 = brm(
  response ~ 0 + (1|condition),
  family = lognormal(),
  data = ABC,
  refresh = 0,
  prior = c(
    prior(student_t(3, 0, 1), class = sd),
    prior(student_t(3, 0, 1), class = sigma))
  )

I “think” the correct thing to do in this case is take the exponent of the estimate? If so, this is the plot…

m3 %>%
  spread_draws(r_condition[condition,]) %>%
  mutate(e_r_condition = exp(r_condition)) %>% 
  median_qi(e_r_condition, .width = c(.95, .66)) %>% 
  ggplot(aes(y = condition, x = e_r_condition, xmin = .lower, xmax = .upper)) +
  geom_pointinterval()

Like using the exgaussian() family (in post 1) using the lognormal() family worked with this simulated data (and my real data), however the error of the estimates seems to be on the log scale.

Does this seem like a good way to proceed, or are there other options for restricting the group-level estimates to positive values? Thanks!

[Solomon]

Ah. Yes, I was answering your question as if you were referring to a population parameter, prior of class = b. Putting a lognormal on a prior of class = sd won’t work. As you said, there’s already a lower bound of zero on those parameters. If you want to constrain the group-level effects to be non-negative, I think you’ll want to use a link function such as the log link, kinda like if you were doing a Poisson regression. So, yes, I believe your family = lognormal() approach would be one reasonable approach.

[jonnations]

Great, thanks! Your mention of the log link led me to what I believe is the best solution: using the log link in the gaussian family, then taking the exponent of the results. Much obliged.

m4 = brm(
  response ~ 0 + (1|condition),
  family = gaussian(link = "log"),
  data = ABC,
  refresh = 0,
  prior = c(
    prior(student_t(3, 0, 1), class = sd),
    prior(student_t(3, 0, 1), class = sigma))

These posterior distributions mirror the overly-large distributions from the m3 figure above (where I took the exponent of the posteriors from the lognormal family).
Next:

m4 %>%
  spread_draws(r_condition[condition,]) %>%
  mutate(e_r_condition = exp(r_condition)) %>% 
  median_qi(e_r_condition, .width = c(.95, .66)) %>% 
  ggplot(aes(y = condition, x = e_r_condition, xmin = .lower, xmax = .upper)) +
  geom_pointinterval()

And BINGO! Posteriors that represent the actual means from the simulated data, with values > 0

Here’s a figure using this method from the empirical example shown in the first post

[Solomon]

Nice! I’m just barely getting fluent with methods like these. It’s really cool to seem them work out in the wild.