Package {goldilocks}

Title:	Adaptive Trial Designs for Survival and Binary Endpoints
Version:	0.6.0
Description:	Implements Goldilocks adaptive trial designs for time-to-event and fixed-time binary endpoints. Outcomes are generated with a piecewise exponential model, with conjugate Gamma priors used for predictive imputation. Final analyses may use log-rank or Cox tests, Bayesian piecewise-exponential inference, frequentist risk differences, or Bayesian beta-binomial inference. The method closely follows Broglio and colleagues <doi:10.1080/10543406.2014.888569> and supports simulation of design operating characteristics.
License:	GPL-3
Encoding:	UTF-8
URL:	https://graemeleehickey.github.io/goldilocks/, https://github.com/graemeleehickey/goldilocks, https://CRAN.R-project.org/package=goldilocks
BugReports:	https://github.com/graemeleehickey/goldilocks/issues
Depends:	R (≥ 4.1.0), survival
Imports:	dplyr, parallel, pbmcapply, PWEALL, Rcpp, rlang, stats
Suggests:	bench, covr, DiagrammeR, testthat (≥ 3.0.0), knitr, rmarkdown
LinkingTo:	BH, Rcpp
Language:	en-US
RoxygenNote:	7.3.3
Config/testthat/edition:	3
VignetteBuilder:	knitr
NeedsCompilation:	yes
Packaged:	2026-07-29 16:25:42 UTC; hickeg3
Author:	Graeme L. Hickey [aut, cre], Ying Wan [aut], Thevaa Chandereng [aut] (bayesDP code as a template), Becton, Dickinson and Company [cph], Tim Kacprowski [ctb] (For code from fastlogrank R package.)
Maintainer:	Graeme L. Hickey <graemeleehickey@gmail.com>
Repository:	CRAN
Date/Publication:	2026-07-29 17:00:11 UTC

goldilocks: Adaptive Trial Designs for Survival and Binary Endpoints

Description

Implements Goldilocks adaptive trial designs for time-to-event and fixed-time binary endpoints. Outcomes are generated with a piecewise exponential model, with conjugate Gamma priors used for predictive imputation. Final analyses may use log-rank or Cox tests, Bayesian piecewise-exponential inference, frequentist risk differences, or Bayesian beta-binomial inference. The method closely follows Broglio and colleagues doi:10.1080/10543406.2014.888569 and supports simulation of design operating characteristics.

Author(s)

Maintainer: Graeme L. Hickey graemeleehickey@gmail.com (ORCID)

Authors:

• Ying Wan ying.wan@bd.com

• Thevaa Chandereng tc3123@cumc.columbia.edu (ORCID) (bayesDP code as a template)

Other contributors:

• Becton, Dickinson and Company [copyright holder]

• Tim Kacprowski t.kacprowski@tu-braunschweig.de (For code from fastlogrank R package.) [contributor]

See Also

Useful links:

• https://graemeleehickey.github.io/goldilocks/

• https://github.com/graemeleehickey/goldilocks

• https://CRAN.R-project.org/package=goldilocks

• Report bugs at https://github.com/graemeleehickey/goldilocks/issues

Simulate exact continuous-time enrollment

Description

Simulate enrollment times from a Poisson process with a piecewise-constant rate.

Usage

enrollment(lambda = 1, N_total, lambda_time = NULL)

Arguments

Details

Major behavior change from goldilocks 0.5.0 and earlier. Versions through 0.5.0 generated Poisson counts in unit-time bins. enrollment() returned rebased integer bin times, after which sim_comp_data() added independent uniform jitter and sorted the result. The current implementation instead generates the continuous arrival times directly from the exact piecewise-constant Poisson process. Consequently:

• seeded simulations do not reproduce enrollment or downstream trial results obtained with version 0.5.0 or earlier;

• enrollment times and operating-characteristic estimates can change, particularly when rates are low or a rate change is not an integer time;

• the schedule API now contains internal knots only: change lambda_time = 0 to lambda_time = NULL for a constant rate, and change lambda_time = c(0, t1, t2) to lambda_time = c(t1, t2) for a piecewise rate.

enrollment() treats time zero as first patient in, matching the time origin used throughout goldilocks. The first returned enrollment time is therefore fixed at zero. The remaining N_total - 1 arrivals form a continuous-time non-homogeneous Poisson process whose rate is constant between the supplied internal knots. Thus, lambda is measured in enrollments per unit of lambda_time; for example, when time is measured in months, lambda = 5 means five enrollments per month on average.

