BayesPET: Bayesian Prediction of Event Times for Blinded Randomized Controlled Trials

Description

The package implements Bayesian models for predicting the calendar time at which a target number of events is reached in clinical trials. The methodology applies to both blinded and unblinded settings and jointly models enrollment, event-time, and censoring processes.

Details

The package provides tools for trial data simulation, model fitting using 'Stan' via the 'rstan' interface, and event time prediction under a wide range of trial designs, including varying sample sizes, enrollment patterns, treatment effects, and event or censoring time distributions. The package is intended to support interim monitoring, operational planning, and decision-making in clinical trial development. For details, please refer to Fu et al. (2025) doi:10.1002/sim.70310.

Author(s)

Xinyi He xinyi.he@uth.tmc.edu, Jingyan Fu, and Ying Yuan

References

Fu, J., Zhao, D., Skanji, D., Liu, H., Tang, R. S., and Yuan, Y. (2025). Bayesian Prediction of Event Times Using Mixture Model for Blinded Randomized Controlled Trials. Statistics in Medicine, 44(28–30), e70310. doi:10.1002/sim.70310

Stan Development Team (2023). RStan: The R Interface to Stan. R package version 2.32.7. https://mc-stan.org

Solve baseline survival parameters by matching the marginal median survival time

Description

In practice, it is sometimes more convenient to work with a marginal (population-level) median survival time. This function solves for exactly one unknown baseline quantity - shape, scale, or median - in a proportional hazards (PH) survival model with either a Weibull or log-logistic baseline distribution, given the other two.

When covariates are provided, median is interpreted as the marginal (population-level) median survival time defined implicitly by

\mathbb{E}_{\boldsymbol{Z}}\{S(\mathrm{median}\mid \boldsymbol{Z})\} = 0.5

and the unknown quantity is obtained by Monte Carlo integration over the covariate distribution.

Usage

convert_median(
  distribution = "Weibull",
  shape,
  median,
  scale,
  cov_type = NULL,
  cov_dist = NULL,
  beta = NULL,
  S = 20000,
  seed.solveparam = 123,
  interval = if (is.null(scale)) c(1e-15, 100) else c(1e-06, 20000),
  tol = .Machine$double.eps^0.25,
  maxiter = 1000
)

Arguments

Details

The proportional hazards (PH) model assumes that, conditional on a covariate vector \boldsymbol{Z}, the hazard function satisfies

h(t \mid \boldsymbol{Z}) = h_0(t)\exp(\boldsymbol{Z}^\top \boldsymbol{\beta})

for t \ge 0, where h_0(t) is the baseline hazard and \boldsymbol{\beta} is a vector of regression coefficients.

The baseline model is specified by distribution through the baseline survival function S_0(t):

• "Weibull": baseline survival is

  S_0(t) = \exp\{-\lambda_0 t^{\rho}\}

  for t \ge 0, where \rho > 0 is the shape parameter and \lambda_0 > 0 is the scale parameter. An equivalent representation used by several standard implementations is

  S_0(t) = \exp\{-(t/\sigma_0)^{\rho}\}

  for t \ge 0, where \sigma_0 > 0 is a reparameterization satisfying \lambda_0 = \sigma_0^{-\rho}. This latter form is used by dweibull, phreg, flexsurvreg, and rstan.

• "Loglogistic": baseline survival is

  S_0(t) = \{1 + (t/b)^{a}\}^{-1}

  for t \ge 0, where a>0 is the shape parameter and b>0 is the scale parameter.

Under the PH assumption, the conditional survival satisfies

S(t \mid \boldsymbol{Z}) = S_0(t)^{\exp(\boldsymbol{Z}^\top \boldsymbol{\beta})}.

When no covariates are supplied (cov_type = NULL), the median survival time corresponds to the baseline median and closed-form solutions are available for some distributions.

When covariates are supplied, the median survival time m is interpreted as the marginal (population) median defined implicitly by

\mathbb{E}_{\boldsymbol{Z}}\{S(m \mid \boldsymbol{Z})\} = 0.5,

where the expectation is taken with respect to the covariate distribution implied by cov_type and cov_dist. This expectation is approximated using Monte Carlo simulation with S draws, and the unknown parameter is obtained by numerical root finding via uniroot.

Value

A numeric scalar giving the solved parameter:

• the shape parameter if shape = NULL;

• the baseline scale parameter if scale = NULL;

• the median survival time if median = NULL.

See Also

dweibull rllogis

Examples

# Weibull: convert a desired median to the corresponding scale parameter
## No covariates
convert_median(
  distribution = "Weibull",
  shape = 3,
  median = 5,
  scale = NULL,
  seed = 123
)
## With covariates
convert_median(
  distribution = "Weibull",
  shape = 3,
  median = 5,
  scale = NULL,
  cov_type = c("binary", "continuous"),
  cov_dist = c(0.5, 1),
  beta = c(1, 0.5),
  seed = 123
)

# Log-logistic: convert median to scale when covariates enter the model
convert_median(
  distribution = "Loglogistic",
  shape = 1.5,
  median = 3,
  scale = NULL,
  cov_type = c("binary", "continuous"),
  cov_dist = c(0.5, 1),
  beta = c(1, 0.5),
  seed = 123
)

Example trial datasets for fitting Stan models and predicting event times

Description

data_example is a simulated set of trial datasets designed to illustrate a complete analysis workflow, including fitting enrollment and survival and censoring models using interim data, and predicting the calendar time at which a prespecified number of events is expected to be observed.

