pwe.npp.lognc: Estimate the logarithm of the normalizing constant for...

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

Estimate the logarithm of the normalizing constant for normalized power prior (NPP) for one data set

Description

Uses Markov chain Monte Carlo (MCMC) and bridge sampling to estimate the logarithm of the normalizing constant of a piecewise exponential (PWE) model under the NPP for a fixed value of the power prior parameter a_0 \in (0, 1) for one data set. The initial priors are independent normal priors on the regression coefficients and half-normal priors on the baseline hazard parameters.

Usage

pwe.npp.lognc(
  formula,
  histdata,
  breaks,
  a0,
  beta.mean = NULL,
  beta.sd = NULL,
  base.hazard.mean = NULL,
  base.hazard.sd = NULL,
  bridge.args = NULL,
  iter_warmup = 1000,
  iter_sampling = 1000,
  chains = 4,
  ...
)

Arguments

formula	a two-sided formula giving the relationship between the response variable and covariates. The response is a survival object as returned by the survival::Surv(time, event) function, where event is a binary indicator for event (0 = no event, 1 = event has occurred). The type of censoring is assumed to be right-censoring.
histdata	a data.frame giving the historical data.
breaks	a numeric vector specifying the time points that define the boundaries of the piecewise intervals. The values should be in ascending order, with the final value being greater than or equal to the maximum observed time.
a0	a scalar between 0 and 1 giving the (fixed) power prior parameter for the historical data.
beta.mean	a scalar or a vector whose dimension is equal to the number of regression coefficients giving the mean parameters for the initial prior on regression coefficients. If a scalar is provided, beta.mean will be a vector of repeated elements of the given scalar. Defaults to a vector of 0s.
beta.sd	a scalar or a vector whose dimension is equal to the number of regression coefficients giving the sd parameters for the initial prior on regression coefficients. If a scalar is provided, same as for beta.mean. Defaults to a vector of 10s.
base.hazard.mean	a scalar or a vector whose dimension is equal to the number of intervals giving the location parameters for the half-normal priors on the baseline hazards of the PWE model. If a scalar is provided, same as for beta.mean. Defaults to 0.
base.hazard.sd	a scalar or a vector whose dimension is equal to the number of intervals giving the scale parameters for the half-normal priors on the baseline hazards of the PWE model. If a scalar is provided, same as for beta.mean. Defaults to 10.
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

The function returns a vector giving the value of a0, the estimated logarithm of the normalizing constant, the minimum estimated bulk effective sample size of the MCMC sampling, and the maximum Rhat.

References

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)
    ## take subset for speed purposes
    E1684 = E1684[1:100, ]
    ## replace 0 failure times with 0.50 days
    E1684$failtime[E1684$failtime == 0] = 0.50/365.25
    E1684$cage = as.numeric(scale(E1684$age))
    nbreaks = 3
    probs   = 1:nbreaks / nbreaks
    breaks  = as.numeric(
      quantile(E1684[E1684$failcens==1, ]$failtime, probs = probs)
    )
    breaks  = c(0, breaks)
    breaks[length(breaks)] = max(10000, 1000 * breaks[length(breaks)])
    pwe.npp.lognc(
      formula = survival::Surv(failtime, failcens) ~ treatment + sex + cage + node_bin,
      histdata = E1684,
      breaks = breaks,
      a0 = 0.5,
      bridge.args = list(silent = TRUE),
      chains = 1, iter_warmup = 500, iter_sampling = 1000
    )
  }
}

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