Write the internal knots as 0 < \tau_1 < \cdots < \tau_K, with \tau_0 = 0 implicit. The enrollment intensity is

\lambda(t) = \lambda_j, \qquad \tau_{j-1} \le t < \tau_j,

for j = 1, \ldots, K, and \lambda(t) = \lambda_{K+1} after the final knot. The final rate therefore continues for as long as needed to reach N_total; there is no finite accrual horizon in this function.

Arrivals are generated exactly by the time-rescaling theorem. If E_2, \ldots, E_N are independent unit-rate exponential variables and S_i = \sum_{k=2}^i E_k, then the enrollment times after first patient in are

T_1 = 0, \qquad T_i = \Lambda^{-1}(S_i),

where \Lambda(t) = \int_0^t \lambda(u)\,du is the cumulative enrollment intensity. This construction gives independent Poisson counts on disjoint calendar intervals, with expected count \int_a^b \lambda(u)\,du over (a,b]. For a constant rate, successive enrollment gaps are independent Exponential(lambda) variables.

The rate-change times are measured from first patient in, not from an earlier site-opening or trial-activation date. If operational delays before first patient in are important, they must be modelled separately before using the returned relative times.

For example, lambda = c(0.3, 0.7, 0.9, 1.2) with lambda_time = c(5, 10, 15) specifies average enrollment rates of 0.3 over [0,5), 0.7 over [5,10), 0.9 over [10,15), and 1.2 from time 15 onward. Fractional knots such as lambda_time = 2.5 are handled exactly; no unit-time binning or post-hoc jitter is used.

Value

A non-decreasing numeric vector of N_total continuous enrollment times, measured from first patient in and expressed in the same time unit used for lambda_time. The first value is zero.

Examples

# Constant enrollment: first patient at zero, then exponential gaps.
enrollment(lambda = 0.7, N_total = 10)

# Three internal rate changes define four enrollment intervals.
enrollment(
  lambda = c(0.3, 0.7, 0.9, 1.2),
  N_total = 50,
  lambda_time = c(5, 10, 15)
)

# Fractional change times are supported exactly.
enrollment(
  lambda = c(0.25, 1),
  N_total = 20,
  lambda_time = 2.5
)

goldilocks

Description

The goal of goldilocks is to implement the Goldilocks Bayesian adaptive design proposed by Broglio et al. (2014), for one- and two-arm trials. Outcomes are generated with an underlying piecewise-exponential event-time model. Final analyses may retain the time-to-event outcome (method = "logrank", "cox", or "bayes-surv") or reduce complete follow-up to event status at a fixed endpoint time (method = "riskdiff" or "bayes-bin").

The method can be used for a confirmatory trial to select a sample size based on accumulating data. During accrual, predictive probabilities are used to determine whether the current sample size is sufficient, whether continuing accrual would be futile, or whether enrollment should continue. The algorithm explicitly accounts for completion of planned follow-up before the primary analysis. Broglio et al. (2014) refer to this as a Goldilocks trial design, as it is constantly asking the question, “Is the sample size too big, too small, or just right?”

References

Broglio KR, Connor JT, Berry SM. Not too big, not too small: a Goldilocks approach to sample size selection. Journal of Biopharmaceutical Statistics, 2014; 24(3): 685–705.

Plot an enrollment projection

Description

Draws the expected cumulative enrollment curve for a Goldilocks trial design, together with optional random enrollment trajectories and projected interim and maximum-sample-size milestones.

Usage

plot_enrollment(
  x = NULL,
  lambda = NULL,
  N_total = NULL,
  lambda_time = NULL,
  interim_look = NULL,
  end_of_study = NULL,
  n_sim = 20L,
  seed = NULL,
  time_unit = NULL,
  xlab = NULL,
  ylab = "Cumulative number of enrolled patients",
  main = NULL,
  annotate = TRUE,
  projection_col = "#276E9B",
  simulation_col = "#777777",
  milestone_col = "#C8682A"
)

Arguments

Details

The blue projection is 1 + \Lambda(t), where \Lambda(t) is the cumulative intensity of the piecewise-constant Poisson enrollment process. The first patient is fixed at time zero, consistently with enrollment(). A milestone's projected time solves 1 + \Lambda(t) = N. With a constant enrollment rate this is also the mean arrival time, (N - 1) / \lambda. With a piecewise rate it is an expected-count projection rather than the mean of the corresponding arrival-time distribution.