Context. Consider a two-arm time-to-event trial that plans to enroll N = 200 participants, with equal randomization between the experimental and control arms. An interim analysis is conducted after E_\mathrm{cutoff} = 75 events have occurred, and the final analysis is planned at E_\mathrm{target} = 150 events. Based on the data available at the interim analysis, the goal is to predict the calendar time (measured from the trial start) at which the cumulative number of events will reach E_\mathrm{target}.

The bundle contains:

example_enroll

Enrollment information (interarrival times and enrollment status).

example_eventcensor

Treatment assignment, follow-up time, event indicator, censoring indicator, and covariates.

Usage

data(data_example)

Format

A named list with two data frames:

example_enroll

A data frame with columns:

No

Integer subject identifier.

interarrivaltime

Interarrival time.

enrollstatus

Enrollment status indicator (1 = enrolled, 0 = administratively censored).

example_eventcensor

A data frame with columns:

No

Integer subject identifier.

trt

Treatment assignment indicator (1 = experimental, 0 = control). This column is used for unblinded analyses and can be ignored or set to NA when fitting blinded event time models.

time

Observed follow-up time (event or censoring time, positive).

eventstatus

Event indicator (1 = event, 0 = right-censored).

censorstatus

Random censoring indicator (1 = random censoring observed, 0 = no random censoring, including administrative censoring or observed event).

X1

Binary covariate (0/1).

X2

Continuous covariate.

Source

Simulated data.

Examples

data(data_example)
names(data_example)
head(data_example$example_enroll)
head(data_example$example_eventcensor)

Fit a Weibull model for random censoring times

Description

Fits a Weibull proportional hazards model for random censoring times and returns posterior draws of the Weibull shape parameter \rho_c>0, the baseline scale parameter \lambda_{0c}>0, and, when covariates are included, the covariate log hazard ratios \boldsymbol{\beta}_c.

Usage

fit_censor(
  t_obs,
  status_censor,
  cov = NULL,
  hyperparams_censor = list(),
  chains = 4,
  iter = 4000,
  seed = 2,
  refresh = 0,
  mc.cores = 1,
  warmup = floor(iter/2),
  control = list(adapt_delta = 0.95),
  return_fit = FALSE,
  quiet = TRUE
)

Arguments

Details

Let T_{ci} denote the censoring time for subject i. Let \boldsymbol{Z}_i denote the covariate vector for subject i (the i-th row of cov when provided). The parameterization has baseline density function

f_{0c}(t) = \lambda_{0c}\,\rho_c\, t^{\rho_c-1}\exp(-\lambda_{0c}t^{\rho_c}), \quad t \ge 0.

The model assumes the survival function

S_c(t \mid \boldsymbol{Z}_i) = \exp\{-\lambda_{0c} \, t^{\rho_c} \exp(\boldsymbol{Z}_i^\top \boldsymbol{\beta}_c)\}, \quad t \ge 0.

The corresponding hazard function is

h_c(t \mid \boldsymbol{Z}_i) = \lambda_{0c}\,\rho_c\,t^{\rho_c-1}\exp(\boldsymbol{Z}_i^\top \boldsymbol{\beta}_c), \quad t \ge 0.

Note: In 'Stan', the Weibull distribution is parameterized by shape parameter \rho_c and scale parameter \sigma, with baseline density function

f_{0c}(t) = (\rho_c/\sigma)\, (t/\sigma)^{\rho_c-1}\, \exp(-(t/\sigma)^{\rho_c}), \quad t \ge 0.

To match the hazard-scale parameterization above, we set

\sigma_i = \{\lambda_{0c}\exp(\boldsymbol{Z}_i^\top \boldsymbol{\beta}_c)\}^{-1/\rho_c}.

The random censoring indicator status_censor follows the convention:

• 0: random censoring time observed at t_obs (density contribution),

• 1: administratively right-censored at t_obs (survival contribution).

Value

A list with the following components:

• rho_c: Posterior draws of the Weibull shape parameter \rho_c.

• lambda_0c: Posterior draws of the baseline Weibull scale parameter \lambda_{0c}.

• beta_c: Posterior draws of the covariate log hazard ratios \boldsymbol{\beta}_c, or NULL if no covariates are included.

• fit: The 'rstan' stanfit object (only if return_fit = TRUE).

See Also

Other BayesPET model fitting: fit_enroll, fit_event_blind, fit_event_unblind, fit_models,
print.BayesPET_fit

Examples

data(data_example)
example_eventcensor<-data_example$example_eventcensor

## ---- fit censoring model ----
## Reduced number of chains and iterations compared to defaults
## to keep the example computationally manageable.
fit <- fit_censor(
  t_obs = example_eventcensor$time,
  status_censor = example_eventcensor$censorstatus,
  cov = example_eventcensor[,6:7],
  chains = 2,
  iter = 2000, quiet = FALSE,
  seed = 2, return_fit = TRUE
)

summary(fit$rho_c)
summary(fit$lambda_0c)
summary(fit$beta_c)
print(fit$fit)

Fit enrollment model

Description

Fits an exponential model to enrollment interarrival times and returns posterior draws of the enrollment rate \mu.

Usage

fit_enroll(
  status_enroll,
  t_enroll,
  hyperparams_enroll = list(),
  chains = 4,
  iter = 4000,
  seed = 1,
  refresh = 0,
  warmup = floor(iter/2),
  mc.cores = 1,
  control = list(adapt_delta = 0.95),
  return_fit = FALSE,
  quiet = TRUE
)

