Skip to contents

Bayesian binary outcome designs

Source: vignettes/bayesian-binary.Rmd
bayesian-binary.Rmd

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 kinds of prior can therefore enter a Bayesian binary design. prior_surv is the Gamma prior on piecewise-exponential hazards used for interim predictive imputation, and prior_surv_final optionally overrides it for final imputation. prior_bin is the Beta prior on the binary endpoint event probability in each arm and is used to analyze every completed or imputed dataset. Thus, interim predictive probabilities are affected by prior_surv and prior_bin; a final imputed analysis uses prior_surv_final and prior_bin; and a non-imputed final binary analysis is affected by prior_bin alone.

The imputation and completed-data analysis are separate models. The piecewise-exponential Gamma model uses event times and available event-free follow-up to predict pending endpoint statuses. The beta-binomial model then analyzes the completed event indicators at end_of_study. prior_bin is not derived from the binary-horizon distribution implied by prior_surv, and goldilocks does not treat the two stages as one joint Bayesian model. This modular combination may be intentional, but its operating characteristics can depend on both prior specifications whenever outcomes require imputation. Choosing binary_imputation = "bernoulli" rather than "event-time" does not remove the distinction: both options use the same piecewise-exponential prediction model.

The saved analysis specification records prior_surv, prior_surv_final, prior_bin, and end_of_study, together with the Gamma prior parameters by analysis stage, arm, and interval. These values should be included in the simulation report because both prior families can affect interim decisions.

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
prior_bin <- 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.

The examples use prop_loss = 0. If dropout is added, prop_loss specifies the CDF of an independent exponential dropout time at end_of_study, not an exact fraction of missing binary outcomes. Events occurring before dropout remain observed. Consequently, excluding early-censored subjects can bias binary event probabilities even though censoring is independent for survival analysis. Assess final imputation and its event-time-model assumptions when evaluating binary designs with dropout.

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_surv = hazard_prior,
  prior_bin = prior_bin,
  bin_method = "quadrature",
  block = 2,
  rand_ratio = c(control = 1, treatment = 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_immediate_success
#> 1    0.9178978 -0.1129032        0.85             0                      0
#>   stop_expected_success trial_success     stopping_reason decision_time
#> 1                     0         FALSE maximum_sample_size        22.911
#>   accrual_stop_time analysis_ready_time planned_completion_time
#> 1            10.911              22.911                  22.911
#>   followup_person_time peak_active_followup
#> 1             1193.166                  106

For method = "bayes-bin", est_final is the posterior mean binary effect: the treatment event probability minus the control event probability. post_prob_ha 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 two approaches imply the same endpoint distribution. Finite Monte Carlo estimates can nevertheless differ, so comparisons should use adequate numbers of trials and imputations. The following small comparison uses a common seed to make the numerical results reproducible:

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 can also be more numerically stable for very small survival probabilities because it works from the remaining cumulative hazard rather than subtracting two nearly zero survival probabilities.

Whichever imputation approach is selected, the completed beta-binomial analysis depends only on the event count and denominator in each arm. With bin_method = "mc", each predictive replicate receives an independent set of N_mcmc beta-posterior draws. Thus, the number of predictive imputations and the number of posterior draws represent distinct sources of Monte Carlo error.

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_surv = hazard_prior,
  prior_bin = prior_bin,
  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_immediate_success
#> 1    0.9980375 0.1346154           1             0                      0
#>   stop_expected_success trial_success  stopping_reason decision_time
#> 1                     1          TRUE expected_success      17.96506
#>   accrual_stop_time analysis_ready_time planned_completion_time
#> 1           5.96506            17.96506                17.96506
#>   followup_person_time peak_active_followup
#> 1             564.1026                   49

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" approximates the posterior event probability or treatment-control difference by a normal distribution with matching mean and variance.
  • 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.

Approximation error can move the posterior probability across prob_ha, changing success and interim stopping decisions. Increasing N_mcmc does not correct normal-approximation error. Use bin_method = "quadrature" or sufficiently precise Monte Carlo Beta draws when the approximation is unsuitable.

Operating characteristics

As with the survival methods, operating characteristics should be evaluated by simulation under clinically relevant null and alternative scenarios. The seed argument makes the Monte Carlo study reproducible, including when ncores is greater than 1.

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_surv = hazard_prior,
  prior_bin = prior_bin,
  bin_method = "quadrature",
  block = 2,
  rand_ratio = c(control = 1, treatment = 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,
  "null: treatment event probability 35%" = out_null
))
effect_by_scenario <- c(
  "target: treatment event probability 25%" = 0.25,
  "null: treatment event probability 35%" = 0.35
)
oc$true_treatment_event_probability <- unname(effect_by_scenario[oc$scenario])

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 estimated operating characteristics.

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().