If x supplies a stored design, explicitly supplied design arguments override the corresponding stored values. This makes it possible, for example, to compare a fitted design with a different enrollment rate.

Value

Invisibly, a list containing the evaluated design, the projection data frame, the milestones data frame, and the simulated enrollment-time vectors in simulations.

Examples

plot_enrollment(
  lambda = 20,
  N_total = 600,
  interim_look = 400,
  end_of_study = 12,
  n_sim = 20,
  seed = 20260727,
  time_unit = "months"
)

# Piecewise enrollment rates are supported.
plot_enrollment(
  lambda = c(8, 20),
  lambda_time = 6,
  N_total = 200,
  interim_look = c(100, 150),
  n_sim = 5,
  seed = 1,
  time_unit = "months"
)

Plot predictive-probability decision maps

Description

Draws one decision map per interim look across simulated trials. The horizontal axis is the predictive probability of success if enrollment continues to the maximum sample size, and the vertical axis is the predictive probability of success if enrollment stops at the current sample size. Shaded regions and dashed lines show the futility, continuation, and expected-success rules. Identical points are aggregated; point size indicates their frequency.

Usage

plot_sim_decisions(x)

Arguments

Value

The simulation traces, invisibly.

Examples

traces <- data.frame(
  trial = 1:6,
  look = rep(c(1, 2), each = 3),
  planned_N = rep(c(40, 60), each = 3),
  ppp_stop_now = c(0.96, 0.4, 0.2, 0.92, 0.5, 0.1),
  success_threshold = rep(c(0.95, 0.9), each = 3),
  ppp_success_at_max = c(0.8, 0.5, 0.02, 0.85, 0.4, 0.01),
  futility_threshold = 0.05,
  decision = c(
    "stop_expected_success", "continue", "stop_futility",
    "stop_expected_success", "continue", "stop_futility"
  )
)
plot_sim_decisions(traces)

Plot operating characteristics across simulation scenarios

Description

Draws operating-characteristic curves across a series of true treatment-effect scenarios. The first panel shows final success and stopping probabilities. The second panel shows mean enrolled sample size.

Usage

plot_sim_ocs(x, effect, xlab = "True treatment effect")

Arguments

Value

x, invisibly.

Examples

operating_characteristics <- data.frame(
  scenario = c("null", "small", "target"),
  effect = c(1, 0.85, 0.7),
  power = c(0.025, 0.55, 0.9),
  stop_success = c(0.01, 0.35, 0.75),
  stop_futility = c(0.7, 0.25, 0.05),
  stop_max_N = c(0.29, 0.4, 0.2),
  mean_N = c(180, 250, 210),
  sd_N = c(45, 60, 55),
  stop_and_fail = c(0.001, 0.01, 0.02)
)
plot_sim_ocs(
  operating_characteristics,
  effect = "effect",
  xlab = "True hazard ratio"
)

Plot stopping outcomes from trial simulations

Description

Draws a stacked bar chart of stopping outcomes by enrolled sample size, with colours distinguishing expected-success, futility, and maximum-sample-size outcomes. The type argument controls whether the function draws marginal, conditional, or cumulative bars, or a flowchart through successive interim looks. Bar-chart subtitles state the denominator used by the selected view. The input can be the sims element returned by sim_trials() or the complete sim_trials() result.

Usage

plot_sim_stopping(
  x,
  type = c("marginal", "conditional", "cumulative", "flowchart")
)

Arguments

Details

The marginal view uses terminal sample sizes observed in N_enrolled. When the complete result from sim_trials(return_trace = TRUE) is supplied, the conditional, cumulative, and flowchart views also include sample sizes recorded in traces, so reached looks at which no trial stopped still appear. The flowchart requires the N_max column and is rendered with DiagrammeR::grViz().

Value

For bar-chart types, the simulation result data frame, invisibly. For type = "flowchart", a DiagrammeR grViz htmlwidget.

Plot predictive probabilities and enrollment at interim looks

Description

Draws three base-R panels showing the predictive probability of success if accrual stops now, the predictive probability of success at the maximum sample size, and enrollment and observed events by treatment arm. Thresholds and early stopping decisions are marked on the probability panels.