Arguments

Details

Let t_i denote the interarrival time between consecutive enrollments. The enrollment process is modeled using an exponential distribution with rate parameter \mu, so that

t_i \sim \mathrm{Exponential}(\mu),

where \mu > 0 represents the average enrollment rate. The corresponding density is

f(t) = \mu e^{-\mu t}, \quad t \ge 0.

Administrative censoring of enrollment is handled through the indicator status_enroll, which follows the convention:

• 1: the interarrival time is fully observed at t_enroll (density contribution);

• 0: the interarrival time is administratively censored at t_enroll (survival contribution).

Value

A list with components:

• mu: posterior draws of \mu,

• fit: The 'rstan' stanfit fit object (only if return_fit = TRUE).

See Also

Other BayesPET model fitting: fit_censor, fit_event_blind, fit_event_unblind, fit_models,
print.BayesPET_fit

Examples

data(data_example)
example_enroll<-data_example$example_enroll

out <- fit_enroll(
  status_enroll = example_enroll$enrollstatus,
  t_enroll = example_enroll$interarrivaltime,
  seed = 1, return_fit = TRUE,
  quiet = FALSE
)

summary(out$mu)
print(out$fit)

Fit a Weibull event-time model with unknown treatment assignments

Description

Fits a Weibull event-time model formulated as a mixture over unobserved treatment assignments. The function returns posterior draws of the Weibull shape parameter \rho_e, the baseline scale parameter \lambda_{0e}, the treatment log hazard ratio \eta, and when covariates are included, the covariate log hazard ratios \boldsymbol{\beta}_e. Posterior draws of the latent treatment indicators x_i are also returned.

Usage

fit_event_blind(
  t_event,
  status_event,
  p_trt,
  cov = NULL,
  hyperparams_event = list(),
  chains = 4,
  iter = 4000,
  seed = 123,
  refresh = 0,
  warmup = floor(iter/2),
  mc.cores = 1,
  control = list(adapt_delta = 0.95),
  return_fit = FALSE,
  quiet = TRUE
)

Arguments

Details

Let T_{ei} denote the event time for subject i. Let \boldsymbol{Z}_i denote the corresponding covariate vector, which is the i th row of cov when provided. Treatment assignment is represented by a latent indicator x_i, where x_i = 0 denotes control and x_i = 1 denotes experimental treatment. The latent treatment indicators are jointly inferred with the model parameters from the observed event time data.

The treatment effect parameter \eta represents the log hazard ratio comparing experimental treatment to control, and \boldsymbol{\beta}_e denotes the vector of covariate regression coefficients (log hazard ratios) in the proportional hazards model.

The baseline event time distribution follows a Weibull model with density

f_{0e}(t) = \lambda_{0e}\,\rho_e\, t^{\rho_e - 1} \exp(-\lambda_{0e} t^{\rho_e}), \quad t \ge 0.

Conditional on treatment and covariates, the hazard function is

h_e(t \mid x_i, \boldsymbol{Z}_i) = \lambda_{0e}\,\rho_e\, t^{\rho_e - 1} \exp(\eta x_i + \boldsymbol{Z}_i^\top \boldsymbol{\beta}_e), \quad t \ge 0,

and the corresponding survival function is

S_e(t \mid x_i, \boldsymbol{Z}_i) = \exp\{-\lambda_{0e}\, t^{\rho_e} \exp(\eta x_i + \boldsymbol{Z}_i^\top \boldsymbol{\beta}_e)\}, \quad t \ge 0.

For reference, 'Stan' uses an equivalent Weibull representation based on a shape parameter \rho_e and a scale parameter \sigma, with baseline density

f_{0e}(t) = (\rho_e / \sigma)\, (t / \sigma)^{\rho_e - 1} \exp\{-(t / \sigma)^{\rho_e}\}, \quad t \ge 0.

The two formulations are related through

\sigma_i = \{\lambda_{0e}\exp(\eta x_i + \boldsymbol{Z}_i^\top \boldsymbol{\beta}_e)\}^{-1/\rho_e}.

Because x_i is not observed in blinded randomized trials, it is treated as a latent variable, and the observed event-time data marginally follow a mixture of two Weibull distributions corresponding to the latent treatment groups.

To avoid label switching in posterior inference, the treatment effect parameter \eta is assigned a Normal prior truncated to [0, \infty), restricting \eta to be nonnegative.

Value

A list with the following components:

• eta: Posterior draws of the treatment log hazard ratio \eta.

• rho_e: Posterior draws of the Weibull shape parameter \rho_e.

• lambda_0e: Posterior draws of the baseline Weibull scale parameter \lambda_{0e}.

• beta_e: Posterior draws of the covariate log hazard ratios \boldsymbol{\beta}_e, or NULL if no covariates are included.

• x: Posterior draws of latent treatment indicators x_i.

• fit: The 'rstan' stanfit object (only if return_fit = TRUE).

See Also

Other BayesPET model fitting: fit_censor, fit_enroll, fit_event_unblind, fit_models,
print.BayesPET_fit

Examples

data(data_example)
example_eventcensor<-data_example$example_eventcensor
# Use 2 chains and iter = 1000 here to reduce runtime for the example;
# use more chains in real analyses.
fit_e_blind <- fit_event_blind(
  t_event = example_eventcensor$time,
  status_event = example_eventcensor$eventstatus,
  cov = example_eventcensor[,6:7],
  p_trt = 0.5,
  chains = 2,
  iter = 1000, seed = 123,
  return_fit = TRUE
)

