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