Package {goldilocks}

Title:	Adaptive Trial Designs for Survival and Binary Endpoints
Version:	1.0.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, Cox, or restricted mean survival time tests, Bayesian piecewise-exponential inference, frequentist risk differences, or Bayesian beta-binomial inference. The method closely follows Broglio and colleagues (2014) <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, survRM2, 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-09-09 11:38:15 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-09-09 14:40:02 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, Cox, or restricted mean survival time tests, Bayesian piecewise-exponential inference, frequentist risk differences, or Bayesian beta-binomial inference. The method closely follows Broglio et al. (2014) 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

Simulates 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. From version 0.6.0, enrollment times are drawn directly from the exact continuous-time 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 lambda_time argument now contains internal change times 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. The value of the intensity at an isolated knot does not change the Poisson process. When assigning a realized enrollment time to an interval, goldilocks follows the survival counting-process convention: intervals are open on the left and closed on the right, so an arrival exactly at \tau_j belongs to the interval ending at \tau_j. The first patient at time zero is a fixed origin and is handled separately.

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 positive times in (0,5], 0.7 over (5,10], 0.9 over (10,15], and 1.2 after time 15. 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
)

Evaluate an externally observed interim data cut

Description

Applies the same posterior predictive decision calculation used by survival_adapt() to subject-level data observed at one interim look. The function does not simulate a trial or modify data. Participants who could still be enrolled before the maximum sample size are represented according to N_total, the observed arm counts, and rand_ratio.

Usage

evaluate_interim(
  data,
  data_cut,
  look,
  N_total,
  end_of_study,
  cutpoints = NULL,
  prior_surv = c(0.1, 0.1),
  prior_bin = c(1, 1),
  bin_method = "mc",
  rand_ratio = c(control = 1, treatment = 1),
  single_arm = FALSE,
  alternative = "greater",
  h0 = 0,
  Fn = 0.05,
  Sn = 0.9,
  prob_ha = 0.95,
  N_impute = 500,
  N_mcmc = 1000,
  mc_conf_level = 0.95,
  empty_interval = c("prior", "propagate", "error"),
  method = "logrank",
  binary_imputation = c("event-time", "bernoulli"),
  seed = NULL,
  Qn = 1,
  rmst_tau = end_of_study,
  prior_surv_final = prior_surv
)

Arguments

Details

The arm-specific maximum enrollment is N_total * rand_ratio / sum(rand_ratio). Potential future accrual in each arm is its maximum enrollment minus its observed enrollment. Randomization block information and concealed future assignment order are unnecessary: the maximum-sample predictive calculation requires only the number of potential future participants in each arm.

The data$status values distinguish observed follow-up: "event" for an observed endpoint event, "complete" for event-free completion of end_of_study, "pending" for a subject still under follow-up, and "censored" for permanent early censoring. Pending and censored outcomes are predictively imputed conditional on time.

With method = "bayes-surv", prior_surv generates predictive outcomes and prior_surv_final is the analysis prior used to test each hypothetical completed trial for success, at both the current and maximum sample sizes. Supply the same final prior as in the prespecified trial design. The default prior_surv_final = prior_surv uses one prior for both roles. Other methods do not use prior_surv_final in this interim calculation; Bayesian binary completed-data analyses use prior_bin.

Qn, Sn, and Fn are scalar thresholds for this look. Immediate success is declared when the estimated probability of completed-data success among the current participants is strictly greater than Qn. Otherwise, expected success stops accrual for planned follow-up when that probability is strictly greater than Sn. If neither success rule applies, futility is declared when the corresponding maximum-sample probability is strictly less than Fn. Set Fn = 0 or NULL to disable futility. Exact one-sided Monte Carlo bounds are returned as diagnostics and do not drive any decision.

The function requires treatment assignments to perform the arm-specific posterior and completed-data analyses. In a blinded trial, an independent unblinded statistician or service should join the treatment assignments, run this function, and return the aggregate decision without distributing the subject-level input.

Value

An object of class goldilocks_interim, containing:

