--- 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 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. 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 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. ## 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_surv = hazard_prior, prior_bin = prior_bin, 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_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 ``` 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_surv = hazard_prior, prior_bin = prior_bin, 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()`.