summary(fit_e_blind$eta)
print(fit_e_blind$fit)

Fit a Weibull event-time model with known treatment assignments

Description

Fits a Weibull event time model in which treatment assignments are observed. The function returns posterior draws of the Weibull shape parameter \rho_e > 0, the baseline scale parameter \lambda_{0e} > 0, the treatment effect coefficient (log hazard ratio) \eta, and when covariates are included, the covariate log hazard ratios \boldsymbol{\beta}_e.

Usage

fit_event_unblind(
  t_event,
  status_event,
  treatment_ind,
  cov = NULL,
  hyperparams_event = list(),
  chains = 4,
  iter = 4000,
  seed = 123,
  refresh = 0,
  warmup = floor(iter/2),
  mc.cores = 1,
  control = list(adapt_delta = 0.95),
  return_fit = FALSE,
  quiet = TRUE
)

Arguments

Details

Let T_{ei} denote the event time for subject i. Let \boldsymbol{Z}_i denote the corresponding covariate vector, which is the i th row of cov when provided. Treatment assignment is represented by a known indicator x_i, where for example x_i = 0 denotes control and x_i = 1 denotes experimental treatment. The treatment effect parameter \eta represents the log hazard ratio comparing experimental treatment to control, and \boldsymbol{\beta}_e denotes the vector of covariate regression coefficients (log hazard ratios) in the proportional hazards model.

The baseline event time distribution follows a Weibull model with density

f_{0e}(t) = \lambda_{0e}\,\rho_e\, t^{\rho_e - 1} \exp(-\lambda_{0e} t^{\rho_e}), \quad t \ge 0.

Conditional on treatment and covariates, the hazard function is

h_e(t \mid x_i, \boldsymbol{Z}_i) = \lambda_{0e}\,\rho_e\, t^{\rho_e - 1} \exp(\eta x_i + \boldsymbol{Z}_i^\top \boldsymbol{\beta}_e), \quad t \ge 0,

and the corresponding survival function is

S_e(t \mid x_i, \boldsymbol{Z}_i) = \exp\{-\lambda_{0e}\, t^{\rho_e} \exp(\eta x_i + \boldsymbol{Z}_i^\top \boldsymbol{\beta}_e)\}, \quad t \ge 0.

For reference, 'Stan' uses an equivalent Weibull representation based on a shape parameter \rho_e and a scale parameter \sigma, with baseline density

f_{0e}(t) = (\rho_e / \sigma)\, (t / \sigma)^{\rho_e - 1} \exp\{-(t / \sigma)^{\rho_e}\}, \quad t \ge 0.

The two formulations are related through

\sigma_i = \{\lambda_{0e}\exp(\eta x_i + \boldsymbol{Z}_i^\top \boldsymbol{\beta}_e)\}^{-1/\rho_e}, \quad t \ge 0.

Value

A list with the following components:

• eta: Posterior draws of the treatment log hazard ratio \eta.

• rho_e: Posterior draws of the Weibull shape parameter \rho_e.

• lambda_0e: Posterior draws of the baseline Weibull scale parameter \lambda_{0e}.

• beta_e: Posterior draws of the covariate log hazard ratios \boldsymbol{\beta}_e, or NULL if no covariates are included.

• fit: The 'rstan' stanfit object (only if return_fit = TRUE).

See Also

Other BayesPET model fitting: fit_censor, fit_enroll, fit_event_blind, fit_models,
print.BayesPET_fit

Examples

data(data_example)
example_eventcensor<-data_example$example_eventcensor
# Use chains = 1 here to reduce runtime for the example;
# use more chains in real analyses.
fit_e_unblind <- fit_event_unblind(
  t_event = example_eventcensor$time,
  status_event = example_eventcensor$eventstatus,
  treatment_ind = example_eventcensor$trt,
  cov = example_eventcensor[,6:7],
  chains = 1, iter = 2000, seed = 123,
  return_fit = TRUE
)

summary(fit_e_unblind$eta)
print(fit_e_unblind$fit)

Fit enrollment, event-time, and censoring models to clinical trial data and return posterior draws model parameters

Description

Fits the enrollment, event-time, and censoring models to trial data and returns posterior draws of model parameters.

Usage

fit_models(
  data.enroll,
  data.eventcensor,
  blinded = TRUE,
  p_trt = NULL,
  hyperparams_enroll = list(),
  hyperparams_event = list(),
  hyperparams_censor = list(),
  chains = 4,
  iter = 4000,
  mc.cores = 1,
  warmup = floor(iter/2),
  seed = list(123),
  refresh = 0,
  control = list(list(adapt_delta = 0.95)),
  return_fit = FALSE,
  quiet = TRUE
)

Arguments

Details

This function fits three submodels: an enrollment model, a censoring model, and an event-time model conditional on the given trial data. If treatment assignments are known, the event-time model is fit using fit_event_unblind; otherwise, a blinded event-time model is fit using
fit_event_blind. Technical details of the likelihoods, priors and parameterizations are documented in fit_enroll, fit_censor, fit_event_unblind, and fit_event_blind.

Value

An object of class "BayesPET_fit", a named list containing posterior draws and related information from the fitted models, with elements:

• blinded: Logical; indicates whether the analysis is blinded,

• mu: Posterior draws of the enrollment rate \mu.

• rho_c: Posterior draws of the censoring-model Weibull shape parameter \rho_c.

• lambda_0c: Posterior draws of the censoring-model baseline Weibull scale parameter \lambda_{0c}.

