Calculating model averaged estimates using loo

Hi,

I am trying to use pseudo-BMA as described in “Using Stacking to Average Bayesian Predictive Distributions” by Yao et al 2018.

Following R code is directly from “loo_model_weights” function help page

library(rstan)

# generate fake data from N(0,1).
N <- 100
y <- rnorm(N, 0, 1)

# Suppose we have three models: N(-1, sigma), N(0.5, sigma) and N(0.6,sigma).
stan_code <- "
  data {
    int N;
    vector[N] y;
    real mu_fixed;
  }
  parameters {
    real<lower=0> sigma;
  }
  model {
    sigma ~ exponential(1);
    y ~ normal(mu_fixed, sigma);
  }
  generated quantities {
    vector[N] log_lik;
    for (n in 1:N) log_lik[n] = normal_lpdf(y[n]| mu_fixed, sigma);
  }"

mod <- stan_model(model_code = stan_code)
fit1 <- sampling(mod, data=list(N=N, y=y, mu_fixed=-1))
fit2 <- sampling(mod, data=list(N=N, y=y, mu_fixed=0.5))
fit3 <- sampling(mod, data=list(N=N, y=y, mu_fixed=0.6))
model_list <- list(fit1, fit2, fit3)
log_lik_list <- lapply(model_list, extract_log_lik)

# optional but recommended
r_eff_list <- lapply(model_list, function(x) {
  ll_array <- extract_log_lik(x, merge_chains = FALSE)
  relative_eff(exp(ll_array))
})

# stacking method:
wts1 <- loo_model_weights(
  log_lik_list,
  method = "stacking",
  r_eff_list = r_eff_list,
  optim_control = list(reltol=1e-10)
)
print(wts1)

From this, I obtain the weights for each model.

My question is: How can I calculate the model averaged posterior estimates, in this example “model averaged” estimate of the sigma parameter?

(Although this entry (Implementing Model Averaging) looks similar, unfortunately still I did not understand it.)

Thank you very much for your help,

Best regards,

Burak Kürsad Günhan

Please use stacking instead of Pseudo-BMA, as it’s much better.

# stacking method:
wts1 <- loo_model_weights(
  log_lik_list,
  method = "stacking",
  r_eff_list = r_eff_list,
  optim_control = list(reltol=1e-10)
)

The code example you listed is using stacking instead of Pseudo-BMA.

Stacking (and pseudo-BMA) can be used to combine predictive distributions, but it’s sensible for averaging parameters only in some special cases. It’s likely that averaging sigmas doesn’t make sense. Instead you should sample from the combined predictive distribution and focus on the predictions.

Sorry that I confused with “stacking” and “pseudo BMA”.

I see, it is about combining predictive distributions.

That means, I need to have “y_rep” in generated quantities of my Stan model, and I need to fit the models again, but this time using “weighted” number of iterations, and then combine all “y_rep” at the end. And this is my combined posterior distribution for “y_rep”, is that right? Is this what you explained here (Implementing Model Averaging), (but for rstanarm mdoels) right?

Thank you again.

1 Like

No problem, there are many things in the paper.

Yes! That is a correct description and correct link.

Aki

Great, thank you!

1 Like