Usage

plot_trial_trace(x)

Arguments

Value

The trace data frame, invisibly.

Cumulative distribution function of the PWE for a vectorized hazard rate parameter

Description

Extends PWEALL::pwe() to allow for vectorization over the hazard rates.

Usage

ppwe(hazard, end_of_study, cutpoints = NULL)

Arguments

Value

A vector of (0, 1) probabilities from evaluation of the PWE cumulative distribution function. Length of the vector matches the number of rows of the hazard matrix parameter.

Print an adaptive trial trace result

Description

Prints the final trial summary and reports how many interim looks were completed when survival_adapt returns a goldilocks_trial object.

Usage

## S3 method for class 'goldilocks_trial'
print(x, ...)

Arguments

Value

The input object, invisibly.

Estimate plausible piecewise constant hazard rates from summary summary event proportions

Description

Given estimates of the event probability at one or more fixed times, the corresponding piecewise hazard rates can be determined through closed-form formulae. This utility function can be useful when simulating trial datasets with plausible event rates.

Usage

prop_to_haz(probs, cutpoints = NULL, endtime)

Arguments

Details

Given J-1 interior cutpoints, then there are J intervals defined as: [s_0, s_1), [s_1, s_2), \dots, [s_{J-1}, s_{J}), with conditions that s_0 = 0 and s_J = \infty. Each interval corresponds to constant hazard \lambda_j.

Value

Vector of constant hazard rates for each time piece defined by cutpoints.

Examples

lambda <- prop_to_haz(0.15, endtime = 36) # 15% probability at 36-months
all.equal(pexp(36, lambda), 0.15)

# 15% probability at 12-months, and 30% at 24-months
prop_to_haz(c(0.15, 0.30), 12, 24)
PWEALL::pwe(12, prop_to_haz(c(0.15, 0.30), 12, 24), c(0, 12))$dist
PWEALL::pwe(24, prop_to_haz(c(0.15, 0.30), 12, 24), c(0, 12))$dist

Impute piecewise exponential time-to-event outcomes

Description

Imputation of time-to-event outcomes using the piecewise constant hazard exponential function conditional on observed exposure.

Usage

pwe_impute(time, hazard, cutpoints = NULL, maxtime = NULL)

Arguments

Details

If a subject is event-free at time s < t, then the conditional probability is

F_{T | s}(t | s) = P(T \le t | T > s) = \frac{F(t) - F(s)}{1 - F(s)}

where F(\cdot) is the cumulative distribution function of the piecewise exponential (PWE) distribution. Equivalently, F(t) = 1 - S(t), where S(t) is the survival function. If U \sim Unif(0, 1), then we can generate an event time (conditional on being event free up until s) as

F^{-1}(U(1 - F(s)) + F(s))

If s = 0, this is equivalent to a direct unconditional sample from the PWE distribution.

Value

A data frame with simulated follow-up times (time) and respective event indicator (event, 1 = event occurred, 0 = censoring).

Examples

pwe_impute(time = c(3, 4, 5), hazard = c(0.002, 0.01), cutpoints = 12)
pwe_impute(time = c(3, 4, 5), hazard = c(0.002, 0.01), cutpoints = 12,
           maxtime = 36)
pwe_impute(time = 19.621870008, hazard = c(2.585924e-02, 3.685254e-09),
           cutpoints = 12, maxtime = 36)

Simulate piecewise exponential time-to-event outcomes

Description

Simulate time-to-event outcomes using the piecewise constant hazard exponential function.

Usage

pwe_sim(n = 1, hazard = 1, cutpoints = NULL, maxtime = NULL)

Arguments

Details

See pwe_impute() for details.

Value

A data frame with simulated follow-up times (time) and respective event indicator (event, 1 = event occurred, 0 = censoring).

Examples

pwe_sim(10, hazard = c(0.005, 0.001), cutpoints = 3, maxtime = 36)
y <- pwe_sim(n = 1, hazard = c(2.585924e-02, 3.685254e-09),
             cutpoints = 12)

Randomization allocation

Description

Generates a randomized treatment assignment sequence for control and treatment arms with different randomization ratios and block sizes.

Usage

randomization(N_total, block = 2, allocation = c(1, 1))

Arguments

Details