• beta_c: Posterior draws of the censoring-model covariate log hazard ratios \boldsymbol{\beta}_c, or NULL if no covariates are included.

• eta: Posterior draws of the treatment log hazard ratio \eta.

• rho_e: Posterior draws of the event-model Weibull shape parameter \rho_e.

• lambda_0e: Posterior draws of the event-model baseline Weibull scale parameter \lambda_{0e}.

• beta_e: Posterior draws of the event-model covariate log hazard ratios \boldsymbol{\beta}_e, or NULL if no covariates are included.

• treatment_ind: Observed treatment assignments x_i \in \{0,1\} (only returned if blinded = FALSE).

• x: Posterior draws of latent treatment assignments x_i \in \{0,1\} (only returned if blinded = TRUE).

• fit: A list with components enroll, censor, and event containing the underlying 'rstan' stanfit objects (present only if return_fit = TRUE).

The default print method displays a concise overview.

See Also

Other BayesPET model fitting: fit_censor, fit_enroll, fit_event_blind, fit_event_unblind,
print.BayesPET_fit

Examples

data(data_example)
example_enroll <- data_example$example_enroll
example_eventcensor <- data_example$example_eventcensor

# Unblinded analysis
## Use 2 chains and iter = 2000 here to reduce runtime for the example;
## use more chains in real analyses.
fit.unblind <- fit_models(
  data.enroll = example_enroll,
  data.eventcensor = example_eventcensor,
  blinded = FALSE,
  chains = 2, iter = 2000, seed = list(123),
  return_fit = TRUE, mc.cores = 1, quiet = FALSE
)

# Blinded analysis
example_eventcensor.blind <- example_eventcensor
example_eventcensor.blind$trt <- NA
## Use 2 chains and iter = 2000 here to reduce runtime for the example;
## use more chains in real analyses.
fit.blind <- fit_models(
  data.enroll = example_enroll,
  data.eventcensor = example_eventcensor.blind,
  blinded = TRUE, p_trt = 0.5,
  chains = 2, iter = 2000, seed = list(123),
  return_fit = TRUE, mc.cores = 1, quiet = FALSE
)

print(fit.unblind)
summary(fit.unblind$eta)

print(fit.blind)
summary(fit.blind$eta)

Generate two-arm trial data with enrollment, event, and censoring processes, and return data formatted for event-time prediction.

Description

Simulates data from a two-arm clinical trial with a time-to-event endpoint. The data generating process incorporates staggered enrollment, event times, and random censoring, with event and censoring distributions specified as Weibull or log-logistic. Treatment and covariate effects are incorporated through a proportional hazards structure.

Usage

generate_data(
  N,
  E_target,
  E_cutoff,
  p_trt,
  cov_type,
  cov_dist,
  logHR.trt = NULL,
  enroll_rate,
  dist.event,
  dist.censor,
  blinded = TRUE,
  event.scale = NULL,
  event.shape = NULL,
  censor.scale = NULL,
  censor.shape = NULL,
  beta.event,
  beta.censor,
  event.scale_trt = NULL,
  event.shape_trt = NULL,
  beta.event_trt = if (is.null(logHR.trt)) beta.event else NULL,
  assess_window = 0,
  seed = 123
)

Arguments

Details

Subjects are randomized independently to the experimental arm with probability p_trt. Baseline covariates are generated independently based on cov_type and cov_dist. Binary covariates follow a Bernoulli distribution, while continuous covariates follow a normal distribution with mean zero and standard deviation determined by cov_dist.

Interarrival times between successive enrollments are drawn from an exponential distribution with rate enroll_rate with model details documented in fit_enroll. Calendar enrollment times are obtained by cumulative summation of these interarrival times.

Event times and random censoring times are generated from Weibull or log-logistic baseline distributions, as specified by dist.event and dist.censor. For the Weibull model, the baseline survival function is parameterized as

S_0(t) = \exp\{-\lambda_0 t^{\rho}\}, \quad t \ge 0,

where \rho > 0 is the shape and \lambda_0 > 0 is the baseline hazard scale. For the log-logistic model, the baseline survival function is

S_0(t) = \{1 + (t / b)^a\}^{-1}, \quad t \ge 0,

where a > 0 and b > 0 denote the shape and scale parameters, respectively. Parameter calibration via marginal median survival can be performed using convert_median prior to simulation.

Covariate effects are incorporated through a proportional hazards structure for both the event and censoring processes. When logHR.trt is provided, the treatment effect is modeled through a proportional hazards formulation. When logHR.trt is NULL, the two treatment arms are allowed to differ through separate baseline parameters and covariate effects. The random censoring mechanism does not depend on treatment assignment.

Value

A list with elements:

• data.enroll: A data frame of observed enrollment information up to the interim data cut. Columns: subject index No, subject enrollment calendar time enrolltime, enrollment interarrival time interarrivaltime, enrollment status enrollstatus (1 = enrolled, 0 = administratively censored enrollment process).

• data.eventcensor: A data frame of observed survival outcomes at the interim cut. Columns include subject index No, treatment assignment indicator trt (NA if blinded = TRUE), observed time time (administratively censored at interim cut), event status eventstatus (1 = event, 0 = right censored), random censoring status censorstatus (1 = random censoring before the interim data cut; 0 = otherwise), followed by covariates.

