Comparing variance or sd of parameters in brms with negative binomial family

I’m relatively new to brms. I’m dealing with one of the studies by my student.

Question: The experiment has different conditions (treatment vs control), type (four levels of exposure in sequence; the four levels are counter-balanced across subjects), and phase (baseline, experiment, post experiment). I have count data (number of vocalizations) as response variable. I’m using negative binomial distribution as family as the response data has 0s and also very high values (like 250). Here’s the prior and model

priors_fit <- c(
  prior(normal(1, 1), class = Intercept),      
  prior(normal(0, 1), class = b),            
  prior(exponential(2), class = sd),           
  prior(gamma(4, 0.1), class = shape)
)



model_fit <- brm(
  formula = vocal ~ condition * type * phase + sex + (1 + phase | ID),
  data = w_data,
  family = negbinomial(),
  prior = priors_fit,
  chains = 4,
  iter = 4000,
  cores = 4,
  file = "model_fit" 
)
  1. Is this model structure correct?
  2. Though I could extracted how population and group level parameters vary in different conditions, types or phase, I’m more interested in comparing the variance or standard deviation across types. I have looked at this Estimating Distributional Models with brms but I’m not sure how to do it for negative binomial distribution. I tried doing similar but the models didn’t converge. Here’s one of the model

variance_model <- brm(
  bf(vocal ~ condition * type * phase + sex + (1 + phase | ID),
     shape ~ type), 
  family = negbinomial(),
  data = w_data,
  prior = c(
    prior(normal(1, 1), class = Intercept),      
    prior(normal(0, 1), class = b),            
    prior(exponential(2), class = sd),           
    prior(gamma(4, 0.1), class = "b", dpar = "shape",lb = 0)
  ),
  control = list(adapt_delta = 0.95),
  cores = 4,
  chains = 4,
  iter = 4000
)

I think, I’m doing something wrong. Maybe the priors are not correct for negative binomial or the fit. Ideally what I want is, comment something like: the variance of response across four levels of types are different (or similar) and which two levels are similar or different.

The idea is that each level of exposure in type can produce different variance of response because of uncertainty perceived by the subjects.

Can someone also guide me to some resources that deal with similar kind of problems? And the interpretation of the output for the model?

Thank you!

  1. I don’t know that I can answer your question in a direct way, but your model formula seems reasonable given your description.
  2. I agree a distributional NB model is an excellent choice for your data. When you fit a distributional NB model, keep in mind brm() will automatically switch to the log link for the shape model. The gamma() prior would be reasonable for shape on the identity scale, but it’s not a good idea for the log scale. I recommend you use a normal or student-t prior instead. There are ways to force brm() to use the identity scale instead, but I think the default is wise in this case.

It looks like this is your first post. Welcome to the Stan forums, @woodpecker!

Thank you! It works; though the model has few (5-6) divergent transitions but good R hat and ESS.

On interpretation: I believe higher the shape parameter, more concentrated the values. Would it make sense to take the inverse of shape parameter to show dispersion? Like obtain the posterior draws and get 1/exp(b_shape_Intercept+b_shape_factor1), 1/exp(b_shape_Intercept+b_shape_factor2) etc. and plot the density to show the differences

For NB as parameterized in brms, greater overdispersion means shape goes to zero, and shape goes to infinity for equidispersion. That is, the higher the shape value, the closer the data are to the Poisson ideal.

Kind of makes sense. Do you think taking an inverse will make it more intuitive to visualize?

Seems like that’s more of an issue of what you and your colleagues would find the most informative. Personally, I wouldn’t go through the effort.

1 Like

In the PC prior framework a prior is placed on the reciprocal or the square root of it. I think the brms default prior for the shape is a reasonable gamma prior that puts more mass on higher values. I’d think a normal or student-t centered at 0 is going to prefer much higher overdispersion.

1 Like

Yeah. I tried various priors, student_t(3,0,3) worked for the data I have.

brms default approximates a PC prior, see Default prior for negative-binomial shape parameter

But as @Solomon mentioned that default is on the identity scale.

1 Like

Sorry about that. I recall it being something different! Great notebook.

Ah, it was changed in 2024 based on my suggestion A better default prior for shape parameter of negative binomial · Issue #1614 · paul-buerkner/brms · GitHub

1 Like