Model averaging/weighting via stacking or pseudo-BMA weighting
Model averaging via stacking of predictive distributions, pseudo-BMA weighting or pseudo-BMA+ weighting with the Bayesian bootstrap. See Yao et al. (2018), Vehtari, Gelman, and Gabry (2017), and Vehtari, Simpson, Gelman, Yao, and Gabry (2019) for background.
loo_model_weights(x, ...)
## Default S3 method:
loo_model_weights(
x,
...,
method = c("stacking", "pseudobma"),
optim_method = "BFGS",
optim_control = list(),
BB = TRUE,
BB_n = 1000,
alpha = 1,
r_eff_list = NULL,
cores = getOption("mc.cores", 1)
)
stacking_weights(lpd_point, optim_method = "BFGS", optim_control = list())
pseudobma_weights(lpd_point, BB = TRUE, BB_n = 1000, alpha = 1)x |
A list of pointwise log-likelihood matrices or |
... |
Unused, except for the generic to pass arguments to individual methods. |
method |
Either |
optim_method |
If |
optim_control |
If |
BB |
Logical used when |
BB_n |
For pseudo-BMA+ weighting only, the number of samples to use for
the Bayesian bootstrap. The default is |
alpha |
Positive scalar shape parameter in the Dirichlet distribution
used for the Bayesian bootstrap. The default is |
r_eff_list |
Optionally, a list of relative effective sample size
estimates for the likelihood |
cores |
The number of cores to use for parallelization. This defaults to
the option
|
lpd_point |
If calling |
loo_model_weights() is a wrapper around the stacking_weights() and
pseudobma_weights() functions that implements stacking, pseudo-BMA, and
pseudo-BMA+ weighting for combining multiple predictive distributions. We can
use approximate or exact leave-one-out cross-validation (LOO-CV) or K-fold CV
to estimate the expected log predictive density (ELPD).
The stacking method (method="stacking"), which is the default for
loo_model_weights(), combines all models by maximizing the leave-one-out
predictive density of the combination distribution. That is, it finds the
optimal linear combining weights for maximizing the leave-one-out log score.
The pseudo-BMA method (method="pseudobma") finds the relative weights
proportional to the ELPD of each model. However, when
method="pseudobma", the default is to also use the Bayesian bootstrap
(BB=TRUE), which corresponds to the pseudo-BMA+ method. The Bayesian
bootstrap takes into account the uncertainty of finite data points and
regularizes the weights away from the extremes of 0 and 1.
In general, we recommend stacking for averaging predictive distributions, while pseudo-BMA+ can serve as a computationally easier alternative.
A numeric vector containing one weight for each model.
Vehtari, A., Gelman, A., and Gabry, J. (2017a). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing. 27(5), 1413–1432. doi:10.1007/s11222-016-9696-4 (journal version, preprint arXiv:1507.04544).
Vehtari, A., Simpson, D., Gelman, A., Yao, Y., and Gabry, J. (2019). Pareto smoothed importance sampling. preprint arXiv:1507.02646
Yao, Y., Vehtari, A., Simpson, D., and Gelman, A. (2018) Using stacking to average Bayesian predictive distributions. Bayesian Analysis, advance publication, doi:10.1214/17-BA1091. (online).
The loo package vignettes, particularly Bayesian Stacking and Pseudo-BMA weights using the loo package.
loo() for details on leave-one-out ELPD estimation.
constrOptim() for the choice of optimization methods and control-parameters.
relative_eff() for computing r_eff.
## Not run:
### Demonstrating usage after fitting models with RStan
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)
# can also pass a list of psis_loo objects to avoid recomputing loo
loo_list <- lapply(1:length(log_lik_list), function(j) {
loo(log_lik_list[[j]], r_eff = r_eff_list[[j]])
})
wts2 <- loo_model_weights(
loo_list,
method = "stacking",
optim_control = list(reltol=1e-10)
)
all.equal(wts1, wts2)
# pseudo-BMA+ method:
set.seed(1414)
loo_model_weights(loo_list, method = "pseudobma")
# pseudo-BMA method (set BB = FALSE):
loo_model_weights(loo_list, method = "pseudobma", BB = FALSE)
# calling stacking_weights or pseudobma_weights directly
lpd1 <- loo(log_lik_list[[1]], r_eff = r_eff_list[[1]])$pointwise[,1]
lpd2 <- loo(log_lik_list[[2]], r_eff = r_eff_list[[2]])$pointwise[,1]
lpd3 <- loo(log_lik_list[[3]], r_eff = r_eff_list[[3]])$pointwise[,1]
stacking_weights(cbind(lpd1, lpd2, lpd3))
pseudobma_weights(cbind(lpd1, lpd2, lpd3))
pseudobma_weights(cbind(lpd1, lpd2, lpd3), BB = FALSE)
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.