Bayesian binary outcome designs

library(goldilocks)

goldilocks can analyze completed binary endpoints using method = "bayes-bin". This is useful when the final decision is based on the event indicator by a fixed endpoint time, rather than the event time itself. The trial simulator still uses time-to-event hazards to model not-yet-observed outcomes at interim looks, but it can either impute a future event time or draw the binary endpoint status directly. The final analysis reduces each completed dataset to binary event status by end_of_study.

Two distinct priors can therefore enter a Bayesian binary design. The prior argument is a Gamma prior on the piecewise-exponential hazards. It is used only for the event-time model that generates predictive imputations at interim looks (and at the final analysis when imputed_final = TRUE). The bin_prior argument is a Beta prior on the binary endpoint event probability in each arm. It is used for the Bayesian binary analysis of every completed or imputed dataset. Thus, interim predictive probabilities are affected by both priors, whereas a non-imputed final binary analysis is affected by bin_prior alone.

For the examples below, we use the default weakly informative \(\operatorname{Gamma}(0.1, 0.1)\) hazard prior and a uniform \(\operatorname{Beta}(1, 1)\) binary-endpoint prior:

hazard_prior <- c(0.1, 0.1)  # Gamma shape and rate for predictive imputation
bin_prior <- c(1, 1)         # Beta shapes for the binary endpoint analysis

Two practical consequences follow:

method = "bayes-bin" only supports one-sided alternatives: alternative = "less" or alternative = "greater".
The final binary analysis requires complete endpoint status. Subjects who are censored before end_of_study must be imputed (imputed_final = TRUE) or, if imputed_final = FALSE, excluded from the final binary analysis.

Two-arm design

Suppose the control event probability by 12 months is 35%, and the treatment is expected to reduce this to 25%. We use a beta-binomial final analysis with a uniform \(\operatorname{Beta}(1, 1)\) prior in each arm. Since lower event probability is beneficial, we use alternative = "less" and compare the posterior distribution of

\[p_{\text{treatment}} - p_{\text{control}}\]

against h0 = 0.

end_of_study <- 12
hc <- prop_to_haz(0.35, endtime = end_of_study)
ht <- prop_to_haz(0.25, endtime = end_of_study)

two_arm_args <- list(
  hazard_treatment = ht,
  hazard_control = hc,
  cutpoints = NULL,
  N_total = 120,
  lambda = 10,
  lambda_time = NULL,
  interim_look = 80,
  end_of_study = end_of_study,
  prior = hazard_prior,
  bin_prior = bin_prior,
  bin_method = "quadrature",
  block = 2,
  rand_ratio = c(1, 1),
  prop_loss = 0,
  alternative = "less",
  h0 = 0,
  Fn = 0.05,
  Sn = 0.90,
  prob_ha = 0.95,
  N_impute = 20,
  method = "bayes-bin",
  imputed_final = FALSE
)

out_two_arm <- do.call(survival_adapt, two_arm_args)
out_two_arm
#>   prob_threshold margin alternative N_treatment N_control N_enrolled N_max
#> 1           0.95      0        less          60        60        120   120
#>   post_prob_ha  est_final ppp_success stop_futility stop_expected_success
#> 1    0.9178978 -0.1129032        0.85             0                     0

The output has the same structure as other survival_adapt() analyses. For method = "bayes-bin", est_final is the posterior mean binary effect: the treatment event probability minus the control event probability. The post_prob_ha value is the posterior probability that this difference is below h0 when alternative = "less".

Choosing the binary imputation approach

The default binary_imputation = "event-time" approach samples a future event time from the piecewise-exponential model, conditional on the subject remaining event-free through the available follow-up time \(T\). The sampled time is then converted to event or no event at the endpoint \(T^*\).

With binary_imputation = "bernoulli", the package skips the unused event time and calculates the endpoint probability directly:

\[ \begin{aligned} p &= \Pr(T_{\text{event}} \leq T^* \mid T_{\text{event}} > T) \\ &= \frac{S(T) - S(T^*)}{S(T)} \\ &= 1 - \exp\left\{-[H(T^*) - H(T)]\right\}. \end{aligned} \]

It then draws \(X \sim \operatorname{Bernoulli}(p)\). For subjects who are not yet enrolled, \(T=0\). Observed events are not imputed. Each completed dataset still uses a sampled posterior hazard, so both approaches retain uncertainty in the predictive piecewise-exponential model.

The approaches have the same endpoint distribution. They can nevertheless produce different Monte Carlo realizations if their random-number streams differ, so operating characteristics should be compared using enough trials and imputations. The following small seeded comparison uses the same design and random-number stream:

compare_binary_imputation <- function(imputation) {
  set.seed(2101)
  fit <- do.call(
    survival_adapt,
    c(two_arm_args, list(binary_imputation = imputation))
  )
  fit[c("ppp_success", "post_prob_ha", "est_final")]
}

