# Model averaging in brms

**URL:** <https://discourse.mc-stan.org/t/model-averaging-in-brms/22755>\
**Category:** General\
**Tags:** loo\
**Created:** [May 28, 2021, 8:27am UTC](https://discourse.mc-stan.org/t/model-averaging-in-brms/22755 "2021-05-28T08:27:39Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![msfarhadinia](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/msfarhadinia/32/11806_2.png) [@msfarhadinia](https://discourse.mc-stan.org/u/msfarhadinia)\
**Post date:** [May 28, 2021, 8:27am UTC](https://discourse.mc-stan.org/t/model-averaging-in-brms/22755/1 "2021-05-28T08:27:39Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [June 2, 2021, 7:39am UTC](https://discourse.mc-stan.org/t/model-averaging-in-brms/22755/2 "2021-06-02T07:39:53Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![msfarhadinia](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/msfarhadinia/32/11806_2.png) [@msfarhadinia](https://discourse.mc-stan.org/u/msfarhadinia)\
**Post date:** [June 8, 2021, 5:10pm UTC](https://discourse.mc-stan.org/t/model-averaging-in-brms/22755/3 "2021-06-08T17:10:42Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [June 15, 2021, 1:47pm UTC](https://discourse.mc-stan.org/t/model-averaging-in-brms/22755/4 "2021-06-15T13:47:30Z")

</div>

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.