--- title: "Bayesian binary outcome designs" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Bayesian binary outcome designs} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) set.seed(5107) ``` ```{r setup, message=FALSE} 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: ```{r priors} 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`. ```{r two_arm_design} 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 ``` 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: ```{r compare_binary_imputation} 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") ) ``` 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`. ```{r single_arm_design} 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 ``` 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. ```{r oc, eval=FALSE} 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()`.