rbind(
  `conditional event time` = compare_binary_imputation("event-time"),
  `direct Bernoulli status` = compare_binary_imputation("bernoulli")
)
#>                         ppp_success post_prob_ha  est_final
#> conditional event time            0    0.4042967 0.02380952
#> direct Bernoulli status           0    0.4042967 0.02380952

The direct Bernoulli calculation is also more stable in extreme tails because it evaluates the remaining cumulative hazard through -expm1(-(H(T^*) - H(T))), rather than subtracting two survival probabilities that may both round to zero.

Single-arm design

For a single-arm design, set hazard_control = NULL. The comparator is an external benchmark event probability supplied through h0, often called a performance goal (PG) or objective performance criterion (OPC).

Here the benchmark event probability is 30%, and the target event probability is 20%. With alternative = "less", success means the posterior probability that the event probability is below 30% exceeds prob_ha.

benchmark <- 0.30
target <- 0.20
ht_single <- prop_to_haz(target, endtime = end_of_study)

out_single_arm <- survival_adapt(
  hazard_treatment = ht_single,
  hazard_control = NULL,
  cutpoints = NULL,
  N_total = 80,
  lambda = 8,
  lambda_time = NULL,
  interim_look = 50,
  end_of_study = end_of_study,
  prior = hazard_prior,
  bin_prior = bin_prior,
  bin_method = "quadrature",
  prop_loss = 0,
  alternative = "less",
  h0 = benchmark,
  Fn = 0.05,
  Sn = 0.90,
  prob_ha = 0.95,
  N_impute = 20,
  method = "bayes-bin",
  imputed_final = FALSE
)

out_single_arm
#>   prob_threshold margin alternative N_treatment N_control N_enrolled N_max
#> 1           0.95    0.3        less          50         0         50    80
#>   post_prob_ha est_final ppp_success stop_futility stop_expected_success
#> 1    0.9980375 0.1346154           1             0                     1

In this setting est_final is the posterior mean event probability in the single arm, and post_prob_ha is the posterior probability that this event probability is below the benchmark.

Choosing `bin_method`

The posterior probability can be calculated in three ways:

bin_method = "mc" draws from the beta posterior directly. Use N_mcmc to control the number of Monte Carlo draws.
bin_method = "normal" uses a normal approximation to the posterior mean or treatment-control difference.
bin_method = "quadrature" uses numerical integration for the two-arm posterior difference, and the closed-form beta CDF for single-arm designs.

The Monte Carlo method has simulation error, controlled by N_mcmc. Quadrature is deterministic (up to numerical integration for a two-arm design) and is a useful high-accuracy default when it is computationally feasible. The normal approximation is fastest in a representative two-arm benchmark, but should be used with care when sample sizes are small or event probabilities are near 0 or 1.

Operating characteristics

As with the survival methods, operating characteristics should be evaluated by simulation. The seed argument makes the simulation reproducible, including when ncores is greater than 1. With backend = "auto", sim_trials() uses forked workers on Unix-like platforms and PSOCK workers on Windows.

out_power <- sim_trials(
  N_trials = 1000,
  hazard_treatment = ht,
  hazard_control = hc,
  cutpoints = NULL,
  N_total = 120,
  lambda = 10,
  lambda_time = NULL,
  interim_look = 80,
  end_of_study = end_of_study,
  prior = hazard_prior,
  bin_prior = bin_prior,
  bin_method = "quadrature",
  block = 2,
  rand_ratio = c(1, 1),
  prop_loss = 0,
  alternative = "less",
  h0 = 0,
  Fn = 0.05,
  Sn = 0.90,
  prob_ha = 0.95,
  N_impute = 20,
  method = "bayes-bin",
  imputed_final = FALSE,
  return_trace = TRUE,
  ncores = 2,
  seed = 5107
)

out_null <- update(out_power, hazard_treatment = hc, seed = 5108)

oc <- summarise_sims(list(
  "target: treatment event probability 25%" = out_power$sims,
  "null: treatment event probability 35%" = out_null$sims
))
oc$true_treatment_event_probability <- c(0.25, 0.35)

oc
plot_sim_ocs(
  oc,
  effect = "true_treatment_event_probability",
  xlab = "True treatment event probability"
)
plot_sim_stopping(out_power)
plot_sim_decisions(out_power)

These plots separate three complementary questions: how operating characteristics change with the true binary event probability, where enrollment stops within one scenario, and how the pair of predictive probabilities drives each interim decision. Set return_trace = TRUE only for scenarios where the decision map is needed; it does not change the compact sims summary.

The binary endpoint model is still calibrated through the event-time simulator. When the binary endpoint corresponds to event status by end_of_study, choose hazards that reproduce clinically meaningful endpoint probabilities through prop_to_haz() and ppwe().