Package 'goldilocks'

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

lambda	finite positive enrollment rates per unit time. Supply one rate for each interval defined by lambda_time, so length(lambda) must equal length(lambda_time) + 1.
N_total	positive integer total sample size.
lambda_time	NULL, or a finite, positive, strictly increasing vector of interior times at which the enrollment rate changes. The initial boundary at time zero is implicit and must not be supplied. Use NULL for a constant enrollment rate.

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 time-to-event endpoint trials, both one- and two-arm, with an underlying piecewise exponential hazard model. The method can be used for a confirmatory trial to select a trial's sample size based on accumulating data. During accrual, frequent sample size selection analyses are made and predictive probabilities are used to determine whether the current sample size is sufficient or whether continuing accrual would be futile. The algorithm explicitly accounts for complete follow-up of all patients before the primary analysis is conducted. 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 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

x	A simulation result returned by sim_trials() with return_trace = TRUE, or its traces data frame.

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

x	A data frame returned by summarise_sims(), with one row per simulation scenario.
effect	Numeric treatment-effect values corresponding to the rows of x, or a single character string naming a numeric column in x.
xlab	Character label for the treatment-effect axis.

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

x

A simulation result data frame or the list returned by sim_trials.

type

Character string specifying the percentages to plot. "marginal" shows the percentage of all simulated trials ending at each sample size; its bars sum to 100 percent across sample sizes. "conditional" shows the percentage stopping at each look among trials still active at the start of that look. "cumulative" shows the status of all simulated trials after each look; every bar sums to 100 percent and includes trials continuing to the next look. "flowchart" starts with all simulated trials and branches at each look into futility, continued enrollment, and expected-success nodes labelled with trial counts.

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

x	A goldilocks_trial object or an interim trace data frame.

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

hazard	matrix. A matrix of hazard rate parameters with number of columns one greater than the length of the cutpoints vector. The number of rows can be anything, and is typically dictated by the number of MCMC draws.
end_of_study	finite positive time at which the cumulative event probability is evaluated. It must be greater than every cutpoint.
cutpoints	finite, positive, strictly increasing vector of interior times at which the hazard rate changes. The number of hazard rates must be one greater than the number of cutpoints. Use NULL for a constant hazard.

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

x	A goldilocks_trial object.
...	Additional arguments passed to print.data.frame.

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

probs	vector. Cumulative event probabilities at each cutpoint and at endtime, in that order. Its length must be one greater than the number of cutpoints. With no cutpoints, supply a single probability at endtime.
cutpoints	finite, positive, strictly increasing vector of interior times at which the baseline hazard changes. Default is NULL, which corresponds to a simple (non-piecewise) exponential model.
endtime	scalar. Time at which final element in probs corresponds to. Typically this would be the study endpoint time.

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

time	vector. The observed time for patient that have had no event or passed maxtime.
hazard	vector. Finite non-negative constant hazard rates for exponential failures. If the final rate is zero, maxtime must be supplied so that subjects without an event can be administratively censored.
cutpoints	finite, positive, strictly increasing vector of interior times at which the hazard rate changes. The number of hazard rates must be one greater than the number of cutpoints. Use NULL for a constant hazard.
maxtime	scalar. Optional administrative censoring time. When supplied, it must be later than every cutpoint.

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

n	integer. The number of random samples to generate. Default is n = 1.
hazard	vector. Finite non-negative constant hazard rates for exponential failures. If the final rate is zero, maxtime must be supplied so that subjects without an event can be administratively censored.
cutpoints	finite, positive, strictly increasing vector of interior times at which the hazard rate changes. The number of hazard rates must be one greater than the number of cutpoints. Use NULL for a constant hazard.
maxtime	scalar. Optional administrative censoring time. When supplied, it must be later than every cutpoint.

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

N_total	integer. Total sample size for randomization allocation.
block	vector. Block size for randomization. Note that it needs to be a multiple of the sum of allocation.
allocation	vector. The randomization allocation in the order c(control, treatment).

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