• decision: a one-row decision summary;

• probabilities: current- and maximum-sample predictive probabilities;

• monte_carlo: estimates, standard errors, and diagnostic bounds;

• diagnostics: observed status counts, potential accruals, warnings, imputation diagnostics, and resolved Gamma prior and posterior parameters by arm and interval;

• trace: a one-row decision trace compatible with plot_trial_trace() and summarise_trial_trace();

• metadata: the evaluated design, resolved prior design, package version, time-origin, data-cut, and random-number policy. For method = "bayes-surv", both Gamma priors are retained in metadata$design and metadata$prior_design; posterior diagnostics describe the predictive model based on the observed interim data. For method = "bayes-bin", metadata$design retains the normalized imputation prior (prior_surv), completed-data analysis prior (prior_bin), and imputation horizon (end_of_study).

Predictive and analysis priors

For method = "bayes-surv", prior_surv_final is used during interim calculations as well as at the actual final analysis. The two arguments specify different roles, not simply different calendar stages:

Within one interim predictive replicate, first update prior_surv with the observed events and exposure, draw hazards, and generate outstanding outcomes. Then start a fresh analysis posterior using prior_surv_final and the completed dataset's events and exposure. Compare its posterior probability of the alternative with prob_ha. The proportion of replicates that pass is the predictive probability used by Qn, Sn, and Fn.

To incorporate external evidence in prediction while using a weak analysis prior, explicitly supply an informative prior_surv and the chosen weak prior_surv_final. Omitting prior_surv_final uses prior_surv for both roles; the package does not automatically weaken the analysis prior. The predictive prior can still affect the selected sample size and stopping decision, so calibrate the design using both prespecified priors.

This table describes Bayesian survival analysis. For method = "bayes-bin", completed-data success tests at interim and final use prior_bin; prior_surv_final governs only optional final imputation. Frequentist completed-data tests use no analysis prior. evaluate_interim() performs the two interim calculations; use the same prior arguments as in the simulated design.

Examples

interim_data <- data.frame(
  id = 1:6,
  treatment = c(0, 1, 0, 1, 0, 1),
  enrollment = c(0, 1, 2, 3, 4, 5),
  time = c(6, 5, 4, 3, 2, 1),
  event = c(1, 0, 0, 1, 0, 0),
  status = c("event", "pending", "pending", "event", "pending", "pending")
)

evaluate_interim(
  data = interim_data,
  data_cut = 6,
  look = 1,
  N_total = 10,
  end_of_study = 12,
  rand_ratio = c(control = 1, treatment = 1),
  alternative = "less",
  N_impute = 5,
  seed = 2026
)

Goldilocks Bayesian adaptive trial designs

Description

Implements the Goldilocks Bayesian adaptive design proposed by Broglio et al. (2014) for single-arm and randomized 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", "rmst", or "bayes-surv") or reduce complete follow-up to event status at a fixed endpoint time (method = "riskdiff-wald", "riskdiff-fm", 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, stopping-accrual-for-expected-success, and immediate-success rules. Identical points are aggregated; point size indicates their frequency.

Usage

plot_sim_decisions(x)

Arguments

Details

The upper shaded region represents ppp_stop_now > immediate_success_threshold; the middle region represents ⁠success_threshold < ppp_stop_now <= immediate_success_threshold⁠. Futility applies below or at the expected-success threshold when the maximum-sample-size prediction is strictly below its threshold. A legacy trace without immediate_success_threshold is treated as if that threshold were one, which disables immediate success.

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),
  immediate_success_threshold = 0.99,
  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 immediate success, stopping accrual for 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.

Calculate endpoint event probabilities from piecewise hazards

Description

Calculates the cumulative event probability at a fixed follow-up time for one or more sets of piecewise-constant hazard rates.

Usage

ppwe(hazard, end_of_study, cutpoints = NULL)

Arguments

Details