• truesurvival: Full underlying data without administrative censoring: No, trt, t_event (true underlying event time), t_randcensor (true underlying random censoring time), t_event.obs (underlying follow-up time without administrative censoring),
  t_event.obswithintervalassess (underlying follow-up time without administrative censoring but applied with assessment windows), status (1 = event before random censoring), enrollmenttime, plus covariates.

• event.interim.obs: Observed number of events at the interim data cut. If the prespecified interim event cutoff E_cutoff cannot be reached, this equals the maximum number of events observed.

• event.max: Number of events that would occur without administrative censoring (i.e., after accounting only for random censoring).

• cuttime.true: The true calendar time at which the cumulative number of observed events reaches the target event count. If the target number of events cannot be reached, this is the calendar time of the last observed event or censoring.

• event.final.obs: The latent true number of events at cuttime.true.

Examples

## --- Weibull event/censoring with a common PH treatment effect ---
data.weibull <- generate_data(
  N = 80, E_target = 50, E_cutoff = 25, p_trt = 0.5,
  cov_type = c("binary", "continuous"),
  cov_dist = c(0.5, sqrt(2)),
  beta.event  = c(0.2, 0.2),
  beta.censor = c(0, 0),
  logHR.trt = log(0.5),
  enroll_rate = 50/3, beta.event_trt = NULL,
  dist.event = "Weibull", dist.censor = "Weibull",
  event.scale = 1/5^3, event.shape = 3,
  censor.scale = 10^(-6), censor.shape = 6,
  blinded = TRUE,
  assess_window = 2,
  seed = 1
)
names(data.weibull)

## --- Log-logistic event/censoring with a common PH treatment effect ---
data.logl <- generate_data(
  N = 80, E_target = 50, E_cutoff = 25, p_trt = 0.5,
  cov_type = c("binary", "continuous"),
  cov_dist = c(0.5, sqrt(2)),
  beta.event  = c(0.2, 0.2),
  beta.censor = c(0, 0),
  logHR.trt = log(0.5),
  enroll_rate = 50/3, beta.event_trt = NULL,
  dist.event = "Loglogistic", dist.censor = "Loglogistic",
  event.scale = 6, event.shape = 6,
  censor.scale = 20, censor.shape = 4,
  blinded = TRUE,
  assess_window = 2,
  seed = 1
)
summary(data.logl$truesurvival$t_event)
### true underlying event time without administrative censoring

## --- Weibull arm-specific models (logHR.trt = NULL) ---
data.weibull.nonPH <- generate_data(
  N = 80, E_target = 50, E_cutoff = 25, p_trt = 0.5,
  cov_type = c("binary", "continuous"),
  cov_dist = c(0.5, sqrt(2)),
  beta.event  = c(0.2, 0.2),
  beta.censor = c(0, 0),
  logHR.trt = NULL,
  enroll_rate = 50/3,
  dist.event = "Weibull", dist.censor = "Weibull",
  event.scale = 1/5^3, event.shape = 3,     # control
  event.scale_trt = 1/6^3, event.shape_trt = 3, # experiment
  beta.event_trt = c(0.15, 0.2),
  censor.scale = 10^(-6), censor.shape = 6,
  blinded = TRUE,
  assess_window = 2,
  seed = 1
)
data.weibull.nonPH$cuttime.true

Generate operating characteristics for event prediction

Description

Generate operating characteristics for multiple simulated trials

Usage

get_oc(
  N,
  E_target,
  E_cutoff,
  p_trt,
  cov_type,
  cov_dist,
  logHR.trt = NULL,
  enroll_rate,
  dist.event,
  dist.censor,
  blinded = TRUE,
  event.scale = NULL,
  event.shape = NULL,
  censor.scale = NULL,
  censor.shape = NULL,
  beta.event,
  beta.censor,
  event.scale_trt = NULL,
  event.shape_trt = NULL,
  beta.event_trt = if (is.null(logHR.trt)) beta.event else NULL,
  assess_window = 0,
  seed = 123,
  hyperparams_enroll = list(),
  hyperparams_event = list(),
  hyperparams_censor = list(),
  chains = 4,
  iter = 4000,
  warmup = floor(iter/2),
  refresh = 0,
  control = list(list(adapt_delta = 0.95)),
  nsim = 1000,
  nsim.max = 3 * nsim,
  n_workers = 1,
  ...
)

Arguments

Details

This function first simulates a two-arm time-to-event trial using generate_data and constructs interim datasets (data.enroll and data.eventcensor). It then fits the enrollment, event-time, and random censoring models and generates posterior predictive draws of the calendar time at which the cumulative number of events reaches E_target using predict_eventtime.

Only simulated datasets in which the target number of events E_target is reachable are retained for analysis. Up to nsim.max datasets are generated in order to obtain at most nsim valid replicates. As a result, the retained replicates correspond to a subset of the attempted simulations, and their associated seeds are recorded explicitly in the returned object.

Value

An object of class "BayesPET_oc" containing operating characteristics for event-time prediction, based on simulated trial replicates for which the target number of events E_target is reachable. The object is a named list with components:

• replicate: A data frame with one row per retained (valid) simulated trial replicate. Columns include:

  • median: Median of the posterior predictive draws of the calendar time at which E_target events are reached.

  • cuttime.true: True calendar time at which the cumulative number of observed events in the simulated trial first reaches E_target.

  • difference: Absolute prediction error, abs(median - cuttime.true).

• n_valid: Number of retained replicates.

• n_attempt: Total number of simulated datasets attempted.

• nsim_target: Target number of replicates requested.

• nsim_max: Maximum number of datasets that may be generated.

