Quadratic effects in brms

Interfaces brms
,
1

Hi all,

I am a PhD student using brms to run mixed models to understand the effect of parental factors on the hatching success of eggs. On of my variables is age, which has been shown to have a non-linear relationship with hatchabiity in species other than the one I am studying. As a result, I want to include the age effects as quadratic effects in my model. My model currently looks like this:

model_hatch <- brm(Hatched ~ 
	#quadratic fixed effects
                     Sire.age.years + Dam.age.years + 
                     I(Dam.age.years^2) + I(Sire.age.years^2) +
	#linear fixed effects
                     Dam.Inbreeding + Sire.Inbreeding + Pair.Inbreeding + 
                     Sire.Inbreeding*Sire.age.years + Dam.Inbreeding*Dam.age.years +
	#random effects
                     (1| Sire.Studbook.ID) + (1| Dam.Studbook.ID)+
                     (1|Lay.Year) + (1| Location),
                           data = data_dev_hatch_sc,
                           family = "bernoulli",
                           prior = priors_hatch,
                           iter = 50000,
                           control=list(adapt_delta=0.99),
                           cores = 8,
                           sample_prior = TRUE)

Based on other research my priors are as follows:

priors_hatch <- c(set_prior("normal(-5.8,5)", class = "b", coef = "Dam.age.years"),
                  set_prior("normal(-5.8,5)", class = "b", coef = "Dam.age.years"),
                  set_prior("normal(-5.8,5)", class = "b", coef = "Sire.age.years"),
                  set_prior("normal(0,5)", class = "b", coef = "Sire.Inbreeding"),
                  set_prior("normal(-5,5)", class = "b", coef = "Dam.Inbreeding"),
                  set_prior("normal(-0.5,5)", class = "b", coef = “Pair.Inbreeding"))

I would like to set priors for these quadratic effects that are different from their linear counterparts but so far am struggling to do so. Is it possible to do this within brms?

Thanks

  • Operating System: macOS Monterey
  • brms Version: 2.18.0
2

What exactly is the issue? - that you can’t work out how to refer to the name of the quadratic part in the priors? If that’s a case a quick fix might be to explicitly code a new variable in your dataframe that is the age squared and give it it’s own name. That was you can just explicitly refer to it in the refression formula and the priors

3

Say you’ve saved your model formula as my_formula. You can learn how brms refers to the various model coefficients for setting priors by executing this code:

get_prior(
  data = data_dev_hatch_sc,
  family = "bernoulli",
  formula = my_formula
)
1 Like
4

Thanks for this - I have tried using a column with squared values but I get very different results when using that new column to the original model (with same seed and no priors).

5

Cheers - when I run this on the original model, the coefficients listed for the quadratic effects are

IDam.age.yearsE2 
ISire.age.yearsE2

It’s when I try and put these into set_prior that I get the problem.

6

Do either of the following code/prior formats work for you? Both these models run for me with a specific prior on the polynomial, maybe you see if there are some differences with how you’ve done the code that might cause an issue:

test_tibble <- tibble(var_a = rnorm(100),
                      var_b = rnorm(100))


get_prior(var_a ~ var_b + I(var_b^2),
          test_tibble,
          gaussian())

priors <- 
  c(set_prior("normal(0 , 1)", class = "Intercept"),
    set_prior("normal(0 , 1)", class = "b", coef = "Ivar_bE2"),
    set_prior("normal(0 , 1)", class = "b", coef = "var_b"))

brm(formula = var_a ~ var_b + I(var_b^2),
     data = test_tibble,
     prior = priors,
     family = "gaussian")

get_prior(var_a ~ poly(var_b, 2),
          test_tibble,
          gaussian())

priors <- 
  c(set_prior("normal(0 , 1)", class = "Intercept"),
    set_prior("normal(0 , 1)", class = "b", coef = "polyvar_b21"),
    set_prior("normal(0 , 1)", class = "b", coef = "polyvar_b22"))

brm(formula = var_a ~ poly(var_b, 2),
    data = test_tibble,
    prior = priors,
    family = "gaussian")
2 Likes
7

Using this fixed it. Thanks for your help!