The cumulative probability depends on interval durations, so the value assigned to an isolated cutpoint has no effect. PWEALL represents its generating hazard with pieces closed on the left and open on the right. When goldilocks assigns realized event times to analysis intervals, it instead uses the survival counting-process convention, open on the left and closed on the right.

Value

A numeric vector of event probabilities in ⁠[0, 1]⁠, with one value for each row of hazard.

Print a calendar-time operating-characteristic summary

Description

Print a calendar-time operating-characteristic summary

Usage

## S3 method for class 'goldilocks_calendar_summary'
print(x, digits = 1, ...)

Arguments

Value

x, invisibly.

Print an externally evaluated interim analysis

Description

Prints the predictive probabilities and decision from an interim analysis returned by evaluate_interim().

Usage

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

Arguments

Value

x, invisibly.

Print a Goldilocks adaptive trial result

Description

Prints the final trial summary and reports how many interim looks were completed for a goldilocks_trial object returned by survival_adapt().

Usage

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

Arguments

Value

The input object, invisibly.

Derive piecewise-constant hazard rates from cumulative event probabilities

Description

Converts cumulative event probabilities at specified follow-up times into the corresponding piecewise-constant hazard rates. This is useful for expressing data-generation assumptions in terms of clinically interpretable event probabilities.

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. This is the PWEALL representation of the continuous generating hazard. Changing the value at an isolated cutpoint does not alter the cumulative probabilities calculated here. When observed event times are assigned to analysis intervals, goldilocks uses ⁠(s_{j-1}, s_j]⁠, matching the survival counting-process convention, so an event at s_j belongs to the interval ending there.

Value

A numeric vector of non-negative hazard rates, with one value for each interval defined by cutpoints and endtime.

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

Draws an event time from a piecewise-exponential distribution conditional on a subject remaining event-free through the observed follow-up time, with optional administrative censoring.

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.

PWEALL represents the generating hazard with pieces closed on the left and open on the right. Its cumulative distribution is continuous at every cutpoint, so this endpoint choice does not affect imputation. For assigning realized event times to analysis intervals, goldilocks uses the survival counting-process convention, open on the left and closed on the right; an event exactly at a cutpoint belongs to the interval ending there.

Value

A data frame with one row per subject and columns time, the imputed event or censoring time, and event, coded 1 for an event and 0 for administrative 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

Simulates event times from a piecewise-exponential distribution, with optional administrative censoring at a fixed follow-up time.

Usage

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

Arguments

Details

PWEALL represents the generating hazard with pieces closed on the left and open on the right. Because the event-time distribution is continuous, the value of the hazard at an isolated cutpoint does not alter the cumulative hazard, distribution, or generated samples. When realized event times are later assigned to analysis intervals, goldilocks follows the survival counting-process convention, open on the left and closed on the right, so an event exactly at a cutpoint belongs to the interval ending there. See pwe_impute() for the conditional sampling details.

Value

A data frame with one row per simulated subject and columns time, the event or censoring time, and event, coded 1 for an event and 0 for administrative 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)

Generate a block-randomized treatment sequence

Description

Generates a randomized treatment assignment sequence for control and treatment arms using a fixed allocation ratio and one or more permitted block sizes.

Usage

randomization(N_total, block = 2, allocation = c(control = 1, treatment = 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 two-arm randomization ratios, but they must be given in integer form. For every value b in block, the required relationship is b %% sum(allocation) == 0; see the equal- and unequal-allocation examples below.

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(control = 2, treatment = 3)
)

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

# 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(control = 1, treatment = 2)
)

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

# randomization() is a two-arm helper; a multi-arm allocation is rejected.
try(randomization(
  N_total = 60,
  block = 6,
  allocation = c(1, 1, 1)
), silent = TRUE)

Simulate complete trial data under piecewise-exponential event rates

Description

Simulates enrollment, treatment allocation, event or censoring times, and loss to follow-up for a single-arm or randomized two-arm trial. Event times follow a piecewise-exponential distribution within each arm.

Usage

