R/loo_model_weights.R
loo_model_weights.Rd
Model averaging via stacking of predictive distributions, pseudoBMA weighting or pseudoBMA+ 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, ...) # S3 method for default 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 loglikelihood 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 pseudoBMA+ 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 
A numeric vector containing one weight for each model.
loo_model_weights()
is a wrapper around the stacking_weights()
and
pseudobma_weights()
functions that implements stacking, pseudoBMA, and
pseudoBMA+ weighting for combining multiple predictive distributions. We can
use approximate or exact leaveoneout crossvalidation (LOOCV) or Kfold 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 leaveoneout
predictive density of the combination distribution. That is, it finds the
optimal linear combining weights for maximizing the leaveoneout log score.
The pseudoBMA 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 pseudoBMA+ 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 pseudoBMA+ can serve as a computationally easier alternative.
Vehtari, A., Gelman, A., and Gabry, J. (2017a). Practical Bayesian model evaluation using leaveoneout crossvalidation and WAIC. Statistics and Computing. 27(5), 14131432. doi:10.1007/s1122201696964 (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/17BA1091. (online).
The loo package vignettes, particularly Bayesian Stacking and PseudoBMA weights using the loo package.
loo()
for details on leaveoneout ELPD estimation.
constrOptim()
for the choice of optimization methods and controlparameters.
relative_eff()
for computing r_eff
.
if (FALSE) { ### 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=1e10) ) 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=1e10) ) all.equal(wts1, wts2) # pseudoBMA+ method: set.seed(1414) loo_model_weights(loo_list, method = "pseudobma") # pseudoBMA 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) }