Prior upper bound with an interaction?

Interfaces brms
1

I’m fitting models where it’s important to constrain some parameters to have an upper bound of zero. I’ve been able to do that after some help from this post from a while ago. However, if I have an interaction with a categorical variable, this can sometimes cause one or more levels of the interaction to end up with a positive coefficient. This ultimately seems to be based on what the reference level is and whether the interaction brings a coefficient up or down from the reference level.

One solution is to change the order of the factor levels and manually select the most appropriate reference level. However, there could end up being problems if the best reference level isn’t the same for all effects that need to be constrained below zero. Is there any way to set the interaction prior such that its upper bound does not exceed the (negative) reference level coefficient value? I guess more broadly, the question is whether you can set a prior’s upper bound based on the value of another coefficient?

The test code below shows an example of this issue and how it can be fixed by reordering the factor levels. Any help would be greatly appreciated!

library(brms)
library(tidyverse)

#Create a dataset that includes positive and negative slopes based on group
n <- 100
int <- c(-3, -4)
X1 <- (c(-0.5, 0.1))
sig <- 0.5

Y <- numeric()
Yreal <- numeric()
X <- numeric()
Group <- character()

for (i in 1:length(X1)){
  
  x <- rnorm(n, 10, 2)
  X <- c(X,x)  
  Group <- c(Group, rep(LETTERS[i], n))
  
  y <- int[i] + X1[i]*x
  
  Yreal <- c(Yreal, y)
  Y <- c(Y, y + rnorm(n, 0, sig))
  
}

df <- tibble(Group = Group,
             X = X,
             Y = Y,
             Yreal = Yreal)

ggplot(df, aes(x = X, y = Y, color = Group)) +
  geom_point() +
  # stat_smooth(method = lm) +
  geom_line(aes(x = X, y = Yreal, group = Group), inherit.aes = FALSE, linetype = "dashed") +
  labs(subtitle = "Dashed lines are actual relationships")


#Fit model with X prior upper bound of zero
test.mod1 <- brm(
  formula = bf(Y ~ intercept +
                 x * X +
                 group * Group +
                 xgroup * X * Group,
               intercept ~ 1,
               x ~ 1,
               group ~ 1,
               xgroup ~ 1,
               nl = TRUE),
  data = df,
  family = gaussian(),
  prior = c(prior(normal(0,100), class = "b", nlpar = "intercept", ub = NA),
            prior(normal(0,100), class = "b", nlpar = "x", ub = 0),
            prior(normal(0,100), class = "b", nlpar = "group", ub = NA),
            prior(normal(0,100), class = "b", nlpar = "xgroup", ub = NA)),
  chains = 4,
  iter = 5000,
  cores = 4,
  save_pars = save_pars(all = TRUE)
)

#Group A is the reference level and has a negative coefficient (satisfies prior ub = 0 constraint)
fixef(test.mod1)

#The slope for Group B ends up positive after adding the xgroup_Intercept value
#The xgroup prior would need to be specified as having an upper bound of -x_Intercept
sum(fixef(test.mod1)[c(2,4),1])

plot(test.mod1)

#Plot results
df$pred1 <- apply(posterior_epred(test.mod1), 2, mean)

ggplot(df, aes(x = X, y = Y, color = Group)) +
  geom_point() +
  geom_line(aes(x = X, y = pred1), linewidth = 1.5) +
  geom_line(aes(x = X, y = Yreal, group = Group), inherit.aes = FALSE, linetype = "dashed") +
  labs(subtitle = "Model slopes should never be positive")


#Reverse the order: change the reference level
df$Group <- factor(df$Group, levels = c("B", "A"))

test.mod2 <- brm(
  formula = bf(Y ~ intercept +
                 x * X +
                 group * Group +
                 xgroup * X * Group,
               intercept ~ 1,
               x ~ 1,
               group ~ 1,
               xgroup ~ 1,
               nl = TRUE),
  data = df,
  family = gaussian(),
  prior = c(prior(normal(0,100), class = "b", nlpar = "intercept", ub = NA),
            prior(normal(0,100), class = "b", nlpar = "x", ub = 0),
            prior(normal(0,100), class = "b", nlpar = "group", ub = NA),
            prior(normal(0,100), class = "b", nlpar = "xgroup", ub = NA)),
  chains = 4,
  iter = 5000,
  cores = 4,
  save_pars = save_pars(all = TRUE)
)

#Group B is the reference level and its coefficient is held at ~zero (satisfies prior ub = 0 constraint)
fixef(test.mod2)

#No issue for Group A as xgroup_Intercept makes it more negative
sum(fixef(test.mod2)[c(2,4),1])

plot(test.mod2)

#Plot results
df$pred2 <- apply(posterior_epred(test.mod2), 2, mean)

ggplot(df, aes(x = X, y = Y, color = Group)) +
  geom_point() +
  geom_line(aes(x = X, y = pred2), linewidth = 1.5) +
  geom_line(aes(x = X, y = Yreal, group = Group), inherit.aes = FALSE, linetype = "dashed") +
  labs(subtitle = "Reversing the group order (reference level) keeps slopes <= 0")
  • Operating System: CentOS 7
  • brms Version: 2.19.0