Mixed Multinomial Model built with brms::brm: diagnostic and goodness of fit

Modeling
1

Hi, all.
I’m new at Bayesian modelling and I just fitted a Mixed Multinomial model with brms::brm function.
The outputs look like follow:

To Stan Forum
To Stan Forum840×720 150 KB

Please, help me answer the following questions:
1)
As you can see, most of the estimates have very large values: is this a symptom of an issue?
If yes, how can I fix this issue?
In general, what tools can I use to make the diagnostic of such models (built with brm or other bayesian modelling functions)
2)
Could you please suggest some tools to assess the goodness of fit for this model (this type of models)?
3)
It was hard for me to find 02 functions to fit Mixed Multinomial regression: one frequentist (mixcat::npmlt) and one Bayesian (brms::brm)
Any suggestion for other functions and/or documentation to build such models will be much appreciated.

Thanks in advance for your answers.
Best regards.
Romaine

2

Hi Romaine,

I’m super new with Bayesian Stats myself, but from what I’ve gathered:

  1. your Rhat values are way beyond 1.0. Basically, your model hasn’t “converged” in the classical sense, so I would not interpret any of these values at all.

  2. Check out the loo (leave-one-out) package and bayesplot for diagnostic tools. Before doing that though, I’d ensure the Rhat values are at 1.0. I don’t know your data situation, i.e. whether this can be fixed with more iterations / better priors or needs a different formula altogether.

  3. Not sure what to make of the question…

All my best

3

Hi Petersen.
Thank you for your answers that are helpful for me.
Based on your suggestions, I re-ran the model with iter=10000 instead of iter=1000;
but, the convergence issue persists.
About my the data:
their specificty is that, the dependent variable and the Predictor1 are values of the same main variable collected at two different time-periods.
Could this specificity be a possible cause of the convergence issue that is occuring?
Thanks for your answers.
Best regards.

4

I guess your model is overparameterized. Can you try what happens if you remove (1 | grouping_variable)? This would help to see if there are biggers problems with your data. Can you also show how you provide the formula to brms? The output shows the formula partially, but not in the form you have given it to brms.

5

Hello avehtari.
Thank you for your reply.
Please, find below how I provided the formula to brms:

Blockquote
model ← brm(independant.variable~ Predictor1+Predictor2+Predictor3+Predictor4 + (1 | Grouping.variable),
data = data, family = categorical(),
control = list(adapt_delta = 0.9, max_treedepth = 15),iter = 1000, chains = 4,seed=1207,core=4)

I incorporated the random intercept term to account for the dependance of the data; endeed, the dependent variable and the predictor1 had been collected on the same indicviduals at two different time-period.

When I removed the random part (1|grouping_variable), the output looked like follows:

Blockquote

Model output
Model output828×558 125 KB

Thanks in advance for your help.
Best regards.

6

As you still have the problems without the varying intercept term, and looking at the posterior estimates and intervals, it’s likely that you have complete separability in the data, and you need to define proper priors for the coefficients. brms uses flat prior by default, and that doesn’t work if there is complete separability in the data.

7

Thank you so much for your feedback.
1)
Please, which clues in the output make you think of complete separability in the data?
2)As I am new at brms modelling, I don’t know yet how to define priors (especially on the coefficients); so,
could you please help me define priors (especially on the coefficients)?
Any documentation or cases examples could be very helpful.
Thanks in advance for your assistance.
Best regards.

8

To set custom priors, you can add the following inside the brm function when you fit the model:

prior = c(prior(normal(0,10), class="b", dpar=“mu5059Kg”),
          prior(normal(0,10), class="b", dpar="mu6069Kg”),
          prior(normal(0,10), class="b", dpar="mu70Kg”))

This will set the same normal(0,10) prior on all four of the fixed effects (predictors 1 through 4). If you want to set a different prior on one of these, you could do this:

prior = c(prior(normal(0,10), class="b", dpar=“mu5059Kg”),
          prior(normal(0,5), class="b", coef="Predictor1", dpar="mu5059Kg"),
          prior(normal(0,10), class="b", dpar="mu6069Kg”),
          prior(normal(0,10), class="b", dpar="mu70Kg”))

This will set a prior of normal(0,5) on Predictor1 for outcome mu5059Kg and normal(0,10) for all of the other fixed effects.

You can set priors on the other parameters (e.g., the random effects) in a similar fashion. Note that only the fixed effects have flat priors by default. The other parameters have some other default prior (e.g., student t). To see the default priors for each parameter, use the code below. The output will also show you which arguments need to be set in order to change the prior for each of these parameters.

get_prior(independant.variable ~ Predictor1+Predictor2+Predictor3+Predictor4 + (1 | Grouping.variable),
          data = data, 
          family = categorical())

To see the actual priors used in a model you’ve fit:

prior_summary(model)
9

Joels,
Thank you so much for your answers.
From your answers, I have the following question:
In general, on which clues do I based on to choose both distribution and its corresponding parameters to set priors an a model’s parameters?
Thanks in advance for your answer.