Bayesian piecewise-exponential designs

library(goldilocks)

The “Two-arm randomized trials” vignette uses method = "logrank" and a single constant hazard for each treatment group. That is convenient when the proportional-hazards assumption is reasonable and the underlying event rate is well-approximated by an exponential distribution. In practice, neither assumption always holds. This vignette shows how to set up a Goldilocks design with:

a piecewise-exponential hazard, in which the constant hazard rate changes at one or more cut-points; and
a Bayesian decision rule (method = "bayes-surv") using a posterior probability threshold on the cumulative-failure-probability scale.

We first use a small example to show the analysis of one adaptive trial, then outline the repeated simulations needed to evaluate operating characteristics.

The Bayesian rule here targets the difference in event probabilities at end_of_study. If the prespecified endpoint is instead average event-free time through a fixed horizon, see the RMST vignette. RMST uses the same piecewise-exponential predictive machinery but a different completed-data test, effect scale, and direction of benefit (alternative = "greater").

When piecewise hazards help

Survival in many settings – post-surgical mortality, transplant outcomes, oncology with early treatment toxicity – shows a higher hazard early and a lower steady-state hazard later. A single exponential rate cannot capture both regimes; using the average rate over-states early survival and under-states late survival. A piecewise-exponential model with a single internal cut-point is often a good compromise, fitting one hazard for the early interval and another for everything afterwards.

The package handles piecewise hazards through related generation and analysis arguments:

generation_cutpoints defines the intervals used to generate event times. hazard_treatment and hazard_control supply one constant hazard per generating interval. For example, PWEALL represents generation_cutpoints = 6 with pieces [0, 6) and [6, \infty).
cutpoints defines the intervals used for posterior estimation, predictive imputation, and final piecewise-exponential analysis. When goldilocks assigns a realized event time to these intervals, it follows the survival counting-process convention (0,6] and (6,\infty): an event recorded exactly at month 6 belongs to the first interval.

By default generation_cutpoints = cutpoints, so one vector describes both models as in earlier package versions. The partitions may instead differ, and they need not contain the same number of intervals. The different boundary representations at an isolated cutpoint give the same continuous event-time distribution because the cutpoint itself has probability zero.

For design specification, prop_to_haz() maps cumulative event probabilities at selected times to the corresponding piecewise hazards.

Setting up the design

Suppose the primary endpoint is overall survival at 24 months. Based on prior data, we expect the control arm to have:

30% mortality by 6 months (an “early high-risk” window), and
50% mortality by 24 months overall.

The treatment is hypothesised to attenuate the early-period hazard, leaving the late-period hazard unchanged. Concretely, we target:

18% mortality by 6 months in the treatment arm, and
40% mortality by 24 months overall.

We set the cut-point at 6 months and translate these proportions into piecewise hazards:

cutpoints <- 6         # one change-point at 6 months gives two intervals
end_of_study <- 24

hc <- prop_to_haz(probs = c(0.30, 0.50), cutpoints = cutpoints, endtime = end_of_study)
ht <- prop_to_haz(probs = c(0.18, 0.40), cutpoints = cutpoints, endtime = end_of_study)

round(rbind(control = hc, treatment = ht), 4)
#>             [,1]   [,2]
#> control   0.0594 0.0187
#> treatment 0.0331 0.0174

The first column is the generating hazard during [0, 6) months and the second is the hazard from 6 months onward. In posterior sufficient statistics, the same columns receive observed events from (0,6] and (6,\infty), respectively. Both arms have a higher hazard in the early window. The implied 24-month survival probabilities can be verified with ppwe():

ppwe(hazard       = matrix(hc, nrow = 1),
     cutpoints    = cutpoints,
     end_of_study = end_of_study)
#> [1] 0.5
ppwe(hazard       = matrix(ht, nrow = 1),
     cutpoints    = cutpoints,
     end_of_study = end_of_study)
#> [1] 0.4

These should be 0.50 and 0.40 respectively (modulo rounding).

The Bayesian decision rule

With method = "bayes-surv", survival_adapt() puts independent \operatorname{Gamma}(\alpha, \beta) priors on each piecewise hazard rate (one per interval, per treatment group) and combines them with the observed exposure time and event counts to obtain a closed-form Gamma posterior on each \lambda_j. Posterior draws of \lambda_j are pushed through the piecewise-exponential cumulative incidence function to obtain posterior draws of the cumulative-failure probability p at end_of_study for the treatment and control groups.

