Multinomial model in brms

Hello Paul,
We were wondering if it was possible to estimate parameters using brms
for the following problem and data setup:
We have the following data

For each “Ethnicity” I have a count associated with the number of
copies of an Allele per Locus, call that y_{i,j} where i corresponds
to the Ethnicity and j corresponds to the Allele. We also have n_{i}
where it is simply the sum of the Allele per locus.

We are interested in estimating the following parameters:

  1. correlation of frequencies across Ethnicities,
  2. posterior estimates of Allele frequencies per Ethnicity.

According to your e-mail response:
you using a multinomal model.

So if A is the matrix containing the ethnicities in the rows (rows of the data) and the j allele in the columns, n is the total number of obs per row, and “obs” is a simple row index, you can go for something like

brm(A | trials(n) ~ 1 + <predictors> + (1 | ID | obs), ...)

The (1 | ID | obs) models a latent normal overdisperson (instead of a beta one as in the beta-binomial model)

and further the | ID | ensures that those are correlated across alleles.

My remaining question is what should ... be?

Best regards,
Manuel

1 Like

Hi Paul,

This is Yosuke, a trainee working with Maunel. For the model you’ve suggested, we ended up using this one:

fit <- brm(
    response | trials(n) ~ 1 + (1 | ID | obs),
    family=multinomial(),
    data=A,
    control=list(adapt_delta = 0.95)
)

We are trying to look at correlation across K alleles. It seems like the multinomial model estimates the correlation among K-1 parameters (which is natural in the multinomial setting), which we can extract with vcov.brmsfit(fit, correlation = TRUE).

For the visualization purpose, we are interested in showing correlation across K x K classes.
What would be the recommendation to recover the full correlation matrix?

Thank you very much for your help.

Sincerely yours,
Yosuke

Hi Paul,

We are trying to fit multinomial model using brms and it seems like we are not able to get a stable fit. We would love to get some advise.

Please see our Gist here: https://gist.github.com/yk-tanigawa/aeb3275bd7e93af009d39daed6cb8df9

Thank you very much.

Best,
Yosuke

vcov is not the right approach. The correlations you are looking for are shown in the summary output under group-level effects or alternatively extracted via VarCorr.

If you want estimates for the full correlation matrix you need to use family = multinomial(refcat = NA) but be aware that this implies a very hard to fit model which requires informative priors to run well enough.

Thank you, so much!

Thank you very much! We will have a look.