Complete randomization may not always be ideal due to the chance of drawing a large block assigned to one treatment group, potentially impacting the time to enrollment completion. Therefore, a block randomization allocation may be preferable. The block randomization allocation specification allows for different randomization ratios, but they must be given in integer form. Additionally, the block size should be an integer that is divisible by the sum of the randomization allocation; see the examples.

Value

An integer treatment assignment vector, coded 0 for control and 1 for treatment.

Examples

# Implementing treatment allocation for control to treatment with 1:1.5
# randomization ratio
randomization(N_total = 100, block = 5, allocation = c(2, 3))

# Treatment allocation with 2:1 for control to treatment
randomization(N_total = 70, block = 9, allocation = c(2, 1))

# Treatment allocation for control to treatment with 1:2 for control
# to treatment with multiple block sizes c(3, 9, 6)
randomization(N_total = 100, block = c(3, 9, 6), allocation = c(1, 2))

# For complete randomization set the N_total to block size
randomization(N_total = 100, block = 100, allocation = c(1, 1))

Simulate a complete clinical trial with event data drawn from a piecewise exponential distribution

Description

Simulate a complete clinical trial with event data drawn from a piecewise exponential distribution

Usage

sim_comp_data(
  hazard_treatment,
  hazard_control = NULL,
  cutpoints = NULL,
  N_total,
  lambda = 0.3,
  lambda_time = NULL,
  end_of_study,
  block = 2,
  rand_ratio = c(1, 1),
  prop_loss = 0
)

Arguments

Details

Enrollment is simulated directly in continuous time by enrollment(). The first patient is placed at time zero and all subsequent enrollment times are measured from first patient in. No uniform jitter is added in sim_comp_data().

lambda_time and cutpoints both contain internal change times, but they describe different clocks. lambda_time describes changes in the trial's calendar-time enrollment rate measured from first patient in. cutpoints describes changes in an individual subject's event hazard measured from that subject's enrollment. They need not have the same values or length. All time quantities supplied to a simulation should nevertheless use one common unit, such as days or months.

Value

A data frame with 1 row per subject and columns:

• time: Time of event or censoring time.

• treatment: Treatment assignment, coded 1L for the treatment arm and 0L for the control arm. Single-arm designs have treatment = 1L for every subject.

• event: Indicator of whether event occurred (1L if occurred and 0L if right-censored).

• enrollment: Time of patient enrollment relative to the time the trial enrolled the first patient. The package treats enrollment and randomization as occurring at the same time.

• id: Identification number for each patient.

• loss_to_fu: Indicator of whether the patient was lost to follow-up during observation.

Simulate one or more clinical trials subject to known design parameters and treatment effect

Description

Simulate multiple clinical trials with fixed input parameters, and tidily extract the relevant data to generate operating characteristics.

Usage

sim_trials(
  hazard_treatment,
  hazard_control = NULL,
  cutpoints = NULL,
  N_total,
  lambda = 0.3,
  lambda_time = NULL,
  interim_look = NULL,
  end_of_study,
  prior_surv = c(0.1, 0.1),
  prior_bin = c(1, 1),
  bin_method = "mc",
  block = 2,
  rand_ratio = c(1, 1),
  prop_loss = 0,
  alternative = "greater",
  h0 = 0,
  Fn = 0.05,
  Sn = 0.9,
  prob_ha = 0.95,
  N_impute = 10,
  N_mcmc = 10,
  N_trials = 10,
  method = "logrank",
  imputed_final = FALSE,
  empty_interval = c("propagate", "prior", "error"),
  return_trace = FALSE,
  ncores = 1L,
  backend = c("auto", "fork", "psock", "sequential"),
  seed = NULL,
  binary_imputation = c("event-time", "bernoulli"),
  prior_surv_final = prior_surv
)

Arguments

Details

This is basically a wrapper function for survival_adapt(), whereby we repeatedly run the function for independent trials (all with the same input design parameters and treatment effect).

To use multiple cores (where available), the argument ncores can be increased from the default of 1. The default backend = "auto" uses pbmcapply::pbmclapply() on Unix-like platforms and a PSOCK cluster on Windows, where forked processes are unavailable. Set backend explicitly to compare backends or to require serial execution.

Set seed to make sim_trials() reproducible. When a seed is supplied, sim_trials() first generates one independent "L'Ecuyer-CMRG" stream for each simulated trial, then each call to survival_adapt() runs with its own per-trial stream. This avoids reusing the same random-number stream across workers when ncores > 1, and produces identical seeded results across supported backends. A seeded call restores the caller's RNG state on exit. With seed = NULL, the function uses and advances R's current global RNG state.