• seed0: Base random seed used to generate simulated datasets; equals the input argument seed.

• seeds: Integer vector of seeds corresponding to the retained replicates (one per row of replicate).

• call: Matched function call.

See Also

Other BayesPET operating characteristics: summary.BayesPET_oc()

Examples

## Using nsim = 2, chains = 2, and iter = 2000 to reduce runtime.
## Use larger nsim, chains and iter in real analyses.
oc <- get_oc(
  N = 200, E_target = 150,
  E_cutoff = 75, p_trt = 0.5,
  cov_type = c("binary", "continuous"),
  cov_dist = c(0.5, 2),
  beta.event = c(-0.2, -0.2),
  beta.censor = c(0, 0),
  logHR.trt = log(0.65),
  enroll_rate = 16,
  dist.event = "Weibull", dist.censor = "Weibull",
  event.scale = 1/5^3, event.shape = 3,
  censor.scale = 1/10^6, censor.shape = 6,
  blinded = TRUE,
  assess_window = 2,
  seed = 1,
  chains = 2, iter = 2000,
  nsim = 2,
  n_workers = 1
)

summary(oc)

Plot method for BayesPET prediction objects

Description

Plots an object of class "BayesPET_predtime" by displaying a histogram of posterior predictive draws of the calendar time at which the target number of events is reached. Only finite draws are included. A vertical line indicates the posterior median when it is finite.

Usage

## S3 method for class 'BayesPET_predtime'
plot(
  x,
  breaks = "Sturges",
  xlab = "Predicted calendar time to reach target number of events",
  ...
)

Arguments

Value

Invisibly returns NULL.

See Also

Other BayesPET prediction: predict_eventtime, print.BayesPET_predtime,
summary.BayesPET_predtime

Examples

data(data_example)
## Reduced number of chains and iterations compared to defaults
## to keep the example computationally manageable.
pred <- predict_eventtime(
  N = 200,
  E_target = 150,
  data.enroll = data_example$example_enroll,
  data.eventcensor = data_example$example_eventcensor,
  blinded = TRUE,
  p_trt = 0.5,
  chains = 2,
  iter = 2000,
  assess_window = 2,
  seed.fit = 1,
  seed.pred = 2,
  return_fit = TRUE,
  return_draws = TRUE,
  quiet = TRUE
)

print(pred)
summary(pred)
plot(pred)

Predict the calendar time at which a target number of events is reached from interim analysis data

Description

Fits enrollment, event-time, and random censoring models to data observed at the interim analysis and predicts the calendar time at which the cumulative number of events reaches E_target.

Usage

predict_eventtime(
  N,
  E_target,
  data.enroll,
  data.eventcensor,
  blinded = TRUE,
  p_trt = NULL,
  hyperparams_enroll = list(),
  hyperparams_event = list(),
  hyperparams_censor = list(),
  chains = 4,
  iter = 4000,
  warmup = floor(iter/2),
  seed.fit = list(123),
  refresh = 0,
  control = list(list(adapt_delta = 0.95)),
  mc.cores = 1,
  assess_window = 0,
  seed.pred = 1,
  return_fit = FALSE,
  quiet = TRUE,
  return_draws = FALSE
)

Arguments

Details

The function fits three components in sequence using fit_models: an enrollment model, an event-time model, and a random censoring model, all based on data observed at the interim analysis. It then performs posterior predictive simulation for two groups of subjects:

• \mathcal{B}_1: subjects enrolled before the interim analysis who have not experienced an event or random censoring by the interim data cut;

• \mathcal{B}_2: subjects not yet enrolled by the interim analysis.

For each posterior draw, the function estimates the calendar time at which the cumulative number of events reaches E_target.

Value

An object of class "BayesPET_predtime", which is a named list with components:

• prediction: A numeric vector of length S, where S is the number of posterior draws from the fitted models. Each element is a posterior predictive draw of the calendar time at which the cumulative number of events reaches E_target. Values may be Inf if fewer than E_target events occur in that draw.

• fit: Present only if return_fit = TRUE. A list of 'rstan' stanfit objects with components enroll, censor, and event, corresponding to the enrollment, censoring, and event-time models, respectively.

• draws: Present only if return_draws = TRUE. A list of posterior draws of model parameters produced by the fitted submodels. This includes posterior draws from the enrollment, event-time, and censoring models.

• call: The function call used to generate this object.

Methods include print, summary, and plot.

See Also

Other BayesPET prediction: plot.BayesPET_predtime, print.BayesPET_predtime,
summary.BayesPET_predtime

Examples

data(data_example)
## Reduced number of chains and iterations compared to defaults
## to keep the example computationally manageable.
pred <- predict_eventtime(
  N = 200,
  E_target = 150,
  data.enroll = data_example$example_enroll,
  data.eventcensor = data_example$example_eventcensor,
  blinded = TRUE,
  p_trt = 0.5,
  chains = 2,
  iter = 2000,
  assess_window = 2,
  seed.fit = 1,
  seed.pred = 2,
  return_fit = TRUE,
  return_draws = TRUE,
  quiet = TRUE
)

print(pred)
summary(pred)
plot(pred)

Print method for BayesPET model fitting objects

Description

Displays a concise overview of an object of class "BayesPET_fit", including whether the analysis is blinded, the number of posterior draws, and the model components.

Usage

## S3 method for class 'BayesPET_fit'
print(x, ...)

Arguments

Value

The object x, invisibly.

See Also

Other BayesPET model fitting: fit_censor, fit_enroll, fit_event_blind, fit_models,
fit_event_unblind

