How to estimate sigma across different treatment groups within this model? (brms)

Modeling
1

Hi there. I am trying to model sigma ~ treatment group (wherein there are 2 groups: supplemented and control) within a model estimating the effect of treatment group on count of surviving offspring to fledging. I read online that I should include the sigma model portion outside of the ‘bf’ but the model I’m running doesnt currently give me output for sigma estimates at all.

I need help to modify my model, or additional code I can run, so that I can have output for sigma ~ treatment group.

My current model is as follows:

total_num_fledge ← brm(
bf(total_fledged ~ 1 + ftreatment + brood_size_sc + hatch_date_sc + # brood size is initial brood size, correct one to use in this model
(1|site) + (1|yearsite) + (1|fyear)),
sigma ~ (1|ftreatment), # need to specify outside bf function to include in model with poisson distribution
data = assumed_fledge_data,
prior = c(set_prior(“normal(0,1)”)),
family = poisson(),
warmup = 1000,
iter = 5000,
chains = 2,
cores = 2,
control = list(adapt_delta = 0.99, max_treedepth = 15),
backend = “cmdstanr”)

summary(total_num_fledge)

2

Unless I’m missing something here, there’s no sigma argument for the Poisson distribution because the Poisson only has one parameter (its expected value).

In general the form for these is something like

brm(bf(response ~ predictor, sigma ~ predictor), family = gaussian(), ...)

3

Since your response variable is a count, you might also consider using the negative binomial likelihood (family = negbinomial). Though it doesn’t contain a \sigma parameter, the negative binomial does have an overdispersion parameter called shape within brm(). You could fit a distributional negative binomial model where shape varies by your treatment group.

2 Likes
4

Oh interesting okay. I haven’t used that before. I would just swap out the family and then sigma for shape? I.e., shape ~ 1+ ftreatment? Or something else?

1 Like
5

From my understanding, Poisson would limit sigma to between 0 and 1 but there can still be variation within that. However, its not unlikely that im incorrect as im quite new to mixed models and brms in general

6

Take a look at the parameterization of the response distributions in brms (https://cran.r-project.org/web/packages/brms/vignettes/brms_families.html#binary-and-count-data-models). As @Solomon suggested, you’ll need a more flexible distribution then the Poisson to model varying dispersion.

2 Likes
7

Yes, you got it.

1 Like
8

The Poisson distribution (at least the Poisson implemented in brms) is defined so that the mean is equivalent to the variance. It cannot vary within some range.

Fortunately, you have a few off-the-shelf options in brms to accommodate discrete data like number of offspring. As noted by @Solomon , the negative binomial is a natural choice if your data are over-dispersed relative to the Poisson. A good way to evaluate under/over-dispersion is with ppc_rootogram() from the tidybayes package.

In my experience working with number of offspring, I often find that the number of young fledged is actually under-dispersed relative to the Poisson. If that’s the case with your data, I can offer a couple of alternatives. First, the Conway-Maxwell Poisson distribution is available in brms as an “experimental option” using family = brmsfamily("com_poisson"). This distribution can accommodate both under- and over-dispersion. Please be aware that the current implementation uses a parameterization of the COM Poisson that offers much slower sampling, although that may not be an issue if your dataset is relatively small (hundreds of rows).

Another option would be to use a cumulative logistic regression, which tends work pretty well for under-dispersed data with a limited number of potential outcomes (for example, say between 1 and 4 fledglings). An example from the literature is available here.

One final thought would be a word of caution about the kind of model you propose to fit (or at least my understanding of it), using a grouping variable (1|ftreatment)in the sigma/shape part of the model. How many treatment groups do you have? If it’s just a few, then the posterior distribution on the standard deviation/variance of the distribution of treatments that impact the variance of the response is likely to be poorly, and most probably very poorly, defined. You can mitigate this issue to some extent by using an informative prior, but my suspicion would be that the prior will be doing some very heavy lifting, if the model converges at all.

I have some other suggestions about your call to brm in the original post. First, I would recommend running fewer post-warmup iterations, spread across more chains. The nature of the HMC sampling in Stan is that you need far fewer posterior samples than in alternative MCMC samplers (JAGS, NIMBLE, etc.). The defaults are generally pretty useful (warmup = 1000, iter = 2000, chains = 4), at least as a starting point. Some of the diagnostics will improve as you add more chains.

I’d also encourage you to think more carefully about and define your priors explicitly, particularly for the group-level variables site, yearsite and fyear. To do this, you might want to run some prior predictive checks to make sure that your prior push-forward distributions match your expectations about how many fledglings are possible. You can do this by setting sample_prior = "only".

I know there’s a lot to absorb above. I’m happy to clarify any points of confusion or elaborate if necessary.

2 Likes