Value

A list containing sims, a data frame with one row per simulated trial, and call. When return_trace = TRUE, the list also contains traces, a data frame with one row per completed interim look and a trial identifier. See survival_adapt() for details of the summary and trace columns.

Examples

hc <- prop_to_haz(c(0.20, 0.30), 12, 36)
ht <- prop_to_haz(c(0.05, 0.15), 12, 36)

out <- sim_trials(
  hazard_treatment = ht,
  hazard_control = hc,
  cutpoints = 12,
  N_total = 600,
  lambda = 20,
  lambda_time = NULL,
  interim_look = c(400, 500),
  end_of_study = 36,
  prior_surv = c(0.1, 0.1),
  block = 2,
  rand_ratio = c(1, 1),
  prop_loss = 0.30,
  alternative = "two.sided",
  h0 = 0,
  Fn = 0.05,
  Sn = 0.9,
  prob_ha = 0.975,
  N_impute = 5,
  N_mcmc = 5,
  method = "logrank",
  N_trials = 2,
  ncores = 1,
  backend = "auto",
  seed = 123)

Summarize simulations to get operating characteristics

Description

Summarize simulations to get operating characteristics

Usage

summarise_sims(data)

Arguments

Value

Data frame reporting the operating characteristics, including the power (which will be equal to the type I error in the null case); the proportion of trials that stopped for early expected success, futility, or went to the maximum sample size. The average stopping sample size (and standard deviation) are also recorded. The proportion of trials that stopped early for expected success, yet went to ultimately fail are also reported.

Summarize an interim decision path

Description

Creates a compact one-row summary of the final interim look and stopping decision. Pass the goldilocks_trial object to include final analysis information, or pass its trace element to summarize the path only.

Usage

summarise_trial_trace(x)

Arguments

Value

A one-row data frame.

Simulate and analyze one Goldilocks adaptive trial

Description

Simulate and analyze one Goldilocks adaptive trial

Usage

survival_adapt(
  hazard_treatment,
  hazard_control = NULL,
  cutpoints = NULL,
  N_total,
  lambda = 0.3,
  lambda_time = NULL,
  interim_look = NULL,
  end_of_study,
  prior_surv = c(0.1, 0.1),
  prior_bin = c(1, 1),
  bin_method = "mc",
  block = 2,
  rand_ratio = c(1, 1),
  prop_loss = 0,
  alternative = "greater",
  h0 = 0,
  Fn = 0.05,
  Sn = 0.9,
  prob_ha = 0.95,
  N_impute = 10,
  N_mcmc = 10,
  empty_interval = c("propagate", "prior", "error"),
  method = "logrank",
  imputed_final = FALSE,
  return_trace = FALSE,
  binary_imputation = c("event-time", "bernoulli"),
  prior_surv_final = prior_surv
)

Arguments

Details

Implements the Goldilocks design method described in Broglio et al. (2014). At each interim analysis, two probabilities are computed:

1. The posterior predictive probability of eventual success. This is calculated as the proportion of imputed datasets at the current sample size that would go on to be success at the specified threshold. At each interim analysis it is compared to the corresponding element of Sn, and if it exceeds the threshold, accrual/enrollment is suspended and the outstanding follow-up allowed to complete before conducting the pre-specified final analysis.

2. The posterior predictive probability of final success. This is calculated as the proportion of imputed datasets at the maximum threshold that would go on to be successful. Similar to above, it is compared to the corresponding element of Fn, and if it is less than the threshold, accrual/enrollment is suspended and the trial terminated. Typically this would be a binding decision. If it is not a binding decision, then one should also explore the simulations with Fn = 0.

Hence, at each interim analysis look, 3 decisions are allowed:

1. Stop for expected success

2. Stop for futility

3. Continue to enroll new subjects, or if at maximum sample size, proceed to final analysis.

At each interim (and final) analysis methods as:

• Log-rank test (method = "logrank"). Each (imputed) dataset with both treatment and control arms can be compared using a standard log-rank test. The output is a P-value, and there is no treatment effect reported. The function returns 1 - P, which is reported in post_prob_ha. Whilst not a posterior probability, it can be contrasted in the same manner. For example, if the success threshold is P < 0.05, then one requires post_prob_ha > 0.95. The reason for this is to enable simple switching between Bayesian and frequentist paradigms for analysis. When alternative = "less" or "greater", a one-sided P-value is computed from the log-rank z-statistic.

