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