hazard_treatment	vector. Finite non-negative constant hazard rates under the treatment arm.
hazard_control	vector. Finite non-negative constant hazard rates under the control arm.
cutpoints	finite, positive, strictly increasing interior times at which the baseline hazard changes. The number of hazards for each arm must be one greater than the number of cutpoints. Default is NULL, which corresponds to a simple (non-piecewise) exponential model.
N_total	integer. Maximum sample size allowable
lambda	finite positive enrollment rates per unit time. Supply one rate for each interval defined by lambda_time. See enrollment() for the precise continuous-time process and time-origin convention.
lambda_time	NULL, or finite, positive, strictly increasing internal times at which the enrollment rate changes. The initial boundary at zero is implicit, so length(lambda) must equal length(lambda_time) + 1.
end_of_study	finite study endpoint, strictly greater than the last cutpoint.
block	scalar. Block size for generating the randomization schedule.
rand_ratio	vector. Randomization allocation for the ratio of control to treatment. Integer values mapping the size of the block. See randomization() for more details.
prop_loss	scalar. Overall proportion of subjects lost to follow-up. Subjects are selected at random for LTFU regardless of treatment assignment or event status. Each LTFU subject's observed time is drawn from a Uniform(0, t) distribution, where t is their potential event or censoring time. Since the LTFU time is always less than t, the event has not yet occurred at dropout and the subject is right-censored. Defaults to zero.

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 = c(0.1, 0.1),
  bin_prior = 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")
)

Arguments