sim_comp_data(
  hazard_treatment,
  hazard_control = NULL,
  generation_cutpoints = NULL,
  N_total,
  lambda = 0.3,
  lambda_time = NULL,
  end_of_study,
  block = 2,
  rand_ratio = c(control = 1, treatment = 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 generation_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. generation_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.

PWEALL represents the continuous generating hazard with pieces closed on the left and open on the right. This differs from the package's open-left, closed-right convention for assigning realized times only at the cutpoints themselves, which have probability zero under the continuous model. The cumulative hazard, event-time distribution, and generated simulations are therefore unchanged.

Dropout is independent censoring conditional on treatment arm. For event time T, loss_to_fu is true only when D < \min(T, \tau). Administrative censoring and dropout after an observed event are not counted as loss to follow-up. For example, prop_loss = 0.05 with end_of_study = 12 specifies a 5% dropout CDF at 12 months if the time unit is months; it does not force five losses in a 100-subject trial. Equal dropout probabilities in arms with different event hazards need not yield equal observed dropout proportions.

To express a dropout probability q supplied at a different reference time t_0, use p = 1 - (1-q)^{\tau/t_0} at end_of_study to preserve the same exponential dropout hazard. Treatment discontinuation is not separately modeled and should not be treated as loss to follow-up if endpoint collection continues.

This independent exponential mechanism replaces selection of ⁠ceiling(prop_loss * arm size)⁠ subjects followed by censoring uniformly before each selected subject's potential event or administrative time. Positive prop_loss values therefore change seeded results and design operating characteristics relative to the previous mechanism. Simulations with prop_loss = 0 are unchanged.

Value

A data frame with one row per subject and columns:

• time: Numeric event or censoring time.

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

• event: Numeric event indicator, coded 1 for an event and 0 for right-censoring.

• enrollment: Numeric time of subject enrollment relative to first patient in. The package treats enrollment and randomization as occurring at the same time.

• id: Integer subject identifier.

• loss_to_fu: Logical indicator that dropout occurred before both the event and the administrative follow-up horizon.

Estimate operating characteristics by trial simulation

Description

Repeats survival_adapt() under fixed design and data-generating assumptions, returning trial-level results from which operating characteristics can be estimated.

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(control = 1, treatment = 1),
  prop_loss = 0,
  alternative = "greater",
  h0 = 0,
  Fn = 0.05,
  Sn = 0.9,
  prob_ha = 0.95,
  N_impute = 500,
  N_mcmc = 1000,
  mc_conf_level = 0.95,
  N_trials = 10,
  method = "logrank",
  imputed_final = FALSE,
  empty_interval = c("prior", "propagate", "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,
  generation_cutpoints = cutpoints,
  Qn = 1,
  rmst_tau = end_of_study
)

Arguments

Details

This function is a wrapper for survival_adapt() that repeatedly simulates independent trials under the same design parameters and assumed treatment effect.

To use multiple cores (where available), the argument ncores can be increased from the default of 1. The default backend = "auto" stays sequential for fewer than four trials and otherwise uses no more than one core per two trials. This avoids parallel-processing overhead for small simulation studies. On Unix-like systems parallel trials use forked R processes; on Windows they use PSOCK processes. Set backend explicitly when a particular computational method is required.

Errors raised by an individual survival_adapt() call are isolated so other trials can finish. Failed trials are excluded from sims, recorded in failures with their trial number, error class, and message, and reported together in one warning. If every requested trial fails, sim_trials() stops and attaches the same failure table to the error as failures. With a supplied seed, the original call and failed trial number reproduce the same per-trial random-number stream.

With a supplied seed, each trial receives an independent random-number stream. The resulting trial-level simulations are identical whether they are run sequentially or with a supported parallel method, and the pre-existing R random-number state is restored afterward. With seed = NULL, the current random-number state is used and advanced.

Value

A list containing sims, a data frame with one row per successfully simulated trial; failures, a data frame with columns trial, error_class, and message; 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. Per-trial calendar-time metrics are always retained in sims; traces additionally retain calendar time and active follow-up at each look. See survival_adapt() for details of the summary and trace columns, and summarise_calendar_time() for wide operating-characteristic tables. The returned object also retains the evaluated decision_design and resolved prior_design attributes from survival_adapt(). An rng_metadata attribute records the random-number generator, computational method, and seed policy. A parallel_metadata attribute records the requested and actual computational method and number of cores. An arguments attribute contains a named list of all evaluated argument values, including defaults. Its prop_loss element contains a named value for every simulated arm, and its rand_ratio element is stored in control, treatment order for two-arm designs. Its cutpoints and generation_cutpoints elements retain the analysis and data-generation partitions, respectively. For method = "bayes-bin", it also retains the imputation priors (prior_surv and prior_surv_final), completed-data analysis prior (prior_bin), and imputation horizon (end_of_study). The attribute can be saved with saveRDS() and supplied to a later call with do.call(sim_trials, attr(result, "arguments")).

Predictive and analysis priors

For method = "bayes-surv", prior_surv_final is used during interim calculations as well as at the actual final analysis. The two arguments specify different roles, not simply different calendar stages:

Within one interim predictive replicate, first update prior_surv with the observed events and exposure, draw hazards, and generate outstanding outcomes. Then start a fresh analysis posterior using prior_surv_final and the completed dataset's events and exposure. Compare its posterior probability of the alternative with prob_ha. The proportion of replicates that pass is the predictive probability used by Qn, Sn, and Fn.

To incorporate external evidence in prediction while using a weak analysis prior, explicitly supply an informative prior_surv and the chosen weak prior_surv_final. Omitting prior_surv_final uses prior_surv for both roles; the package does not automatically weaken the analysis prior. The predictive prior can still affect the selected sample size and stopping decision, so calibrate the design using both prespecified priors.

This table describes Bayesian survival analysis. For method = "bayes-bin", completed-data success tests at interim and final use prior_bin; prior_surv_final governs only optional final imputation. Frequentist completed-data tests use no analysis prior. evaluate_interim() performs the two interim calculations; use the same prior arguments as in the simulated design.

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(control = 1, treatment = 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 operating characteristics on the calendar-time scale

Description

Summarizes trial duration, follow-up burden, and (when retained) interim-look timing without adding any new simulation arguments. Calendar time is measured from the first patient's enrollment at time zero. analysis_ready_time is the time at which the last enrolled subject's observed event or censoring becomes available; it does not include external data-cleaning or database-lock delays.

Pass a complete result from sim_trials() to retain requested, analyzed, and failed simulation counts. Interim timing requires simulations run with return_trace = TRUE. A simulation data frame can also be supplied, but only the trial-duration table can then be calculated.

Usage

summarise_calendar_time(data)

Arguments

Value

An object of class goldilocks_calendar_summary, containing two wide data frames:

• trial_duration: one row per scenario and stopping reason, plus an overall row. It reports simulation denominators, stopping counts, sample size, accrual-stop time, analysis-ready time, planned completion time, total person-time under follow-up, and peak concurrent follow-up. Newer results also report the time at which the terminal decision was made.

• interim_timing: one row per scenario and interim look. It reports how often the look was reached, its calendar time, and the number of subjects actively under follow-up. This table has zero rows when traces are not available.

Continuous quantities are reported as means, Monte Carlo standard errors, and the 10th, 50th, and 90th percentiles. Percentages use the requested number of simulations as their denominator, so excluded failures remain visible.

Estimate operating characteristics from trial simulations

Description

Estimates success, stopping, and sample-size operating characteristics from a collection of simulated trials and quantifies their Monte Carlo uncertainty.

Usage

summarise_sims(data, max_mcse = NULL)

Arguments

Value

A 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 declared immediate success, stopped accrual for expected success, stopped for futility, or went to the maximum sample size. stop_success retains its historical meaning of stopping accrual for expected success; stop_any_success combines both success-stopping decisions. The average stopping sample size (and standard deviation) are also recorded. The proportion of trials that stopped accrual for expected success, yet went on to fail, is also reported. Each probability and mean is accompanied by its Monte Carlo standard error and 95% Monte Carlo confidence limits, with columns ending in ⁠_mcse⁠, ⁠_mc_lower⁠, and ⁠_mc_upper⁠. Probability intervals use the Wilson method; the mean sample-size interval uses a t distribution.

These intervals describe uncertainty from using a finite number of simulated trials under fixed design and data-generating assumptions. They are not clinical confidence intervals, treatment-effect intervals, or measures of model uncertainty.

The output always reports n_used, the number of successfully analyzed simulations used by the operating-characteristic estimands. When complete sim_trials() results are supplied, it also reports requested, analyzed, and failed counts, the failure rate, computational method, and seed. For a raw simulation data frame, requested and failed counts are unknown and are reported as NA. Details of the call, random-number generation, parallel computation, timing, failures, and design assumptions are retained in the simulation_metadata attribute.

Summarize an interim decision path

Description

Creates a compact one-row summary of the final interim look and stopping decision. Pass a 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

Simulates one single-arm or randomized two-arm trial under a Goldilocks sample-size design. At each planned interim look, posterior predictive probabilities determine whether to declare immediate success, stop accrual for expected success, stop for futility, or continue toward the maximum sample size.

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(control = 1, treatment = 1),
  prop_loss = 0,
  alternative = "greater",
  h0 = 0,
  Fn = 0.05,
  Sn = 0.9,
  prob_ha = 0.95,
  N_impute = 500,
  N_mcmc = 1000,
  mc_conf_level = 0.95,
  empty_interval = c("prior", "propagate", "error"),
  method = "logrank",
  imputed_final = FALSE,
  return_trace = FALSE,
  binary_imputation = c("event-time", "bernoulli"),
  prior_surv_final = prior_surv,
  generation_cutpoints = cutpoints,
  Qn = 1,
  rmst_tau = end_of_study
)

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 satisfy the completed-data success criterion. At each interim analysis this proportion is first compared to the corresponding element of Qn. If it is strictly greater than Qn, the trial stops and declares immediate success. Otherwise, if it is strictly greater than Sn, accrual/enrollment is suspended and the outstanding follow-up is allowed to complete before conducting the pre-specified final analysis.

2. The posterior predictive probability of success at the maximum sample size. This is calculated as the proportion of imputed datasets at the maximum sample size that satisfy the completed-data success criterion. It is compared to the corresponding element of Fn, and if it is below 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, each interim look has four possible decisions, applied in this order:

1. Stop and declare immediate success when P_{n,l} > Q_l.

2. Stop accruing for expected success and follow when S_l < P_{n,l} \le Q_l.

3. Stop for futility when P_{n_{max},l} < F_l.

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

The following completed-data analysis methods are available at interim and final analyses:

• 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 and outcomes are missing, the Cox model is fitted 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.

• Restricted mean survival time (method = "rmst"). Estimates the area under each Kaplan-Meier survival curve from zero through the prespecified rmst_tau. The effect is treatment minus control RMST, measured in time units. The Wald test uses the sum of the independent arm-specific Greenwood variances and reports 1 - P. All three alternatives and nonzero margins are supported. For an adverse event, a positive difference favors treatment. Observed censored subjects remain in the analysis. An arm whose last follow-up precedes rmst_tau with positive estimated survival makes the analysis non-estimable; the horizon is never reduced automatically. A curve reaching zero earlier is allowed. The test requires positive total variance, including after Rubin pooling for imputed_final = TRUE. This unadjusted two-arm analysis does not assume proportional hazards. Interim prediction and final imputation still depend on the specified piecewise-exponential model. See vignette("rmst", package = "goldilocks").

• Bayesian difference in cumulative event probability (method = "bayes-surv"). Each imputed dataset is used to update the conjugate Gamma analysis prior prior_surv_final, at both interim looks and 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. In a single-arm study, inference concerns the treatment-arm event probability. In a two-arm study, the independent arm-specific posteriors define the posterior distribution of the treatment-minus-control difference. The reported posterior probability is determined by alternative and 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, "prior", leaves an empty interval prior-driven, making the absence of interval data explicit. The legacy "propagate" option borrows sufficient statistics from the nearest non-empty interval within the same treatment arm. It is operationally stable but statistically consequential because adjacent observed data then inform the empty interval's posterior. "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.

  The normal approximation matches posterior moments but can misrepresent tail probabilities when a Beta posterior is skewed, particularly near event-probability boundaries of 0 or 1. This can change success and interim stopping decisions. Use bin_method = "quadrature" or sufficiently precise "mc" when the approximation is unsuitable; increasing N_mcmc does not correct normal-approximation error.

  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-wald" or "riskdiff-fm"). Each complete or predictively imputed dataset is reduced to binary event outcomes at end_of_study. The estimated treatment effect is p_\textrm{treatment} - p_\textrm{control}. "riskdiff-wald" uses the observed arm risks in an unpooled Wald variance. "riskdiff-fm" instead uses maximum likelihood arm risks constrained by the null difference h0 in a Farrington-Manning score variance. The latter remains defined for common sparse tables, including equal-arm all-zero and all-one outcomes. Both methods report 1 - P in post_prob_ha and support all three alternatives. Because they require complete binary outcomes, lost-to-follow-up subjects are excluded when imputed_final = FALSE. With genuinely missing final outcomes, only "riskdiff-wald" supports imputed_final = TRUE: estimates and variances from at least two completed datasets are pooled using Rubin's rules, requiring positive total variance. FM final imputation is rejected rather than converted to a Wald test. With complete final outcomes, either flag preserves the selected test.

• 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, RMST, and the frequentist risk-difference Wald analysis, or on the non-imputed dataset. Interim prediction completes outcomes that are not yet observed, whereas final imputation applies only to subjects right-censored because of loss to follow-up before end_of_study. Design evaluations should prespecify whether the final analysis imputes these outcomes and assess sensitivity to that choice, particularly when appreciable attrition is expected. Simulated dropout is independent of event time within each arm, as described under prop_loss. This supports survival analyses with right-censoring, but does not make complete-case binary analysis unbiased: early events can be observed before dropout, whereas later endpoint outcomes can be missing. Binary designs with dropout should assess model-based final imputation and its assumptions. For Cox regression and RMST the final estimates and variances are pooled with Rubin's rules. It cannot be used with method = "logrank".

When imputation is involved, either at interim analyses or through imputed_final = TRUE, the package uses a two-stage impute-then-analyze procedure. First, the piecewise-exponential model is fitted to the observed time-to-event data and used to complete pending outcomes. Second, each completed dataset is analyzed using the model selected by method.

For method = "bayes-bin", these are deliberately separate models. The imputation model has piecewise hazards with Gamma prior prior_surv (or prior_surv_final during final imputation), whereas the completed-data analysis has an event probability at end_of_study with Beta prior prior_bin. The Beta prior is not derived from the Gamma prior, and the two stages are not one joint Bayesian model. Both prior specifications can therefore affect predictive decisions when outcomes require imputation. The "event-time" and "bernoulli" binary-imputation options use the same piecewise-exponential prediction model and do not change this separation.

For method = "bayes-surv", the second analysis instead forms a fresh piecewise-exponential posterior from the completed data and prior_surv_final. At interim looks the first stage uses prior_surv, allowing predictive borrowing to differ from the final success criterion. For frequentist methods ("logrank", "cox", "rmst", "riskdiff-wald", and "riskdiff-fm"), each completed dataset uses a standard test rather than a posterior. Imputed Cox, RMST, and risk-difference final analyses pool estimates and variances using Rubin's rules.

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 one-row data frame containing the evaluated design and final trial results, 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. For RMST this is treatment minus control RMST through rmst_tau, in time units. 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. It is NA after an immediate-success decision, for which no later analysis is required, or when the optional diagnostic analysis after binding futility cannot be computed.

• 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 RMST and risk-difference analyses. For non-imputed frequentist analyses it is 1 - P from the corresponding test. It is NA after an immediate-success decision, for which no later analysis is required, or when the optional diagnostic analysis after binding futility cannot be computed.

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

• stop_immediate_success: Logical indicator of whether the trial stopped and declared immediate success at an interim look.

• stop_expected_success: Logical indicator of whether the trial stopped accruing for expected success and continued planned follow-up.

• trial_success: Logical indicator of the trial's official success outcome. An immediate-success decision is final at the interim look; binding futility is a final failure.

• stopping_reason: One of "immediate_success", "expected_success", "futility", or "maximum_sample_size".

• decision_time: Calendar time at which the trial decision becomes final. This is the interim look time for immediate success or futility, and the analysis-ready time otherwise.

• accrual_stop_time: Calendar time of the last enrollment in the trial.

• analysis_ready_time: Calendar time at which the last enrolled subject's observed event or censoring becomes available. This excludes external data-cleaning and database-lock delays.

• planned_completion_time: Calendar time at which the last enrolled subject would complete the full planned follow-up.

• followup_person_time: Sum of observed follow-up times across enrolled subjects.

• peak_active_followup: Largest number of enrolled subjects concurrently under follow-up.

Calendar time is measured from the first patient's enrollment at time zero. Times use the same units as lambda_time, cutpoints, generation_cutpoints, and end_of_study.

The returned object has a decision_design attribute containing interim_look, Fn, Sn, Qn, and the Monte Carlo settings. Thresholds in this information are stored as one value per interim look (and have length zero when no interim looks are planned). A prior_design attribute contains the resolved Gamma shape, rate, mean hazard, and standard deviation for every stage, arm, and interval.

Both return forms have an arguments attribute containing a named list of the evaluated argument values, including defaults. It can be saved with saveRDS() and supplied to a later call with do.call(survival_adapt, attr(result, "arguments")). Its prop_loss element contains a named value for every simulated arm. Its rand_ratio element is stored in control, treatment order for two-arm designs. Its cutpoints and generation_cutpoints elements retain the analysis and data-generation partitions, respectively. For method = "bayes-bin", this metadata explicitly retains the imputation priors (prior_surv and prior_surv_final), completed-data analysis prior (prior_bin), and imputation horizon (end_of_study). The separate prior_design attribute gives the resolved Gamma parameters by stage, arm, and interval.

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. prior_diagnostics contains the resolved interim and final priors. posterior_diagnostics reports observed and effective sufficient statistics and conjugate posterior parameters by completed look, arm, and interval. The trace records calendar time, the number of subjects actively under follow-up, enrollment and observed events by arm, predictive probabilities, diagnostic Monte Carlo standard errors and exact bounds, draw counts, thresholds, the decision and reason, empty-interval fallback diagnostics, and warnings raised during that look. It deliberately excludes imputed data sets and posterior draws to keep the output compact.

Predictive and analysis priors

For method = "bayes-surv", prior_surv_final is used during interim calculations as well as at the actual final analysis. The two arguments specify different roles, not simply different calendar stages:

Within one interim predictive replicate, first update prior_surv with the observed events and exposure, draw hazards, and generate outstanding outcomes. Then start a fresh analysis posterior using prior_surv_final and the completed dataset's events and exposure. Compare its posterior probability of the alternative with prob_ha. The proportion of replicates that pass is the predictive probability used by Qn, Sn, and Fn.

To incorporate external evidence in prediction while using a weak analysis prior, explicitly supply an informative prior_surv and the chosen weak prior_surv_final. Omitting prior_surv_final uses prior_surv for both roles; the package does not automatically weaken the analysis prior. The predictive prior can still affect the selected sample size and stopping decision, so calibrate the design using both prespecified priors.

This table describes Bayesian survival analysis. For method = "bayes-bin", completed-data success tests at interim and final use prior_bin; prior_surv_final governs only optional final imputation. Frequentist completed-data tests use no analysis prior. evaluate_interim() performs the two interim calculations; use the same prior arguments as in the simulated design.

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(control = 1, treatment = 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")