Learn R Programming

hdbayes (version 0.2.0)

aft.logml.npp: Log marginal likelihood of an accelerated failure time (AFT) model under normalized power prior (NPP)

Description

Uses bridge sampling to estimate the logarithm of the marginal likelihood of an AFT model under the normalized power prior (NPP).

Usage

aft.logml.npp(post.samples, bridge.args = NULL)

Value

The function returns a list with the following objects

model

"aft_npp"

logml

the estimated logarithm of the marginal likelihood

bs

an object of class bridge or bridge_list containing the output from using bridgesampling::bridge_sampler() to compute the logarithm of the marginal likelihood of the normalized power prior (NPP)

Arguments

post.samples

output from aft.npp() giving posterior samples of an AFT model under the normalized power prior (NPP), with an attribute called 'data' which includes the list of variables specified in the data block of the Stan program.

bridge.args

a list giving arguments (other than samples, log_posterior, data, lb, and ub) to pass onto bridgesampling::bridge_sampler().

References

Duan, Y., Ye, K., and Smith, E. P. (2005). Evaluating water quality using power priors to incorporate historical information. Environmetrics, 17(1), 95–106.

Gronau, Q. F., Singmann, H., and Wagenmakers, E.-J. (2020). bridgesampling: An r package for estimating normalizing constants. Journal of Statistical Software, 92(10).

Examples

Run this code
# \donttest{
  if(requireNamespace("parallel")){
    library(parallel)
    ncores    = 2

    if(requireNamespace("survival")){
      library(survival)
      data(E1684)
      data(E1690)
      ## take subset for speed purposes
      E1684 = E1684[1:100, ]
      E1690 = E1690[1:50, ]
      ## 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)
      formula = survival::Surv(failtime, failcens) ~ treatment + sex + cage + node_bin
    }

    a0 = seq(0, 1, length.out = 11)
    if (instantiate::stan_cmdstan_exists()) {
      ## call created function
      ## wrapper to obtain log normalizing constant in parallel package
      logncfun = function(a0, ...){
        hdbayes::aft.npp.lognc(
          formula = formula, histdata = data_list[[2]], a0 = a0, dist = "weibull",
          ...
        )
      }

      cl = makeCluster(ncores)
      clusterSetRNGStream(cl, 123)
      clusterExport(cl, varlist = c('formula', 'data_list'))
      a0.lognc = parLapply(
        cl = cl, X = a0, fun = logncfun, iter_warmup = 500,
        iter_sampling = 1000, chains = 1, refresh = 0
      )
      stopCluster(cl)
      a0.lognc = data.frame( do.call(rbind, a0.lognc) )

      ## sample from normalized power prior
      d.npp = aft.npp(
        formula = formula,
        data.list = data_list,
        a0.lognc = a0.lognc$a0,
        lognc = a0.lognc$lognc,
        dist = "weibull",
        chains = 1, iter_warmup = 500, iter_sampling = 1000,
        refresh = 0
      )
      aft.logml.npp(
        post.samples = d.npp,
        bridge.args = list(silent = TRUE)
      )
    }
  }
# }

Run the code above in your browser using DataLab