hazard_treatment	vector. Finite non-negative constant hazard rates under the treatment arm.
hazard_control	vector. Finite non-negative constant hazard rates under the control arm.
cutpoints	finite, positive, strictly increasing interior times at which the baseline hazard changes. The number of hazards for each arm must be one greater than the number of cutpoints. Default is NULL, which corresponds to a simple (non-piecewise) exponential model.
N_total	integer. Maximum sample size allowable
lambda	finite positive enrollment rates per unit time. Supply one rate for each interval defined by lambda_time. See enrollment() for the precise continuous-time process and time-origin convention.
lambda_time	NULL, or finite, positive, strictly increasing internal times at which the enrollment rate changes. The initial boundary at zero is implicit, so length(lambda) must equal length(lambda_time) + 1.
interim_look	vector. Sample size for each interim look. Note: the maximum sample size should not be included. For two-arm designs, each interim look must be at least the (largest) block size (see block), ensuring both treatment groups are present at every interim analysis; a smaller look could enroll subjects from one treatment group only, leaving the interim posterior undefined for the missing group.
end_of_study	finite study endpoint, strictly greater than the last cutpoint.
prior	vector. The prior distributions for the piecewise hazard rate parameters are each $Gamma(a_0, b_0)$, where $a_0$ is the shape parameter and $b_0$ is the rate parameter (i.e., the inverse of the scale). This follows R's stats::rgamma() parameterization. The same prior is applied to all piecewise intervals and to both treatment groups. The default non-informative prior distribution used is Gamma(0.1, 0.1), which is specified by setting prior = c(0.1, 0.1).
bin_prior	vector. Prior distribution for the event probability when method = "bayes-bin". The two values are the shape parameters of the Beta(a, b) prior. The same prior is applied to both treatment arms.
bin_method	character. Method used to calculate the posterior probability for method = "bayes-bin", must be one of "mc" (Monte Carlo sampling), "normal" (normal approximation), or "quadrature" (numerical integration). The default is "mc".
block	scalar. Block size for generating the randomization schedule.
rand_ratio	vector. Randomization allocation for the ratio of control to treatment. Integer values mapping the size of the block. See randomization() for more details.
prop_loss	scalar. Overall proportion of subjects lost to follow-up. Subjects are selected at random for LTFU regardless of treatment assignment or event status. Each LTFU subject's observed time is drawn from a Uniform(0, t) distribution, where t is their potential event or censoring time. Since the LTFU time is always less than t, the event has not yet occurred at dropout and the subject is right-censored. Defaults to zero.
alternative	character. The string specifying the alternative hypothesis, must be one of "greater" (default), "less" or "two.sided". One-sided alternatives ("greater" and "less") are supported for method = "bayes-surv" and method = "bayes-bin". All three options are supported for method = "logrank", method = "cox", and method = "riskdiff". For survival outcomes, "less" corresponds to the treatment arm having a lower cumulative incidence (i.e., treatment is beneficial), and "greater" corresponds to the treatment arm having a higher cumulative incidence.
h0	single finite numeric null hypothesis value or margin. Default is h0 = 0. For Bayesian analyses, h0 must lie in ⁠[0, 1]⁠ for a single-arm design and ⁠[-1, 1]⁠ for a two-arm design. When method = "bayes-surv", h0 is the null value of $p_\textrm{treatment} - p_\textrm{control}$. In a single-arm design, h0 is the external benchmark event probability, often referred to as a performance goal (PG) or objective performance criterion (OPC). When method = "bayes-bin", h0 is the null value of $p_\textrm{treatment} - p_\textrm{control}$ for a two-arm design, or the null event probability for a single-arm design. When method = "cox", h0 is the null log hazard ratio for treatment versus control. Use h0 = 0 for the usual hazard ratio of 1 null, or h0 = log(margin) for a non-inferiority margin specified as a hazard ratio. A Cox non-inferiority test should usually use alternative = "less". When method = "riskdiff", h0 is the null value of $p_\textrm{treatment} - p_\textrm{control}$ and must lie in ⁠[-1, 1]⁠. The argument is ignored for method = "logrank" after its finite-value validation; the usual equal-survival null hypothesis is used.
Fn	vector of values between 0 and 1. Each element is the probability threshold to stop at the $i$-th look early for futility. If there are no interim looks (i.e. interim_look = NULL), then Fn is not used in the simulations or analysis. Set Fn = 0 to disable futility monitoring. The length of Fn should be the same as interim_look, else the values are recycled.
Sn	vector of values between 0 and 1. Each element is the probability threshold to stop at the $i$-th look early for expected success. If there are no interim looks (i.e. interim_look = NULL), then Sn is not used in the simulations or analysis. The length of Sn should be the same as interim_look, else the values are recycled.
prob_ha	scalar value between 0 and 1. Probability threshold of alternative hypothesis.
N_impute	integer. Number of imputations for Monte Carlo simulation of missing data. An imputed Cox or risk-difference final analysis requires at least two.
N_mcmc	integer. Number of posterior samples used by method = "bayes-surv" and by method = "bayes-bin" when bin_method = "mc".
N_trials	integer. Number of trials to simulate.
method	character. For an imputed data set (or the final data set after follow-up is complete), whether the analysis should be a log-rank (method = "logrank") test, Cox proportional hazards regression model Wald test (method = "cox"), a fully-Bayesian piecewise-exponential analysis (method = "bayes-surv"), a Bayesian beta-binomial analysis of complete binary outcomes (method = "bayes-bin"), or a frequentist risk-difference Wald test of complete binary outcomes (method = "riskdiff"). See Details section.
imputed_final	logical. Should the final analysis (after all subjects have been followed-up to the study end) be based on imputed outcomes for subjects who were LTFU (i.e. right-censored with time less than end_of_study)? Default is FALSE, which means that the final analysis incorporates right-censoring. With method = "cox" or method = "riskdiff", setting this to TRUE analyzes each imputed dataset and pools the scalar treatment effects and variances using Rubin's rules; this requires N_impute >= 2. Imputed final analyses remain unavailable for method = "logrank".
empty_interval	character. Policy for empty piecewise-exponential intervals in method = "bayes-surv" posterior calculations. An empty interval is an interval with no exposed subjects in a treatment arm at the analysis time. "propagate" (the default, matching earlier package behavior) copies exposure time and event counts from the nearest non-empty interval in the same treatment arm and emits a warning. "prior" leaves the interval at zero exposure time and zero events, so its posterior is driven only by prior. "error" stops when any empty interval is found.
return_trace	logical. Should the compact interim decision trace from every simulated trial be retained? The default, FALSE, preserves the compact output. When TRUE, the returned list also contains a traces data frame with a trial column linking each trace row to the corresponding row of sims.
ncores	positive integer. Number of cores to use for parallel processing. Defaults to 1L (serial execution).
backend	character. Parallel backend. "auto" (the default) uses serial execution for ncores = 1, the existing fork backend on Unix-like platforms, and a PSOCK cluster on Windows. "fork", "psock", and "sequential" select a backend explicitly.
seed	optional integer. Seed used to generate independent per-trial "L'Ecuyer-CMRG" random-number streams. The default, NULL, does not reset the global RNG state, preserving the usual unseeded simulation behavior.
binary_imputation	character. Predictive imputation approach for method = "bayes-bin" or method = "riskdiff". "event-time" (the default) draws a conditional piecewise-exponential event time and reduces it to event status at end_of_study. "bernoulli" draws the endpoint status directly from its conditional event probability. This argument is ignored for time-to-event analysis methods.

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 = 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

