pwe.logml.pp: Log marginal likelihood of a piecewise exponential (PWE)...

In hdbayes: Bayesian Analysis of Generalized Linear Models with Historical Data

Log marginal likelihood of a piecewise exponential (PWE) model under the power prior (PP)

Description

Uses Markov chain Monte Carlo (MCMC) and bridge sampling to estimate the logarithm of the marginal likelihood of a PWE model under the power prior (PP).

The arguments related to MCMC sampling are utilized to draw samples from the power prior (PP). These samples are then used to compute the logarithm of the normalizing constant of the PP using only historical data set.

Usage

pwe.logml.pp(
  post.samples,
  bridge.args = NULL,
  iter_warmup = 1000,
  iter_sampling = 1000,
  chains = 4,
  ...
)

Arguments

post.samples	output from pwe.pp() giving posterior samples of a PWE model under the power prior (PP), 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().
iter_warmup	number of warmup iterations to run per chain. Defaults to 1000. See the argument iter_warmup in sample() method in cmdstanr package.
iter_sampling	number of post-warmup iterations to run per chain. Defaults to 1000. See the argument iter_sampling in sample() method in cmdstanr package.
chains	number of Markov chains to run. Defaults to 4. See the argument chains in sample() method in cmdstanr package.
...	arguments passed to sample() method in cmdstanr package (e.g., seed, refresh, init).

Value

If the power prior parameter (a_0) is equal to zero, then the function will return the same result as the output from pwe.logml.post().

If the power prior parameters (a_0) is non-zero, the function will return a list with the following objects

model

"pwe_pp"

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 normalizing constant of the power prior (PP) using all data sets

bs.hist

an object of class bridge or bridge_list containing the output from using bridgesampling::bridge_sampler() to compute the logarithm of the normalizing constant of the PP using historical data set

min_ess_bulk

the minimum estimated bulk effective sample size of the MCMC sampling

max_Rhat

the maximum Rhat

References

Chen, M.-H. and Ibrahim, J. G. (2000). Power prior distributions for Regression Models. Statistical Science, 15(1).

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

if (instantiate::stan_cmdstan_exists()) {
  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)
    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)])
    d.pp = pwe.pp(
      formula = survival::Surv(failtime, failcens) ~ treatment + sex + cage + node_bin,
      data.list = data_list,
      breaks = breaks,
      a0 = 0.5,
      chains = 1, iter_warmup = 500, iter_sampling = 1000
    )
    pwe.logml.pp(
      post.samples = d.pp,
      bridge.args = list(silent = TRUE),
      chains = 1, iter_warmup = 1000, iter_sampling = 2000
    )
  }
}

hdbayes documentation built on Nov. 21, 2025, 1:07 a.m.