Model averaging in brms

General
1

I am running 9 GAM models with brms, two models got all the model weights, so I am trying to do model averaging between the the top two models.

model1 ← brm(bf(y ~ s(x1)+ (Grouping variable)),
data = data, prior = prior, family = gaussian(),
cores = 4, iter = 4000, warmup = 2000,
control = list(adapt_delta = 0.99))

model2 ← brm(bf(y ~ s(x2)+ (Grouping variable)),
data = data, prior = prior, family = gaussian(),
cores = 4, iter = 4000, warmup = 2000,
control = list(adapt_delta = 0.99))

pred_avg ← pp_average(model1, model2, weights = “loo”, method = “predict”)

Then, to get the posteriors:

posterior_summary(pred_avg, pars = c("^b_", “^sd_”, “sigma”, “^bs_”), probs = c(0.025, 0.975) )

                 Estimate  Est.Error       Q2.5       Q97.5

Estimate -0.4209933 0.24448195 -0.8014096 0.01751906

Est.Error 0.6804755 0.03969645 0.6372006 0.77561411

Q2.5 -1.7658316 0.23672804 -2.1557231 -1.34990816

Q97.5 0.9216051 0.28203470 0.5136583 1.44436048

First question: I need the coef for each covariate to report, which seems missing here. Where am I wrong here?

Second question: How can I plot the results of the final model averaged object?

I used stanplot (pred_avg, type = “hist”) and mcmc_plot(pred_avg), none working with the model averaged object though. Any suggestion?

Third question: I need to plot the predictions, but seems that marginal_effect function again does not work with the model averaged object. What is the solution?

Many thanks for your advice

2

Welcome to the Stan forums! I am having a hard time reading your output of posterior summary. Can you reformat that a bit. You can wrap chunks of code in ``` to make it more readable. pp_average averages samples from the posterior predictive distribution not from the posterior distribution itself. You can try to get the latter via posterior_average but note that you need to make sure that the parameters with the same names are actually comparable across models.

With regard to plotting: You need to use methods of the bayesplot package directly as the averaged posterior is simple a data.frame not a brmsfit object anymore. With regard to conditional_effects I am not sure there is a nice alternative for the models averaged that way.

Next time you ask a brms related question, I recommend adding the “interfaces - brms” tag

1 Like
3

Thanks Paul, very helpful.

I am running 9 GAM models, two of them possess all the loo weight (m1 and m3).

m1 <- brm(bf(Habitat.change ~ s(MCI)+ Region),
            data = data, 
          prior = prior, 
          family = gaussian(), 
          cores = 4, iter = 4000, warmup = 2000,
            control = list(adapt_delta = 0.99))

m3 <- brm(bf(Habitat.change ~ s(Ranger.1000km2)+ Region),
          data = data, 
          prior = prior, 
          family = gaussian(), 
          cores = 4, iter = 4000, warmup = 2000,
          control = list(adapt_delta = 0.99))

Each of them has a single predictor + one grouping variable.
So, I ran the below as you suggested:

posterior_average(m3, m1)

But, as the grouping variable (Region) is common between two models, the model averaging provides only predicted posteriors for Region, not the other two variables.

b_Intercept b_Region1 b_Region2 b_Region3
-1.054897389 0.2553771967 0.212400930 0.6342051572
-0.459273627 0.2133165428 0.217377924 0.5728726600
-0.682074405 0.0861842014 0.580051210 -0.1638639849
-0.522941073 0.0177629059 0.141624344 -0.1380124456
-1.078018820 0.2185330513 0.862905920 0.6814814559
-0.679868545 0.2450694988 0.729577995 0.1356133630
-0.737136057 0.0190874273 -0.177539907 0.2043542725
-0.944302943 0.0409637633 0.732163755 0.5875435218
-0.982298403 0.5107304236 0.480600404 0.4748065664
-0.871294004 0.3490444618 0.782958540 0.3190764625
…
1-10 of 8,000 rows | 1-4 of 7 columns
Show in New Window

Now, I need to report B coef, but seems that I need to go two supplementary ways:

  1. For variables which are not shared between models, do I need to report the outcomes from their respective models, instead of the averaged model?
    For example, to get the estimates for MCI, I have to run summary(m1). Am I right?

  2. Then, for variables which are shared between models entering the model averaging, as Region (grouping variable) here, we can see the predictions for y, rather than the B estimates for parameters when running posterior_average(m3, m1). So, how to get B coef for the parameters which are shared between models (Region in this case)?

  3. How to create a plot for all these three variable B coefs in single plot?

Many thanks
Mohammad

4

  1. I don’t have a good answer for this to be honest. It depends on what the model not estimating that parameter implicitely assumes as fixed value, for example, 0 for a regression parameter of a variable excluded from the model. That way, we would still have to average although one model would have a fixed parameter value. That is nothing brms can do for you and you have to do it manually.

  2. I am sorry, I don’t understand this question.

  3. I don’t have a good recommendation for this, but I would assume many ways to plot those variables do exist since you have the full set of posterior samples for these model averaged parameters.