• Cox proportional hazards regression Wald test (method = "cox"). Similar to the log-rank test, a P-value is calculated and 1 - P is reported in post_prob_ha. When alternative = "two.sided", the standard two-sided Wald P-value is used when h0 = 0. For other values of h0, the Wald test is centered on the specified null log hazard ratio. When alternative = "less" or "greater", a one-sided P-value is derived from the Wald z-statistic relative to h0. The treatment effect (log hazard ratio) is also reported. When imputed_final = TRUE, the Cox model is fitted separately to each of at least two imputed datasets. The log hazard ratios and their within-imputation variances are combined using Rubin's rules; the pooled Wald test uses Rubin's large-sample degrees of freedom. When imputed_final = FALSE, the existing single Cox model is fitted directly to the observed right-censored data.

• Bayesian absolute difference (method = "bayes-surv"). Each imputed dataset is used to update the conjugate Gamma prior (defined by prior_surv at interim looks and prior_surv_final at the final stage), yielding a posterior distribution for the piecewise exponential rate parameters. In turn, the posterior distribution of the cumulative incidence function (1 - S(t), where S(t) is the survival function) evaluated at time end_of_study is calculated. If a single-arm study, then this summarizes the treatment effect, else, if a two-armed study, the independent posteriors are used to estimate the posterior distribution of the difference. A posterior probability is calculated according to the specification of the test type (alternative) and the value of the null hypothesis (h0).

  For piecewise-exponential analyses, an interim or final dataset may contain intervals with no exposed subjects in one treatment arm, especially when later cutpoints occur after the available follow-up at early looks. The empty_interval argument controls this case. The default, "propagate", preserves historical package behavior by borrowing sufficient statistics from the nearest non-empty interval within the same treatment arm. This is operationally stable but statistically consequential because the empty interval's posterior is informed by adjacent observed data. "prior" instead leaves the empty interval prior-driven, making the absence of interval data explicit. "error" is strict and stops the simulation or analysis when an empty interval is encountered.

• Bayesian beta-binomial analysis (method = "bayes-bin"). Each complete or imputed dataset is reduced to binary event outcomes at end_of_study. A conjugate Beta(a, b) prior, specified with prior_bin, is updated with the number of events and non-events in each arm. In a single-arm study, inference is based on the posterior event probability. In a two-arm study, inference is based on p_\textrm{treatment} - p_\textrm{control}. This posterior probability can be calculated using Monte Carlo beta draws (bin_method = "mc"), a normal approximation ("normal"), or numerical quadrature ("quadrature"). Like the risk-difference test, this method requires complete binary outcomes: censored subjects must either be followed to end_of_study, imputed, or excluded when imputed_final = FALSE.

  Two equivalent predictive imputation approaches are available through binary_imputation. With "event-time", the package samples a future event time conditional on the available event-free follow-up and then records whether it falls by end_of_study. With "bernoulli", it calculates the same endpoint probability directly. If T is the observed event-free follow-up, T^* is end_of_study, S(t) is the survival function, and H(t) is the cumulative hazard, that probability is

  \Pr(X = 1 \mid T_\mathrm{event} > T) = \frac{S(T) - S(T^*)}{S(T)} = 1 - \exp\{-[H(T^*) - H(T)]\}.

  A Bernoulli outcome is drawn with this probability. For a subject not yet enrolled, T = 0; observed events are retained unchanged. Because no precise event time is generated, the imputed time is set to end_of_study and only the binary event status is analyzed. Each imputation still uses a sampled posterior hazard draw, so uncertainty in the piecewise-exponential model is retained.

• Frequentist risk difference (method = "riskdiff"). Each complete or imputed dataset is reduced to binary event outcomes at end_of_study. The estimated treatment effect is p_\textrm{treatment} - p_\textrm{control}, with an unpooled binomial variance. A Wald test compares this estimate with h0, and 1 - P is reported in post_prob_ha. All three alternatives are supported. Because the test requires complete binary outcomes, lost-to-follow-up subjects are excluded when imputed_final = FALSE. When imputed_final = TRUE, estimates and within-imputation variances from at least two completed datasets are combined using Rubin's rules.