The decision rule is one-sided. The treatment effect is defined as

\Delta = p_{\text{treatment}} - p_{\text{control}},

i.e. the difference in failure (not survival) probabilities at end_of_study. On the survival scale this is equivalent to S_{\text{control}} - S_{\text{treatment}}, with a negative \Delta corresponding to the treatment having higher survival. With alternative = "less" and a margin h0 (default 0), the trial declares success at the final analysis when

\Pr(\Delta < h_0 \mid \text{data}) \;>\; \texttt{prob\_ha}.

Because a beneficial treatment has a lower failure probability, alternative = "less" is the appropriate choice here. (method = "bayes-surv" does not allow alternative = "two.sided" – it raises an error.)

At each interim look, prior_surv is updated with the observed events and exposure to generate predictive completions for outstanding outcomes. Each completed dataset is then analyzed using prior_surv_final, the prior for the actual final analysis. The default prior_surv_final = prior_surv uses the same prior for both roles, as in this example. A separate informative predictive prior may incorporate external evidence while the final analysis prior remains weak. This requires explicitly supplying prior_surv_final; omitting it also uses the predictive prior in the analysis. The fraction of imputations that would declare success after enrollment continues to the maximum sample size is compared with Fn for the futility rule. Separately, the fraction that would declare success after completing follow-up for the subjects currently enrolled is compared with Sn for the expected-success rule. The default Qn = 1 disables the optional immediate-success rule in this example.

We use an independent weakly informative \operatorname{Gamma}(0.1, 0.1) prior on every hazard component:

prior_surv <- c(0.1, 0.1) # shape and rate for each lambda_j

In this example, leaving prior_surv_final at its default makes it equal to prior_surv. See the observed-interim example for an executable example with an informative predictive prior and a separate diffuse analysis prior. Both current-sample and maximum-sample predictions use that analysis prior to test their completed datasets.

What a simulated trial dataset looks like

Before running the adaptive design, it is helpful to inspect the subject-level data generated by sim_comp_data(). The following example uses the same hazards, follow-up horizon, and monthly time unit as the design below:

set.seed(7195)

example_trial_data <- sim_comp_data(
  hazard_treatment = ht,
  hazard_control = hc,
  generation_cutpoints = cutpoints,
  N_total = 12,
  lambda = 5,
  lambda_time = NULL,
  end_of_study = end_of_study,
  block = 4,
  rand_ratio = c(control = 1, treatment = 1),
  prop_loss = 0.05
)

knitr::kable(head(example_trial_data), digits = 2)

time	treatment	event	enrollment	id	loss_to_fu
24.00	1	0	0.00	1	FALSE
3.61	0	1	0.34	2	FALSE
1.90	0	1	0.35	3	FALSE
5.34	1	1	0.63	4	FALSE
19.77	0	1	0.72	5	FALSE
3.10	0	1	0.80	6	FALSE

Each row represents one simulated subject:

time is follow-up time from enrollment/randomization to the event or censoring, measured in months here.
treatment is the randomized arm: 1 for treatment and 0 for control.
event is 1 when the event was observed and 0 when follow-up was right-censored.
enrollment is trial-calendar time from first patient in to that subject’s enrollment, also measured in months here.
id is the subject’s simulated identifier.
loss_to_fu indicates whether the subject was censored because of simulated loss to follow-up.

The time and enrollment columns therefore use the same unit but different clocks: time is subject-relative follow-up, whereas enrollment is trial-calendar time.

Here prop_loss = 0.05 means \Pr(D\leq24)=0.05 for an independent exponential dropout time with rate -\log(0.95)/24 per month. The observed time is the minimum of the event time, dropout time, and 24 months. An event before dropout is retained, so the observed proportion censored by dropout can be below 5% and can differ between arms despite a common dropout distribution. The number of losses varies across trials; small examples may have none.

A single simulated trial

We simulate one trial under the alternative hypothesis to illustrate the statistical calculations. The interim analysis is scheduled after 60 participants have enrolled; the minimum allowable value is the randomization block size, max(block) = 4. A constant enrollment rate of five participants per month is assumed. With lambda_time = NULL, the first participant is placed at time zero and subsequent inter-arrival times follow an exponential distribution with rate 5. For piecewise accrual, lambda = c(2, 5) and lambda_time = 6 would specify two expected enrollments per month through month 6 and five per month thereafter. Enrollment-rate change-points use trial calendar time, whereas hazard cutpoints use participant follow-up time.

set.seed(7194)