Examples

data(data_example)
example_enroll <- data_example$example_enroll
example_eventcensor <- data_example$example_eventcensor

# Blinded analysis
example_eventcensor.blind <- example_eventcensor
example_eventcensor.blind$trt <- NA
## Use 2 chains and iter = 2000 here to reduce runtime for the example;
## use more chains in real analyses.
fit.blind <- fit_models(
  data.enroll = example_enroll,
  data.eventcensor = example_eventcensor.blind,
  blinded = TRUE, p_trt = 0.5,
  chains = 2, iter = 2000, seed = list(123),
  return_fit = TRUE, mc.cores = 1, quiet = FALSE
)
print(fit.blind)

Print method for BayesPET prediction objects

Description

Displays a brief overview of an object of class "BayesPET_predtime" and lists the components available in the result. For numerical summaries, use summary; for visualization, use plot.

Usage

## S3 method for class 'BayesPET_predtime'
print(x, ...)

Arguments

Value

The object x, invisibly.

See Also

Other BayesPET prediction: plot.BayesPET_predtime, predict_eventtime,
summary.BayesPET_predtime

Examples

data(data_example)
## Reduced number of chains and iterations compared to defaults
## to keep the example computationally manageable.
pred <- predict_eventtime(
  N = 200,
  E_target = 150,
  data.enroll = data_example$example_enroll,
  data.eventcensor = data_example$example_eventcensor,
  blinded = TRUE,
  p_trt = 0.5,
  chains = 2,
  iter = 2000,
  assess_window = 2,
  seed.fit = 1,
  seed.pred = 2,
  return_fit = TRUE,
  return_draws = TRUE,
  quiet = TRUE
)

print(pred)
summary(pred)
plot(pred)

Summary method for BayesPET operating characteristics object

Description

Computes summary measures of prediction accuracy from a "BayesPET_oc" object. The summaries are based on the differences between the predicted median and true calendar times to reach the target number of events.

Usage

## S3 method for class 'BayesPET_oc'
summary(object, thresholds = c(0.25, 0.5, 1, 1.5, 2), ...)

## S3 method for class 'summary.BayesPET_oc'
print(x, digits = 3, ...)

Arguments

Value

summary.BayesPET_oc returns an object of class "summary.BayesPET_oc", a list containing:

• n_valid: Number of valid simulation replicates (where the target event count can be reached).

• n_attempt: Number of datasets generated to obtain the valid replicates.

• success_rate: Proportion of attempted datasets that produced valid replicates.

• thresholds: Threshold values used to evaluate prediction accuracy.

• pr_lt: Proportion of replicates with prediction error less than each threshold.

• mae: Mean absolute prediction error.

• median_ae: Median absolute prediction error.

• rmse: Root mean squared prediction error.

• call: The matched function call.

print.summary.BayesPET_oc returns the object x, invisibly.

See Also

Other BayesPET operating characteristics: get_oc()

Examples

## Using nsim = 2, chains = 2, and iter = 2000 to reduce runtime.
## Use larger nsim, chains and iter in real analyses.
oc <- get_oc(
  N = 200, E_target = 150,
  E_cutoff = 75, p_trt = 0.5,
  cov_type = c("binary", "continuous"),
  cov_dist = c(0.5, 2),
  beta.event = c(-0.2, -0.2),
  beta.censor = c(0, 0),
  logHR.trt = log(0.65),
  enroll_rate = 16,
  dist.event = "Weibull", dist.censor = "Weibull",
  event.scale = 1/5^3, event.shape = 3,
  censor.scale = 1/10^6, censor.shape = 6,
  blinded = TRUE,
  assess_window = 2,
  seed = 1,
  chains = 2, iter = 2000,
  nsim = 2,
  n_workers = 1
)

summary(oc)

Summary method for BayesPET prediction objects

Description

Summarizes an object of class "BayesPET_predtime" by reporting summary statistics of the posterior predictive distribution of the calendar time at which the target number of events is reached.

Usage

## S3 method for class 'BayesPET_predtime'
summary(object, ...)

## S3 method for class 'summary.BayesPET_predtime'
print(x, digits = 4, ...)

Arguments

Value

summary.BayesPET_predtime returns an object of class "summary.BayesPET_predtime", which is a named list containing summary information for the posterior predictive distribution of the target event time. Components include:

• S: Total number of posterior predictive draws.

• n_infinite: Number of draws in which the target number of events is not reached.

• prob_not_reached: Proportion of draws in which the target is not reached.

• q25: Posterior 25th percent quantile of the target event time.

• median: Posterior median of the target event time.

• q75: Posterior 75th percent quantile of the target event time.

• call: The function call used to generate the prediction.

print.summary.BayesPET_predtime returns the object x, invisibly.

See Also

Other BayesPET prediction: plot.BayesPET_predtime, predict_eventtime,
print.BayesPET_predtime

Examples

data(data_example)
## Reduced number of chains and iterations compared to defaults
## to keep the example computationally manageable.
pred <- predict_eventtime(
  N = 200,
  E_target = 150,
  data.enroll = data_example$example_enroll,
  data.eventcensor = data_example$example_eventcensor,
  blinded = TRUE,
  p_trt = 0.5,
  chains = 2,
  iter = 2000,
  assess_window = 2,
  seed.fit = 1,
  seed.pred = 2,
  return_fit = TRUE,
  return_draws = TRUE,
  quiet = TRUE
)

print(pred)
summary(pred)
plot(pred)