data	list (of data frames) or a single data frame. If summarizing a single run of simulations, data will be a data.frame object returned from survival_adapt(). If summarizing multiple simulation scenarios, data will be a list object, with each element being a data.frame object.

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

x	A goldilocks_trial object or an interim trace data frame.

Value

A one-row data frame.

---

Simulate and execute a single adaptive clinical trial design with a time-to-event endpoint

Description

Simulate and execute a single adaptive clinical trial design with a time-to-event endpoint

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 = c(0.1, 0.1),
  bin_prior = 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")
)

Arguments

hazard_treatment	vector. Finite non-negative constant hazard rates under the treatment arm.
hazard_control	vector. Finite non-negative constant hazard rates under the control arm.
cutpoints	finite, positive, strictly increasing interior times at which the baseline hazard changes. The number of hazards for each arm must be one greater than the number of cutpoints. Default is NULL, which corresponds to a simple (non-piecewise) exponential model.
N_total	integer. Maximum sample size allowable
lambda	finite positive enrollment rates per unit time. Supply one rate for each interval defined by lambda_time. See enrollment() for the precise continuous-time process and time-origin convention.
lambda_time	NULL, or finite, positive, strictly increasing internal times at which the enrollment rate changes. The initial boundary at zero is implicit, so length(lambda) must equal length(lambda_time) + 1.
interim_look	vector. Sample size for each interim look. Note: the maximum sample size should not be included. For two-arm designs, each interim look must be at least the (largest) block size (see block), ensuring both treatment groups are present at every interim analysis; a smaller look could enroll subjects from one treatment group only, leaving the interim posterior undefined for the missing group.
end_of_study	finite study endpoint, strictly greater than the last cutpoint.
prior	vector. The prior distributions for the piecewise hazard rate parameters are each $Gamma(a_0, b_0)$, where $a_0$ is the shape parameter and $b_0$ is the rate parameter (i.e., the inverse of the scale). This follows R's stats::rgamma() parameterization. The same prior is applied to all piecewise intervals and to both treatment groups. The default non-informative prior distribution used is Gamma(0.1, 0.1), which is specified by setting prior = c(0.1, 0.1).
bin_prior	vector. Prior distribution for the event probability when method = "bayes-bin". The two values are the shape parameters of the Beta(a, b) prior. The same prior is applied to both treatment arms.
bin_method	character. Method used to calculate the posterior probability for method = "bayes-bin", must be one of "mc" (Monte Carlo sampling), "normal" (normal approximation), or "quadrature" (numerical integration). The default is "mc".
block	scalar. Block size for generating the randomization schedule.
rand_ratio	vector. Randomization allocation for the ratio of control to treatment. Integer values mapping the size of the block. See randomization() for more details.
prop_loss	scalar. Overall proportion of subjects lost to follow-up. Subjects are selected at random for LTFU regardless of treatment assignment or event status. Each LTFU subject's observed time is drawn from a Uniform(0, t) distribution, where t is their potential event or censoring time. Since the LTFU time is always less than t, the event has not yet occurred at dropout and the subject is right-censored. Defaults to zero.
alternative	character. The string specifying the alternative hypothesis, must be one of "greater" (default), "less" or "two.sided". One-sided alternatives ("greater" and "less") are supported for method = "bayes-surv" and method = "bayes-bin". All three options are supported for method = "logrank", method = "cox", and method = "riskdiff". For survival outcomes, "less" corresponds to the treatment arm having a lower cumulative incidence (i.e., treatment is beneficial), and "greater" corresponds to the treatment arm having a higher cumulative incidence.
h0	single finite numeric null hypothesis value or margin. Default is h0 = 0. For Bayesian analyses, h0 must lie in ⁠[0, 1]⁠ for a single-arm design and ⁠[-1, 1]⁠ for a two-arm design. When method = "bayes-surv", h0 is the null value of $p_\textrm{treatment} - p_\textrm{control}$. In a single-arm design, h0 is the external benchmark event probability, often referred to as a performance goal (PG) or objective performance criterion (OPC). When method = "bayes-bin", h0 is the null value of $p_\textrm{treatment} - p_\textrm{control}$ for a two-arm design, or the null event probability for a single-arm design. When method = "cox", h0 is the null log hazard ratio for treatment versus control. Use h0 = 0 for the usual hazard ratio of 1 null, or h0 = log(margin) for a non-inferiority margin specified as a hazard ratio. A Cox non-inferiority test should usually use alternative = "less". When method = "riskdiff", h0 is the null value of $p_\textrm{treatment} - p_\textrm{control}$ and must lie in ⁠[-1, 1]⁠. The argument is ignored for method = "logrank" after its finite-value validation; the usual equal-survival null hypothesis is used.
Fn	vector of values between 0 and 1. Each element is the probability threshold to stop at the $i$-th look early for futility. If there are no interim looks (i.e. interim_look = NULL), then Fn is not used in the simulations or analysis. Set Fn = 0 to disable futility monitoring. The length of Fn should be the same as interim_look, else the values are recycled.
Sn	vector of values between 0 and 1. Each element is the probability threshold to stop at the $i$-th look early for expected success. If there are no interim looks (i.e. interim_look = NULL), then Sn is not used in the simulations or analysis. The length of Sn should be the same as interim_look, else the values are recycled.
prob_ha	scalar value between 0 and 1. Probability threshold of alternative hypothesis.
N_impute	integer. Number of imputations for Monte Carlo simulation of missing data. An imputed Cox or risk-difference final analysis requires at least two.
N_mcmc	integer. Number of posterior samples used by method = "bayes-surv" and by method = "bayes-bin" when bin_method = "mc".
empty_interval	character. Policy for empty piecewise-exponential intervals in method = "bayes-surv" posterior calculations. An empty interval is an interval with no exposed subjects in a treatment arm at the analysis time. "propagate" (the default, matching earlier package behavior) copies exposure time and event counts from the nearest non-empty interval in the same treatment arm and emits a warning. "prior" leaves the interval at zero exposure time and zero events, so its posterior is driven only by prior. "error" stops when any empty interval is found.
method	character. For an imputed data set (or the final data set after follow-up is complete), whether the analysis should be a log-rank (method = "logrank") test, Cox proportional hazards regression model Wald test (method = "cox"), a fully-Bayesian piecewise-exponential analysis (method = "bayes-surv"), a Bayesian beta-binomial analysis of complete binary outcomes (method = "bayes-bin"), or a frequentist risk-difference Wald test of complete binary outcomes (method = "riskdiff"). See Details section.
imputed_final	logical. Should the final analysis (after all subjects have been followed-up to the study end) be based on imputed outcomes for subjects who were LTFU (i.e. right-censored with time less than end_of_study)? Default is FALSE, which means that the final analysis incorporates right-censoring. With method = "cox" or method = "riskdiff", setting this to TRUE analyzes each imputed dataset and pools the scalar treatment effects and variances using Rubin's rules; this requires N_impute >= 2. Imputed final analyses remain unavailable for method = "logrank".
return_trace	logical. Should the interim decision path be returned in addition to the usual final summary? The default, FALSE, returns the historical one-row data frame. When TRUE, the result is a goldilocks_trial object with summary, trace, and call elements.
binary_imputation	character. Predictive imputation approach for method = "bayes-bin" or method = "riskdiff". "event-time" (the default) draws a conditional piecewise-exponential event time and reduces it to event status at end_of_study. "bernoulli" draws the endpoint status directly from its conditional event probability. This argument is ignored for time-to-event analysis methods.

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), 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 bin_prior, 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 = 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")

---