How to specify the priors of a random effect ("sds") for a gamm in brms

Interfaces brms
,
1

Hello,

I am trying to fit a GAMM in brms where I specify different priors. The model returns the following error :

Error: The following priors do not correspond to any model parameter:
sds_s(predator_id, bs ="re") ~ normal(0.60, 0.5)
Function 'get_prior' might be helpful to you.

The reason is I cannot find the proper name to write in the “coef=” argument of “set_prior()”. The variable I cannot specify is the random effect called “predator_id”.

Here is the name of the coefficients when I call “variables(fit)”, with "fit being my model :

[1] "b_Intercept"              "b_Zgame_duration"         "bs_sZcumul_xp_1"
[4] "sds_sZcumul_xp_1"         "sds_spredator_id_1"       "phi"
[7] "Intercept"

Here are my priors :

priors <- c(
  # priors on fixed effects
  set_prior("normal(0, 1)",
            class = "b",
            coef = "Zgame_duration"),
  set_prior("normal(0, 2)",
            class = "b",
            coef = "sZcumul_xp_1"),
  # prior on the intercept
  set_prior("normal(0.05, 0.5)",
            class = "Intercept"),
  # priors on smooth terms
  set_prior("cauchy(0, 2)",
            class = "sds",
            coef = "s(Zcumul_xp)"),
  # prior on sds predator_id
  set_prior("normal(0.60, 0.5)",
            class = "sds",
            coef = paste('s(predator_id, bs =', '"re")', sep = "")),         
  # priors on phi
  set_prior("normal(2, 0.5)",
            class = "phi")
            )

When I run “prior_summary(fit)”, I get :

r$> prior_summary(mod1)
                prior     class                      coef group resp dpar nlpar lb ub
               (flat)         b
         normal(0, 2)         b              sZcumul_xp_1
         normal(0, 1)         b            Zgame_duration
         normal(0, 2) Intercept
       normal(2, 0.5)       phi                                                  0   
 student_t(3, 0, 2.5)       sds                                                  0   
 student_t(3, 0, 2.5)       sds s(predator_id, bs = "re")                        0
         normal(0, 1)       sds              s(Zcumul_xp)                        0
       source
      default
         user
         user
         user
         user
      default
 (vectorized)
         user

I thus assumed that the name I should write in the “coef” argument of “set_prior()” would be “s(predator_id, bs = “re”)”, but it does not work. I tried many different combinations with no success. Could someone please help me specify the prior for this random effect?

Here is the model formula :

model_formula <- brmsformula(
  hunting_success | vint(4) ~
      s(Zcumul_xp) +
      s(predator_id, bs = "re") +
      Zgame_duration
)

Thank you very much for your help!

2

Why do you have that in front of the name of the coef in your code to set priors? Could that be causing some problem?
I haven’t had trouble setting priors on the sds in brms before, so this seems strange. Do you have a reproducible example?

3

@Maxime Presumably you would want to include your group-level effect using this:

... + (1|predator_id) ...

sampled from a gaussian distribution (assuming that predator_id is a factor). Other distributional forms for the group-level effects are possible using brms::brm.

… instead of this:

... + s(predator_id, bs = "re") ... , which is very specific to mgcv::gam

… and then setting priors should work.

4

Hi @jd_c ,

I actually tried many solutions (including this one) and it did not work so I don’t think it stems from using paste(). It was to insert a string with double quote marks inside to reproduce the exact name of the coefficient.

I could provide the model fit (.rds file) but cannot provide a reproducible example because the data is private at the moment.

5

Hi @MilaniC,

Thank you for your answer. Using the common (1 |predator_id) would give the exact same results? I thought it was a bit different though I admit I don’t know what would be the difference.

Maxime

6

Good catch by @MilaniC . @Maxime why would you prefer to use GAM implementation instead of simply a varying intercept (1|predator_id)? What would be the advantage of the GAM implementation? They aren’t the same but I think should have similar results. See here Using random effects in GAMs with mgcv

7

@Maxime The link provided by @jd_c provides an informative overview. There are several approaches to addressing random effects or group-level effects. brms::brm uses the lme4 syntax and that is what you need to apply in this case (which is similar to using mgcv::gamm rather than mgcv::gam).

They are two alternative approaches and so will not in most cases give exactly the same estimates for the group-level effects — but they will be similar and not likely to change the overall model inference.

8

Hi @MilaniC and @jd_c,

Maybe I should have added a bit of context. I am actually interested in modelling player expertise in a video game, which is described by a nonlinear relationship between success and experience. Because we hypothesize that the relationship will attain a plateau at some level of experience when the player is advanced, a linear relationship cannot describe what I am aiming to test.

I am aware that brms:brm uses the lme4 syntax and that the non linear syntax of brms:brm employs mgcv::gamm. My question to @MilaniC was more about the technical differences between the syntax for the random effects, but I still don’t see why I should apply the lme4 syntax over mgcv since I am modelling a nonlinear relationship. In my case, I am applying the “GS” model, defined by Perdersen et al. 2019 as “a global smoother plus group-level smoothers that have the same wiggliness” (GS for Global smoother with individual effects that have a Shared penalty) (see Figure 4).

In any case, I found a way to apply my prior by specifying the prior for all the “sds” class before calling the “sds” prior for “s(Zcumul_xp)” (see below).

Thank you very much for the discussion! I am still open for suggestions if you have any.

priors <- c(
  # priors on fixed effects
  set_prior("normal(0, 1)",
            class = "b",
            coef = "Zgame_duration"),
  set_prior("normal(0, 2)",
            class = "b",
            coef = "sZcumul_xp_1"),
  # prior on the intercept
  set_prior("normal(0.05, 0.5)",
            class = "Intercept"),
  # prior on sds predator_id
  set_prior("normal(0.60, 0.5)",
            class = "sds"),
  # priors on smooth terms
  set_prior("cauchy(0, 2)",
            class = "sds",
            coef = "s(Zcumul_xp)"),      
  # priors on phi
  set_prior("normal(2, 0.5)",
            class = "phi")
            )
9

I made a mistake here :

“In my case, I am applying the “GS” model, defined by Perdersen et al. 2019 as “a global smoother plus group-level smoothers that have the same wiggliness” (GS for Global smoother with individual effects that have a Shared penalty) (see Figure 4).”

I am actually modelling model G in this case. I am also doing GS but it is not the model presented in this post

10

Yeah, it looks like model G could simply use the (1|predator_id) varying intercept and get similar results. Model GS is something entirely different and would require the use of bs=fs for a group-level smooth, which gives a different smooth fit for each level of the factor, which as I understand it, is sort of like varying slopes but smooths instead of linear slopes. (depending on the size of the group, this can take a long time to fit!!)