sample.ensemble: Sample from the ensemble posterior distribution

Description

Given model averaging weights (e.g., from Bayesian model averaging (BMA), pseudo-BMA, or stacking) and a matrix of posterior samples from the candidate models, this function draws samples from the model-averaged posterior distribution. Here, each "model" refers to a unique combination of an outcome model and its associated priors. Posterior draws are randomly selected from the candidate models in proportion to their specified weights, producing samples from the ensemble of posterior distributions.

Usage

sample.ensemble(wts, samples.mtx)

Value

The function returns a numeric vector of ensemble posterior draws, sampled proportionally to the provided model weights. The returned vector has the same length as the number of rows in samples.mtx.

Arguments

wts: a numeric vector of normalized model averaging weights (e.g., from compute.ensemble.weights()). The length of wts must match the number of columns in samples.mtx.
samples.mtx: a matrix of posterior samples. Each column corresponds to samples from a different model, and each row is one posterior draw (e.g., from Markov chain Monte Carlo (MCMC) sampling). All columns must have the same number of samples.

Details

This function is typically used in combination with compute.ensemble.weights(), which computes model averaging weights using methods such as Bayesian model averaging (BMA), pseudo-BMA, pseudo-BMA with the Bayesian bootstrap, or stacking). The input matrix of posterior samples should have one column per candidate model, with each column containing posterior draws from that model.

Examples

Run this code

if (instantiate::stan_cmdstan_exists()) {
  if(requireNamespace("survival")){
    library(survival)
    data(E1684)
    data(E1690)
    ## replace 0 failure times with 0.50 days
    E1684$failtime[E1684$failtime == 0] = 0.50/365.25
    E1690$failtime[E1690$failtime == 0] = 0.50/365.25
    E1684$cage = as.numeric(scale(E1684$age))
    E1690$cage = as.numeric(scale(E1690$age))
    data_list = list(currdata = E1690, histdata = E1684)
    nbreaks = 3
    probs   = 1:nbreaks / nbreaks
    breaks  = as.numeric(
      quantile(E1690[E1690$failcens==1, ]$failtime, probs = probs)
    )
    breaks  = c(0, breaks)
    breaks[length(breaks)] = max(10000, 1000 * breaks[length(breaks)])
    fit.pwe.pp = pwe.pp(
      formula = survival::Surv(failtime, failcens) ~ treatment + sex + cage + node_bin,
      data.list = data_list,
      breaks = breaks,
      a0 = 0.5,
      get.loglik = TRUE,
      chains = 1, iter_warmup = 1000, iter_sampling = 2000
    )
    fit.pwe.post = pwe.post(
      formula = survival::Surv(failtime, failcens) ~ treatment + sex + cage + node_bin,
      data.list = data_list,
      breaks = breaks,
      get.loglik = TRUE,
      chains = 1, iter_warmup = 1000, iter_sampling = 2000
    )
    fit.pwe.commensurate = pwe.commensurate(
      formula = survival::Surv(failtime, failcens) ~ treatment + sex + cage + node_bin,
      data.list = data_list,
      breaks = breaks,
      p.spike = 0.1,
      get.loglik = TRUE,
      chains = 1, iter_warmup = 1000, iter_sampling = 2000
    )
    fit.list = list(fit.pwe.post, fit.pwe.pp, fit.pwe.commensurate)
    samples.mtx = do.call(
     cbind, lapply(fit.list, function(d){
       as.numeric( d[["treatment"]] )
     })
    )
    wts = compute.ensemble.weights(
      fit.list = fit.list,
      type = "pseudobma+"
    )$weights
    sample.ensemble(
     wts = wts, samples.mtx = samples.mtx
    )
  }
}