• Imputed final analysis (imputed_final). The overall final analysis conducted after accrual is suspended and follow-up is complete can be analyzed on imputed datasets for Bayesian methods ("bayes-surv" and "bayes-bin"), Cox regression, and the frequentist risk-difference analysis, or on the non-imputed dataset. Since the imputations/predictions used during the interim analyses assume all subjects are imputed (since loss to follow-up is not yet known), it would seem most appropriate to conduct the trial in the same manner, especially if loss to follow-up rates are appreciable. Note, this only applies to subjects who are right-censored due to loss to follow-up, which we assume is a non-informative process. For Cox regression the final estimates and variances are pooled with Rubin's rules. It cannot be used with method = "logrank".

When method = "bayes-surv" or method = "bayes-bin" and imputation is involved (either at interim analyses or via imputed_final = TRUE), a two-stage posterior procedure is used. First, the posterior distribution of the piecewise hazard rates is estimated from the observed data and used to draw imputed event times for censored subjects. Second, a new posterior is estimated from the combined observed and imputed data: the piecewise-exponential posterior for method = "bayes-surv" or the beta posterior for method = "bayes-bin". This posterior is used for inference. This is consistent with the predictive probability framework described in Broglio et al. (2014), but users should be aware that the imputation model's posterior influences the analysis posterior. For frequentist methods ("logrank", "cox", "riskdiff"), each completed dataset uses a standard test rather than a posterior, so this feedback loop does not arise. Imputed Cox final analyses then pool the completed-data estimates and variances using Rubin's rules. Imputed risk-difference final analyses use the same scalar combining rule.

At each interim look, follow-up times are masked (censored) to reflect the calendar time of the analysis. The package treats enrollment and randomization as occurring at the same time. Subjects enrolled at the exact interim boundary have zero follow-up time. These times are clamped to .Machine$double.eps (approximately 2.2 \times 10^{-16}) so that they contribute negligible but non-zero exposure to the interim posterior. This affects at most one subject per interim look.

Value

With return_trace = FALSE (the default), a data frame containing some input parameters (arguments) as well as statistics from the analysis, including:

• N_treatment: Number of patients enrolled in the treatment arm.

• N_control: Number of patients enrolled in the control arm.

• est_final: Treatment effect estimated at the final analysis. The final analysis occurs when either the maximum sample size is reached and follow-up is complete, or the interim analysis triggered early stopping of enrollment/accrual and follow-up for those subjects is complete.

• post_prob_ha: Posterior probability from the final analysis. If a Bayesian method uses imputed_final = TRUE, this is calculated for each imputed final-analysis dataset and averaged over N_impute imputations. For an imputed Cox analysis it is 1 - P from the Rubin-pooled Wald test. The same interpretation applies to imputed risk-difference analyses. For non-imputed frequentist analyses it is 1 - P from the corresponding test.

• stop_futility: Logical indicator of whether the trial stopped early for futility.

• stop_expected_success: Logical indicator of whether the trial stopped early for expected success.

With return_trace = TRUE, a goldilocks_trial object is returned. Its summary element is the same data frame and its trace element has one row per interim look. The trace records enrollment and observed events by arm, predictive probabilities and their thresholds, the decision taken, and warnings raised during that look. It deliberately excludes imputed data sets and posterior draws to keep the output compact.

References

Broglio KR, Connor JT, Berry SM. Not too big, not too small: a Goldilocks approach to sample size selection. Journal of Biopharmaceutical Statistics, 2014; 24(3): 685–705.

Examples

# RCT with exponential hazard (no piecewise breaks)
# Note: the number of imputations is small to enable this example to run
#       quickly on CRAN tests. In practice, much larger values are needed.
survival_adapt(
 hazard_treatment = -log(0.85) / 36,
 hazard_control = -log(0.7) / 36,
 cutpoints = NULL,
 N_total = 600,
 lambda = 20,
 lambda_time = NULL,
 interim_look = 400,
 end_of_study = 36,
 prior_surv = c(0.1, 0.1),
 block = 2,
 rand_ratio = c(1, 1),
 prop_loss = 0.30,
 alternative = "less",
 h0 = 0,
 Fn = 0.05,
 Sn = 0.9,
 prob_ha = 0.975,
 N_impute = 10,
 N_mcmc = 10,
 method = "bayes-surv")