out <- survival_adapt(
  hazard_treatment = ht,
  hazard_control   = hc,
  cutpoints        = cutpoints,
  generation_cutpoints = cutpoints,
  N_total          = 100,
  lambda           = 5,                # enrollments per month
  lambda_time      = NULL,             # constant enrollment rate
  interim_look     = 60,
  end_of_study     = end_of_study,
  prior_surv       = prior_surv,
  block            = 4,
  rand_ratio       = c(control = 1, treatment = 1),
  prop_loss        = 0.05,
  alternative      = "less",
  h0               = 0,
  Fn               = 0.05,
  Sn               = 0.95,
  prob_ha          = 0.975,
  N_impute         = 50,
  N_mcmc           = 2000,
  method           = "bayes-surv")

out
#>   prob_threshold margin alternative N_treatment N_control N_enrolled N_max
#> 1          0.975      0        less          50        50        100   100
#>   post_prob_ha  est_final ppp_success stop_futility stop_immediate_success
#> 1       0.9525 -0.1687922        0.12             0                      0
#>   stop_expected_success trial_success     stopping_reason decision_time
#> 1                     0         FALSE maximum_sample_size      43.91658
#>   accrual_stop_time analysis_ready_time planned_completion_time
#> 1           20.2688            43.91658                 44.2688
#>   followup_person_time peak_active_followup
#> 1             1560.222                   72

For this trial replicate, post_prob_ha reports the posterior probability of the alternative at the final (or stopped) analysis (post_prob_ha), the posterior mean treatment effect on the cumulative-failure scale (est_final), the predictive probability of success (ppp_success), and indicators for immediate success, expected-success stopping, or futility. Immediate success is disabled in this example by the default Qn = 1.

Notes on the piecewise model

Two practical considerations are worth flagging:

Empty intervals at interim looks. Early interim looks may have no subjects with follow-up reaching the later piecewise intervals. The empty_interval argument controls how these intervals are handled. The default, empty_interval = "prior", leaves such intervals at zero exposure and zero events, making their posteriors prior-driven. The legacy empty_interval = "propagate" option reproduces historical behavior by copying exposure time and event counts from the nearest non-empty interval within the same treatment group and emitting a warning. This is a sensitivity or migration option, not observed evidence about the empty interval. Use empty_interval = "error" to stop the analysis whenever an empty interval is encountered. By the final analysis, all intervals will typically be populated.
Number of cut-points. Each additional cut-point adds two hazard parameters to estimate in a two-arm design (one per treatment group). With limited interim data this can make individual interval posteriors diffuse. In our experience, one or two well-motivated cut-points (e.g., tied to a clinical milestone) is usually sufficient; finer partitions tend to add variance without commensurate bias reduction.

Sensitivity to the cut-point specification

If you suspect a piecewise structure but are unsure where the analysis cut-point should sit, fix hazard_treatment, hazard_control, and generation_cutpoints as the data-generating truth, then vary cutpoints across separate sim_trials() calls. The hazard vectors continue to match the fixed generation partition, while each analysis prior must match the candidate analysis partition. If the operating characteristics are similar, the design is robust to the analysis cut-point specification. If they differ markedly, the cut-point becomes a design decision worth justifying in the protocol.

A simpler – but very different – comparison is the equivalent design that is both simulated and analyzed under a single constant hazard for each treatment group matched to the overall 24-month proportions:

hc_flat <- prop_to_haz(0.50, endtime = end_of_study)   # control, single hazard
ht_flat <- prop_to_haz(0.40, endtime = end_of_study)   # treatment, single hazard

out_flat <- survival_adapt(
  hazard_treatment = ht_flat,
  hazard_control   = hc_flat,
  cutpoints        = NULL,
  N_total          = 100,
  lambda           = 5,
  lambda_time      = NULL,
  interim_look     = 60,
  end_of_study     = end_of_study,
  prior_surv       = prior_surv,
  block            = 4,
  rand_ratio       = c(control = 1, treatment = 1),
  prop_loss        = 0.05,
  alternative      = "less",
  h0               = 0,
  Fn               = 0.05,
  Sn               = 0.95,
  prob_ha          = 0.975,
  N_impute         = 50,
  N_mcmc           = 2000,
  method           = "bayes-surv")

Note that this changes both the simulated data-generating process and the analysis model, so any difference in operating characteristics conflates the two effects. It is most useful when the question is “how would the trial behave if the world really were a single exponential?” rather than “how robust is my analysis cut-point?”.