Search code examples
rmixed-modelsrandom-effectsrstan

Defining prior for both random effects and random effects variance in mixed effect model (R brms)

I want to fit GLMM Poisson model counts. I have 121 subjects (subject), I observe 8 Poisson counts (count) per subject: they correspond to 2 types of events (event) x 4 periods (period). 

'data.frame':   968 obs. of  4 variables:
 $ count  : num  4 0 2 3 3 0 1 0 8 14 ...
 $ subject: num  1 1 1 1 1 1 1 1 2 2 ...
 $ event  : Factor w/ 2 levels "call","visit": 2 2 2 2 2 2 2 2 1 1 ...
 $ period : Factor w/ 4 levels "period_1","period_2",..: 1 2 3 4 1 2 3 4 1 2 ...

What I want:

  • I want to fit GLMM with Bayesian approach, assuming random effects for (1) subject, (2) subject:visit, (3) subject:event.
  • I want to set N(0, 10^6) prior for all fixed effects,
  • I want to set N(0, sigma2_a), N(0, sigma2_b), N(0, sigma2_c) priors for random effects for (1) subject, (2) subject:visit, (3) subject:event, respectively,
  • I want to set uniform priors for sigma2_a, sigma2_b, sigma2_c respectively.

What I managed to get:

  • I believe I am setting the model formula properly, and I am setting the desired priors for fixed effect parameters: 

    ## define some priors          
prior <- c(prior_string("normal(0,10^3)", class = "b"),
       prior_string("normal(0,10^3)", class = "b", coef = "eventvisit"),
       prior_string("normal(0,10^3)", class = "b", coef = "eventvisit:periodperiod_2"),
       prior_string("normal(0,10^3)", class = "b", coef = "eventvisit:periodperiod_3"),
       prior_string("normal(0,10^3)", class = "b", coef = "eventvisit:periodperiod_4"),
       prior_string("normal(0,10^3)", class = "b", coef = "periodperiod_2"),
       prior_string("normal(0,10^3)", class = "b", coef = "periodperiod_3"),
       prior_string("normal(0,10^3)", class = "b", coef = "periodperiod_4"))

## fit model
fit1 <- brm(count ~ event + period + event:period + (1|subject) + (0 + event|subject) + (0 + period|subject),
    data = data.long, family = poisson(),
    prior = prior,
    warmup = 1000, iter = 4000, chains = 4,
    cores = 7)

What I struggle with:

  • How to set N(0, sigma2_a), N(0, sigma2_b), N(0, sigma2_c) priors for random effects for (1) subject, (2) subject:visit, (3) subject:event, respectively,
  • How to set uniform priors for sigma2_a, sigma2_b, sigma2_c respectively.
Solution

  • With regard to what you want: Could you please explain in more detail what you want with subject:visit and subject:event? Currently, I cannot tell whether your model is correctly specified for your purpose.

    With regard to the priors:

    1. If you set prior_string("normal(0,10^3)", class = "b"), there is no need to specifiy that prior again for each coefficient. Setting it globally for the whole prior class is sufficient. Note that, given the scale of your response and predictor variables, normal(0, 1000) is so wide that it has practically no influence. You migh as well just leave it out.
    2. Those normal(0, sigma_x) priors are set automatically for all random effects terms. No need to worry about them.
    3. You can set priors on random effects SDs using class = "sd". For instance prior_string("cauchy(0, 5)", class = "sd"). I explicitely encourage you not to use uniform priors in particular because they have a hard upper bound for a parameter (standard deviation), which has not theoretical upper bound.