simtrial

Clinical trial simulation for time-to-event endpoints

Version

1.0.0

Repository

GitHub

CRAN

View on CRAN

Documentation

llms.txt

Documentation

Raw

FUNCTION: as_gt 
TITLE: Convert summary table to a gt object 
DESCRIPTION: 
Convert summary table to a gt object
 
ARGUMENTS: x A object returned by \code{\link[=summary]{summary()}}.
... Additional parameters (not used).
title Title of the gt table.
subtitle Subtitle of the gt table. 
EXAMPLE: 

# Parameters for enrollment
enroll_rampup_duration <- 4 # Duration for enrollment ramp up
enroll_duration <- 16 # Total enrollment duration
enroll_rate <- gsDesign2::define_enroll_rate(
  duration = c(
    enroll_rampup_duration, enroll_duration - enroll_rampup_duration),
 rate = c(10, 30))

# Parameters for treatment effect
delay_effect_duration <- 3 # Delay treatment effect in months
median_ctrl <- 9 # Survival median of the control arm
median_exp <- c(9, 14) # Survival median of the experimental arm
dropout_rate <- 0.001
fail_rate <- gsDesign2::define_fail_rate(
  duration = c(delay_effect_duration, 100),
  fail_rate = log(2) / median_ctrl,
  hr = median_ctrl / median_exp,
  dropout_rate = dropout_rate)

# Other related parameters
alpha <- 0.025 # Type I error
beta <- 0.1 # Type II error
ratio <- 1 # Randomization ratio (experimental:control)

# Build a one-sided group sequential design
design <- gsDesign2::gs_design_ahr(
  enroll_rate = enroll_rate, fail_rate = fail_rate,
  ratio = ratio, alpha = alpha, beta = beta,
  analysis_time = c(12, 24, 36),
  upper = gsDesign2::gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = alpha),
  lower = gsDesign2::gs_b,
  lpar = rep(-Inf, 3))

# Define cuttings of 2 IAs and 1 FA
ia1_cut <- create_cut(target_event_overall = ceiling(design$analysis$event[1]))
ia2_cut <- create_cut(target_event_overall = ceiling(design$analysis$event[2]))
fa_cut <- create_cut(target_event_overall = ceiling(design$analysis$event[3]))

# Run simulations
simulation <- sim_gs_n(
  n_sim = 3,
  sample_size = ceiling(design$analysis$n[3]),
  enroll_rate = design$enroll_rate,
  fail_rate = design$fail_rate,
  test = wlr,
  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
  weight = fh(rho = 0, gamma = 0.5))

# Summarize simulations
simulation |>
 summary(bound = gsDesign::gsDesign(k = 3, test.type = 1, sfu = gsDesign::sfLDOF)$upper$bound) |>
 simtrial::as_gt()

# Summarize simulations and compare with the planned design
simulation |>
  summary(design = design) |>
  simtrial::as_gt()
 


FUNCTION: check_args 
TITLE: Check argument types, length, or dimension 
DESCRIPTION: 
Check argument types, length, or dimension
 
ARGUMENTS: arg An argument to be checked.
type A character vector of candidate argument type.
length A numeric value of argument length or \code{NULL}.
dim A numeric vector of argument dimension or \code{NULL}. 
EXAMPLE: 
\dontrun{
tbl <- as.data.frame(matrix(1:9, nrow = 3))
simtrial:::check_args(arg = tbl, type = c("data.frame"))

vec <- c("a", "b", "c")
simtrial:::check_args(arg = vec, type = c("character"), length = 3)
}
 


FUNCTION: counting_process 
TITLE: Process survival data into counting process format 
DESCRIPTION: 
Produces a data frame that is sorted by stratum and time.
Included in this is only the times at which one or more event occurs.
The output dataset contains stratum, TTE (time-to-event),
at risk count, and count of events at the specified TTE
sorted by stratum and TTE.
 
ARGUMENTS: x A data frame with no missing values and contain variables:
\itemize{
\item \code{stratum}: Stratum.
\item \code{treatment}: Treatment group.
\item \code{tte}: Observed time.
\item \code{event}: Binary event indicator, \code{1} represents event,
\code{0} represents censoring.
}
arm Value in the input \code{treatment} column that indicates
treatment group value. 
EXAMPLE: 
# Example 1
x <- data.frame(
  stratum = c(rep(1, 10), rep(2, 6)),
  treatment = rep(c(1, 1, 0, 0), 4),
  tte = 1:16,
  event = rep(c(0, 1), 8)
)
counting_process(x, arm = 1)

# Example 2
x <- sim_pw_surv(n = 400)
y <- cut_data_by_event(x, 150) |> counting_process(arm = "experimental")
# Weighted logrank test (Z-value and 1-sided p-value)
z <- sum(y$o_minus_e) / sqrt(sum(y$var_o_minus_e))
c(z, pnorm(z))
 


FUNCTION: create_cut 
TITLE: Create a cutting function 
DESCRIPTION: 
Create a cutting function for use with \code{\link[=sim_gs_n]{sim_gs_n()}}
 
ARGUMENTS: ... Arguments passed to \code{\link[=get_analysis_date]{get_analysis_date()}} 
EXAMPLE: 
# Simulate trial data
trial_data <- sim_pw_surv()

# Create a cutting function that applies the following 2 conditions:
# - At least 45 months have passed since the start of the study
# - At least 300 events have occurred
cutting <- create_cut(
  planned_calendar_time = 45,
  target_event_overall = 350
)

# Cut the trial data
cutting(trial_data)
 


FUNCTION: create_test 
TITLE: Create a cutting test function 
DESCRIPTION: 
Create a cutting test function for use with \code{\link[=sim_gs_n]{sim_gs_n()}}
 
ARGUMENTS: test A test function such as \code{\link[=wlr]{wlr()}}, \code{\link[=maxcombo]{maxcombo()}}, or \code{\link[=rmst]{rmst()}}
... Arguments passed to the cutting test function 
EXAMPLE: 
# Simulate trial data
trial_data <- sim_pw_surv()

# Cut after 150 events
trial_data_cut <- cut_data_by_event(trial_data, 150)

# Create a cutting test function that can be used by sim_gs_n()
regular_logrank_test <- create_test(wlr, weight = fh(rho = 0, gamma = 0))

# Test the cutting
regular_logrank_test(trial_data_cut)

# The results are the same as directly calling the function
stopifnot(all.equal(
  regular_logrank_test(trial_data_cut),
  wlr(trial_data_cut, weight = fh(rho = 0, gamma = 0))
))
 


FUNCTION: cut_data_by_date 
TITLE: Cut a dataset for analysis at a specified date 
DESCRIPTION: 
Cut a dataset for analysis at a specified date
 
ARGUMENTS: x A time-to-event dataset, for example, generated by \code{\link[=sim_pw_surv]{sim_pw_surv()}}.
cut_date Date relative to start of randomization
(\code{cte} from input dataset) at which dataset is to be cut off for analysis. 
EXAMPLE: 
# Use default enrollment and event rates and
# cut at calendar time 5 after start of randomization
sim_pw_surv(n = 20) |> cut_data_by_date(5)
 


FUNCTION: cut_data_by_event 
TITLE: Cut a dataset for analysis at a specified event count 
DESCRIPTION: 
Takes a time-to-event data set and cuts the data at which an
event count is reached.
 
ARGUMENTS: x A time-to-event dataset, for example, generated by \code{\link[=sim_pw_surv]{sim_pw_surv()}}.
event Event count at which data cutoff is to be made. 
EXAMPLE: 
# Use default enrollment and event rates at cut at 100 events
x <- sim_pw_surv(n = 200) |> cut_data_by_event(100)
table(x$event, x$treatment)
 


FUNCTION: early_zero 
TITLE: Zero early weighting function 
DESCRIPTION: 
Zero early weighting function
 
ARGUMENTS: early_period The initial delay period where weights increase;
after this, weights are constant at the final weight in the delay period.
fail_rate Failure rate 
EXAMPLE: 
\dontshow{if (requireNamespace("gsDesign2", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)({ # examplesIf}
library(gsDesign2)

# Example 1: Unstratified ----
sim_pw_surv(n = 200) |>
  cut_data_by_event(125) |>
  wlr(weight = early_zero(early_period = 2))

# Example 2: Stratified ----
n <- 500
# Two strata
stratum <- c("Biomarker-positive", "Biomarker-negative")
prevalence_ratio <- c(0.6, 0.4)

# Enrollment rate
enroll_rate <- define_enroll_rate(
  stratum = rep(stratum, each = 2),
  duration = c(2, 10, 2, 10),
  rate = c(c(1, 4) * prevalence_ratio[1], c(1, 4) * prevalence_ratio[2])
)
enroll_rate$rate <- enroll_rate$rate * n / sum(enroll_rate$duration * enroll_rate$rate)

# Failure rate
med_pos <- 10 # Median of the biomarker positive population
med_neg <- 8 # Median of the biomarker negative population
hr_pos <- c(1, 0.7) # Hazard ratio of the biomarker positive population
hr_neg <- c(1, 0.8) # Hazard ratio of the biomarker negative population
fail_rate <- define_fail_rate(
  stratum = rep(stratum, each = 2),
  duration = c(3, 1000, 4, 1000),
  fail_rate = c(log(2) / c(med_pos, med_pos, med_neg, med_neg)),
  hr = c(hr_pos, hr_neg),
  dropout_rate = 0.01
)

# Simulate data
temp <- to_sim_pw_surv(fail_rate) # Convert the failure rate
set.seed(2023)

sim_pw_surv(
  n = n, # Sample size
  # Stratified design with prevalence ratio of 6:4
  stratum = data.frame(stratum = stratum, p = prevalence_ratio),
  # Randomization ratio
  block = c("control", "control", "experimental", "experimental"),
  enroll_rate = enroll_rate, # Enrollment rate
  fail_rate = temp$fail_rate, # Failure rate
  dropout_rate = temp$dropout_rate # Dropout rate
) |>
  cut_data_by_event(125) |>
  wlr(weight = early_zero(early_period = 2, fail_rate = fail_rate))
\dontshow{}) # examplesIf}
 


FUNCTION: ex1_delayed_effect 
TITLE: Time-to-event data example 1 for non-proportional hazards working group 
DESCRIPTION: 
Survival objects reverse-engineered datasets from published Kaplan-Meier
curves.
Individual trials are de-identified since the data are only
approximations of the actual data.
Data are intended to evaluate methods and designs for trials where
non-proportional hazards may be anticipated for outcome data.
 
ARGUMENTS:  
EXAMPLE: 
library(survival)

data(ex1_delayed_effect)
km1 <- with(ex1_delayed_effect, survfit(Surv(month, evntd) ~ trt))
km1
plot(km1)
with(subset(ex1_delayed_effect, trt == 1), survfit(Surv(month, evntd) ~ trt))
with(subset(ex1_delayed_effect, trt == 0), survfit(Surv(month, evntd) ~ trt))
 


FUNCTION: ex2_delayed_effect 
TITLE: Time-to-event data example 2 for non-proportional hazards working group 
DESCRIPTION: 
Survival objects reverse-engineered datasets from published Kaplan-Meier
curves.
Individual trials are de-identified since the data are only
approximations of the actual data.
Data are intended to evaluate methods and designs for trials where
non-proportional hazards may be anticipated for outcome data.
 
ARGUMENTS:  
EXAMPLE: 
library(survival)

data(ex2_delayed_effect)
km1 <- with(ex2_delayed_effect, survfit(Surv(month, evntd) ~ trt))
km1
plot(km1)
with(subset(ex2_delayed_effect, trt == 1), survfit(Surv(month, evntd) ~ trt))
with(subset(ex2_delayed_effect, trt == 0), survfit(Surv(month, evntd) ~ trt))
 


FUNCTION: ex3_cure_with_ph 
TITLE: Time-to-event data example 3 for non-proportional hazards working group 
DESCRIPTION: 
Survival objects reverse-engineered datasets from published Kaplan-Meier
curves.
Individual trials are de-identified since the data are only
approximations of the actual data.
Data are intended to evaluate methods and designs for trials where
non-proportional hazards may be anticipated for outcome data.
 
ARGUMENTS:  
EXAMPLE: 
library(survival)

data(ex3_cure_with_ph)
km1 <- with(ex3_cure_with_ph, survfit(Surv(month, evntd) ~ trt))
km1
plot(km1)
 


FUNCTION: ex4_belly 
TITLE: Time-to-event data example 4 for non-proportional hazards working group 
DESCRIPTION: 
Survival objects reverse-engineered datasets from published Kaplan-Meier
curves.
Individual trials are de-identified since the data are only
approximations of the actual data.
Data are intended to evaluate methods and designs for trials where
non-proportional hazards may be anticipated for outcome data.
 
ARGUMENTS:  
EXAMPLE: 
library(survival)

data(ex4_belly)
km1 <- with(ex4_belly, survfit(Surv(month, evntd) ~ trt))
km1
plot(km1)
 


FUNCTION: ex5_widening 
TITLE: Time-to-event data example 5 for non-proportional hazards working group 
DESCRIPTION: 
Survival objects reverse-engineered datasets from published Kaplan-Meier
curves.
Individual trials are de-identified since the data are only
approximations of the actual data.
Data are intended to evaluate methods and designs for trials where
non-proportional hazards may be anticipated for outcome data.
 
ARGUMENTS:  
EXAMPLE: 
library(survival)

data(ex5_widening)
km1 <- with(ex5_widening, survfit(Surv(month, evntd) ~ trt))
km1
plot(km1)
 


FUNCTION: ex6_crossing 
TITLE: Time-to-event data example 6 for non-proportional hazards working group 
DESCRIPTION: 
Survival objects reverse-engineered datasets from published Kaplan-Meier
curves.
Individual trials are de-identified since the data are only
approximations of the actual data.
Data are intended to evaluate methods and designs for trials where
non-proportional hazards may be anticipated for outcome data.
 
ARGUMENTS:  
EXAMPLE: 
library(survival)

data(ex6_crossing)
km1 <- with(ex6_crossing, survfit(Surv(month, evntd) ~ trt))
km1
plot(km1)
 


FUNCTION: fh 
TITLE: Fleming-Harrington weighting function 
DESCRIPTION: 
Fleming-Harrington weighting function
 
ARGUMENTS: rho Non-negative number. \verb{rho = 0, gamma = 0} is equivalent to regular logrank test.
gamma Non-negative number. \verb{rho = 0, gamma = 0} is equivalent to regular logrank test. 
EXAMPLE: 
sim_pw_surv(n = 200) |>
  cut_data_by_event(100) |>
  wlr(weight = fh(rho = 0, gamma = 1))
 


FUNCTION: fit_pwexp 
TITLE: Piecewise exponential survival estimation 
DESCRIPTION: 
Computes survival function, density function, -2 * log-likelihood based
on input dataset and intervals for piecewise constant failure rates.
Initial version assumes observations are right censored or events only.
 
ARGUMENTS: srv Input survival object (see \code{\link[survival:Surv]{survival::Surv()}});
note that only 0 = censored, 1 = event for \code{\link[survival:Surv]{survival::Surv()}}.
intervals Vector containing positive values indicating
interval lengths where the exponential rates are assumed.
Note that a final infinite interval is added if any events occur
after the final interval specified. 
EXAMPLE: 
# Use default arguments for delayed effect example dataset (ex1_delayed_effect)
library(survival)

# Example 1
rateall <- fit_pwexp()
rateall

# Example 2
# Estimate by treatment effect
rate1 <- with(subset(ex1_delayed_effect, trt == 1), fit_pwexp(Surv(month, evntd)))
rate0 <- with(subset(ex1_delayed_effect, trt == 0), fit_pwexp(Surv(month, evntd)))

rate1
rate0
rate1$rate / rate0$rate

# Chi-square test for (any) treatment effect (8 - 4 parameters = 4 df)
pchisq(sum(rateall$m2ll) - sum(rate1$m2ll + rate0$m2ll),
  df = 4,
  lower.tail = FALSE
)

# Compare with logrank
survdiff(formula = Surv(month, evntd) ~ trt, data = ex1_delayed_effect)

# Example 3
# Simple model with 3 rates same for each for 3 months,
# different for each treatment after months
rate1a <- with(subset(ex1_delayed_effect, trt == 1), fit_pwexp(Surv(month, evntd), 3))
rate0a <- with(subset(ex1_delayed_effect, trt == 0), fit_pwexp(Surv(month, evntd), 3))
rate1a$rate / rate0a$rate

m2ll0 <- rateall$m2ll[1] + rate1a$m2ll[2] + rate0a$m2ll[2]
m2ll1 <- sum(rate0$m2ll) + sum(rate1$m2ll)

# As a measure of strength, chi-square examines improvement in likelihood
pchisq(m2ll0 - m2ll1, df = 5, lower.tail = FALSE)
 


FUNCTION: get_analysis_date 
TITLE: Derive analysis date for interim/final analysis given multiple conditions 
DESCRIPTION: 
Derive analysis date for interim/final analysis given multiple conditions
 
ARGUMENTS: data A simulated data generated by \code{\link[=sim_pw_surv]{sim_pw_surv()}}.
planned_calendar_time A numerical value specifying the
planned calendar time for the analysis.
target_event_overall A numerical value specifying the
targeted events for the overall population.
target_event_per_stratum A numerical vector specifying the
targeted events per stratum.
max_extension_for_target_event A numerical value specifying the
maximum time extension to reach targeted events.
previous_analysis_date A numerical value specifying the
previous analysis date.
min_time_after_previous_analysis A numerical value specifying the
planned minimum time after the previous analysis.
min_n_overall A numerical value specifying the
minimal overall sample size enrolled to kick off the analysis.
min_n_per_stratum A numerical value specifying the
minimal sample size enrolled per stratum to kick off the analysis.
min_followup A numerical value specifying the
minimal follow-up time after specified enrollment fraction in
\code{min_n_overall} or \code{min_n_per_stratum}. 
EXAMPLE: 
\dontshow{if (requireNamespace("gsDesign2", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)({ # examplesIf}
library(gsDesign2)

alpha <- 0.025
ratio <- 3
n <- 500
info_frac <- c(0.7, 1)
prevalence_ratio <- c(0.4, 0.6)
study_duration <- 48

# Two strata
stratum <- c("Biomarker-positive", "Biomarker-negative")

prevalence_ratio <- c(0.6, 0.4)
# enrollment rate
enroll_rate <- define_enroll_rate(
  stratum = rep(stratum, each = 2),
  duration = c(2, 10, 2, 10),
  rate = c(c(1, 4) * prevalence_ratio[1], c(1, 4) * prevalence_ratio[2])
)
enroll_rate$rate <- enroll_rate$rate * n / sum(enroll_rate$duration * enroll_rate$rate)

# Failure rate
med_pos <- 10 # Median of the biomarker positive population
med_neg <- 8 # Median of the biomarker negative population
hr_pos <- c(1, 0.7) # Hazard ratio of the biomarker positive population
hr_neg <- c(1, 0.8) # Hazard ratio of the biomarker negative population
fail_rate <- define_fail_rate(
  stratum = rep(stratum, each = 2),
  duration = 1000,
  fail_rate = c(log(2) / c(med_pos, med_pos, med_neg, med_neg)),
  hr = c(hr_pos, hr_neg),
  dropout_rate = 0.01
)

# Simulate data
temp <- to_sim_pw_surv(fail_rate) # Convert the failure rate
set.seed(2023)
simulated_data <- sim_pw_surv(
  n = n, # Sample size
  # Stratified design with prevalence ratio of 6:4
  stratum = data.frame(stratum = stratum, p = prevalence_ratio),
  # Randomization ratio
  block = c("control", "control", "experimental", "experimental"),
  enroll_rate = enroll_rate, # Enrollment rate
  fail_rate = temp$fail_rate, # Failure rate
  dropout_rate = temp$dropout_rate # Dropout rate
)

# Example 1: Cut for analysis at the 24th month.
# Here, we only utilize the `planned_calendar_time = 24` argument,
# while leaving the remaining unused arguments as their default value of `NA`.
get_analysis_date(
  simulated_data,
  planned_calendar_time = 24
)

# Example 2: Cut for analysis when there are 300 events in the overall population.
# Here, we only utilize the `target_event_overall = 300` argument,
# while leaving the remaining unused arguments as their default value of `NA`.
get_analysis_date(
  simulated_data,
  target_event_overall = 300
)

# Example 3: Cut for analysis at the 24th month and there are 300 events
# in the overall population, whichever arrives later.
# Here, we only utilize the `planned_calendar_time = 24` and
# `target_event_overall = 300` argument,
# while leaving the remaining unused arguments as their default value of `NA`.
get_analysis_date(
  simulated_data,
  planned_calendar_time = 24,
  target_event_overall = 300
)

# Example 4a: Cut for analysis when there are at least 100 events
# in the biomarker-positive population, and at least 200 events
# in the biomarker-negative population, whichever arrives later.
# Here, we only utilize the `target_event_per_stratum = c(100, 200)`,
# which refers to 100 events in the biomarker-positive population,
# and 200 events in the biomarker-negative population.
# The remaining unused arguments as their default value of `NA`,
# so the analysis date is only decided by the number of events
# in each stratum.
get_analysis_date(
  simulated_data,
  target_event_per_stratum = c(100, 200)
)
# Example 4b: Cut for analysis when there are at least 100 events
# in the biomarker-positive population, but we don't have a requirement
# for the biomarker-negative population. Additionally, we want to cut
# the analysis when there are at least 150 events in total.
# Here, we only utilize the `target_event_overall = 150` and
# `target_event_per_stratum = c(100, NA)`, which refers to 100 events
# in the biomarker-positive population, and there is event requirement
# for the biomarker-negative population.
# The remaining unused arguments as their default value of `NA`,
# so the analysis date is only decided by the number of events
# in the biomarker-positive population, and the total number of events,
# which arrives later.
get_analysis_date(
  simulated_data,
  target_event_overall = 150,
  target_event_per_stratum = c(100, NA)
)
# Example 4c: Cut for analysis when there are at least 100 events
# in the biomarker-positive population, but we don't have a requirement
# for the biomarker-negative population. Additionally, we want to cut
# the analysis when there are at least 150 events in total and after 24 months.
# Here, we only utilize the `planned_calendar_time = 24`,
# `target_event_overall = 150` and
# `target_event_per_stratum = c(100, NA)`, which refers to 100 events
# in the biomarker-positive population, and there is event requirement
# for the biomarker-negative population.
# The remaining unused arguments as their default value of `NA`,
# so the analysis date is only decided by the number of events
# in the biomarker-positive population, the total number of events, and
# planned calendar time, which arrives later.
get_analysis_date(
  simulated_data,
  planned_calendar_time = 24,
  target_event_overall = 150,
  target_event_per_stratum = c(100, NA)
)

# Example 5: Cut for analysis when there are at least 100 events
# in the biomarker positive population, and at least 200 events
# in the biomarker negative population, whichever arrives later.
# But will stop at the 30th month if events are fewer than 100/200.
# Here, we only utilize the `max_extension_for_target_event = 30`,
# and `target_event_per_stratum =  c(100, 200)`, which refers to
# 100/200 events in the biomarker-positive/negative population.
# The remaining unused arguments as their default value of `NA`,
# so the analysis date is only decided by the number of events
# in the 2 strata, and the max extension to arrive at the targeted
# events, which arrives later.
get_analysis_date(
  simulated_data,
  target_event_per_stratum = c(100, 200),
  max_extension_for_target_event = 30
)

# Example 6a: Cut for analysis after 12 months followup when 80%
# of the patients are enrolled in the overall population.
# The remaining unused arguments as their default value of `NA`,
# so the analysis date is only decided by
# 12 months + time when 80% patients enrolled.
get_analysis_date(
  simulated_data,
  min_n_overall = n * 0.8,
  min_followup = 12
)
# Example 6b: Cut for analysis after 12 months followup when 80%
# of the patients are enrolled in the overall population. Besides,
# the analysis happens when there are at least 150 events in total.
# The remaining unused arguments as their default value of `NA`,
# so the analysis date is only decided by the total number of events,
# and 12 months + time when 80% patients enrolled, which arrives later.
get_analysis_date(
  simulated_data,
  target_event_overall = 150,
  min_n_overall = n * 0.8,
  min_followup = 12
)

# Example 7a: Cut for analysis when 12 months after at least 200/160 patients
# are enrolled in the biomarker positive/negative population.
# The remaining unused arguments as their default value of `NA`,
# so the analysis date is only decided by 12 months + time when there are
# 200/160 patients enrolled in the biomarker-positive/negative stratum.
get_analysis_date(
  simulated_data,
  min_n_per_stratum = c(200, 160),
  min_followup = 12
)
# Example 7b: Cut for analysis when 12 months after at least 200 patients
# are enrolled in the biomarker positive population, but we don't have a
# specific requirement for the biomarker negative population.
# The remaining unused arguments as their default value of `NA`,
# so the analysis date is only decided by 12 months + time when there are
# 200 patients enrolled in the biomarker-positive stratum.
get_analysis_date(
  simulated_data,
  min_n_per_stratum = c(200, NA),
  min_followup = 12
)
# Example 7c: Cut for analysis when 12 months after at least 200 patients
# are enrolled in the biomarker-positive population, but we don't have a
# specific requirement for the biomarker-negative population. We also want
# there are at least 80% of the patients enrolled in the overall population.
# The remaining unused arguments as their default value of `NA`,
# so the analysis date is only decided by 12 months + max(time when there are
# 200 patients enrolled in the biomarker-positive stratum, time when there are
# 80% patients enrolled).
get_analysis_date(
  simulated_data,
  min_n_overall = n * 0.8,
  min_n_per_stratum = c(200, NA),
  min_followup = 12
)
\dontshow{}) # examplesIf}
 


FUNCTION: get_cut_date_by_event 
TITLE: Get date at which an event count is reached 
DESCRIPTION: 
Get date at which an event count is reached
 
ARGUMENTS: x A time-to-event dataset, for example, generated by \code{\link[=sim_pw_surv]{sim_pw_surv()}}.
event Event count at which dataset is to be cut off for analysis. 
EXAMPLE: 
\dontshow{if (requireNamespace("dplyr", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)({ # examplesIf}
library(dplyr)

# Use default enrollment and calendar cut date
# for 50 events in the "Positive" stratum
x <- sim_pw_surv(
  n = 200,
  stratum = data.frame(
    stratum = c("Positive", "Negative"),
    p = c(.5, .5)
  ),
  fail_rate = data.frame(
    stratum = rep(c("Positive", "Negative"), 2),
    period = rep(1, 4),
    treatment = c(rep("control", 2), rep("experimental", 2)),
    duration = rep(1, 4),
    rate = log(2) / c(6, 9, 9, 12)
  ),
  dropout_rate = data.frame(
    stratum = rep(c("Positive", "Negative"), 2),
    period = rep(1, 4),
    treatment = c(rep("control", 2), rep("experimental", 2)),
    duration = rep(1, 4),
    rate = rep(.001, 4)
  )
)

d <- get_cut_date_by_event(x |> filter(stratum == "Positive"), event = 50)

y <- cut_data_by_date(x, cut_date = d)
table(y$stratum, y$event)
\dontshow{}) # examplesIf}
 


FUNCTION: maxcombo 
TITLE: MaxCombo test 
DESCRIPTION: 
WARNING: This experimental function is a work-in-progress. The function
arguments will change as we add additional features.
 
ARGUMENTS: data A TTE dataset.
rho Numeric vector. Must be greater
than or equal to zero. Must be the same length as \code{gamma}.
gamma Numeric vector. Must be
greater than or equal to zero. Must be the same length as \code{rho}.
return_variance A logical flag that, if \code{TRUE}, adds columns
estimated variance for weighted sum of observed minus expected;
see details; Default: \code{FALSE}.
return_corr A logical flag that, if \code{TRUE}, adds columns
estimated correlation for weighted sum of observed minus expected;
see details; Default: \code{FALSE}. 
EXAMPLE: 
sim_pw_surv(n = 200) |>
  cut_data_by_event(150) |>
  maxcombo(rho = c(0, 0), gamma = c(0, 1), return_corr = TRUE)
 


FUNCTION: mb_delayed_effect 
TITLE: Simulated survival dataset with delayed treatment effect 
DESCRIPTION: 
Magirr and Burman (2019) considered several scenarios for their
modestly weighted logrank test.
One of these had a delayed treatment effect with a hazard ratio
of 1 for 6 months followed by a hazard ratio of 1/2 thereafter.
The scenario enrolled 200 patients uniformly over 12 months and
cut data for analysis 36 months after enrollment was opened.
This dataset was generated by the \code{\link[=sim_pw_surv]{sim_pw_surv()}} function
under the above scenario.
 
ARGUMENTS:  
EXAMPLE: 
library(survival)

fit <- survfit(Surv(tte, event) ~ treatment, data = mb_delayed_effect)

# Plot survival
plot(fit, lty = 1:2)
legend("topright", legend = c("control", "experimental"), lty = 1:2)

# Set up time, event, number of event dataset for testing
# with arbitrary weights
ten <- mb_delayed_effect |> counting_process(arm = "experimental")
head(ten)

# MaxCombo with logrank, FH(0,1), FH(1,1)
mb_delayed_effect |>
  maxcombo(rho = c(0, 0, 1), gamma = c(0, 1, 1), return_corr = TRUE)

# Generate another dataset
ds <- sim_pw_surv(
  n = 200,
  enroll_rate = data.frame(rate = 200 / 12, duration = 12),
  fail_rate = data.frame(
    stratum = c("All", "All", "All"),
    period = c(1, 1, 2),
    treatment = c("control", "experimental", "experimental"),
    duration = c(42, 6, 36),
    rate = c(log(2) / 15, log(2) / 15, log(2) / 15 * 0.6)
  ),
  dropout_rate = data.frame(
    stratum = c("All", "All"),
    period = c(1, 1),
    treatment = c("control", "experimental"),
    duration = c(42, 42),
    rate = c(0, 0)
  )
)
# Cut data at 24 months after final enrollment
mb_delayed_effect_2 <- ds |> cut_data_by_date(max(ds$enroll_time) + 24)
 


FUNCTION: mb 
TITLE: Magirr and Burman weighting function 
DESCRIPTION: 
Magirr and Burman weighting function
 
ARGUMENTS: delay The initial delay period where weights increase;
after this, weights are constant at the final weight in the delay period.
w_max Maximum weight to be returned.
Set \code{delay = Inf}, \code{w_max = 2} to be consistent with recommendation of
Magirr (2021). 
EXAMPLE: 
sim_pw_surv(n = 200) |>
  cut_data_by_event(100) |>
  wlr(weight = mb(delay = 8, w_max = Inf))
 


FUNCTION: milestone 
TITLE: Milestone test for two survival curves 
DESCRIPTION: 
Milestone test for two survival curves
 
ARGUMENTS: data Data frame containing at least 3 columns:
\itemize{
\item \code{tte} - Time to event.
\item \code{event} - Event indicator.
\item \code{treatment} - Grouping variable.
}
ms_time Milestone analysis time.
test_type Method to build the test statistics.
There are 2 options:
\itemize{
\item \code{"native"}: a native approach by dividing the KM survival difference by its standard derivations,
see equation (1) of Klein, J. P., Logan, B., Harhoff, M., & Andersen, P. K. (2007).
\item \code{"log-log"}: a log-log transformation of the survival, see equation (3) of
Klein, J. P., Logan, B., Harhoff, M., & Andersen, P. K. (2007).
} 
EXAMPLE: 
cut_data <- sim_pw_surv(n = 200) |>
  cut_data_by_event(150)

cut_data |>
  milestone(10, test_type = "log-log")

cut_data |>
  milestone(10, test_type = "naive")
 


FUNCTION: multitest 
TITLE: Perform multiple tests on trial data cutting 
DESCRIPTION: 
WARNING: This experimental function is a work-in-progress. The function
arguments and/or returned output format may change as we add additional
features.
 
ARGUMENTS: data Trial data cut by \code{\link[=cut_data_by_event]{cut_data_by_event()}} or \code{\link[=cut_data_by_date]{cut_data_by_date()}}
... One or more test functions. Use \code{\link[=create_test]{create_test()}} to change
the default arguments of each test function. 
EXAMPLE: 
trial_data <- sim_pw_surv(n = 200)
trial_data_cut <- cut_data_by_event(trial_data, 150)

# create cutting test functions
wlr_partial <- create_test(wlr, weight = fh(rho = 0, gamma = 0))
rmst_partial <- create_test(rmst, tau = 20)
maxcombo_partial <- create_test(maxcombo, rho = c(0, 0), gamma = c(0, 0.5))

multitest(
  data = trial_data_cut,
  wlr = wlr_partial,
  rmst = rmst_partial,
  maxcombo = maxcombo_partial
)
 


FUNCTION: randomize_by_fixed_block 
TITLE: Permuted fixed block randomization 
DESCRIPTION: 
Fixed block randomization. The \code{block} input should repeat each
treatment code the number of times it is to be included within each block.
The final block will be a partial block if \code{n} is not an exact multiple
of the block length.
 
ARGUMENTS: n Sample size to be randomized.
block Vector of treatments to be included in each block. 
EXAMPLE: 
\dontshow{if (requireNamespace("dplyr", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)({ # examplesIf}
library(dplyr)

# Example 1
# 2:1 randomization with block size 3, treatments "A" and "B"
data.frame(x = 1:10) |> mutate(Treatment = randomize_by_fixed_block(block = c("A", "B", "B")))

# Example 2
# Stratified randomization
data.frame(stratum = c(rep("A", 10), rep("B", 10))) |>
  group_by(stratum) |>
  mutate(Treatment = randomize_by_fixed_block())
\dontshow{}) # examplesIf}
 


FUNCTION: rmst_single_arm 
TITLE: Calculate RMST for a single cut-off time point 
DESCRIPTION: 
Calculate RMST for a single cut-off time point
 
ARGUMENTS: time_var A numeric vector of follow up time.
event_var A numeric or integer vector of the status indicator;
0=alive 1=event.
tau A value of pre-defined cut-off time point.
group_label A character of customized treatment group name.
alpha A numeric value of the significant level for RMST
confidence interval. Default is 0.05. 
EXAMPLE: 
data(ex1_delayed_effect)
data_single_arm <- ex1_delayed_effect[ex1_delayed_effect$trt == 1, ]
simtrial:::rmst_single_arm(
  time_var = data_single_arm$month,
  event_var = data_single_arm$evntd,
  tau = 10,
  group_label = "Treatment 1",
  alpha = 0.05
)
 


FUNCTION: rmst_two_arm 
TITLE: Calculate RMST difference 
DESCRIPTION: 
Calculate RMST difference
 
ARGUMENTS: time_var A numeric vector of follow up time.
event_var A numeric or integer vector of the status indicator;
0=alive 1=event.
group_var A vector of treatment groups.
trunc_time A numeric vector of pre-defined cut-off time point(s).
reference Group name of reference group for RMST comparison.
Default is the first group name by alphabetical order.
alpha A numeric value of the significant level for RMST
confidence interval. Default is 0.05. 
EXAMPLE: 
data(ex1_delayed_effect)
with(
  ex1_delayed_effect,
  simtrial:::rmst_two_arm(
    time_var = month,
    event_var = evntd,
    group_var = trt,
    trunc_time = 6,
    reference = "0",
    alpha = 0.05
  )
)
 


FUNCTION: rmst 
TITLE: RMST difference of 2 arms 
DESCRIPTION: 
RMST difference of 2 arms
 
ARGUMENTS: data A time-to-event dataset with a column \code{tte} indicating the
survival time and a column of \code{event} indicating whether it is
event or censor.
tau RMST analysis time.
var_label_tte Column name of the TTE variable.
var_label_event Column name of the event variable.
var_label_group Column name of the grouping variable.
formula (default: \code{NULL}) A formula that indicates the TTE, event, and
group variables using the syntax \verb{Surv(tte, event) ~ group)} (see Details
below for more information). This is an alternative to specifying the
variables as strings. If a formula is provided, the values passed to
\code{var_label_tte}, \code{var_label_event}, and \code{var_label_group} are ignored.
reference A group label indicating the reference group.
alpha Type I error. 
EXAMPLE: 
data(ex1_delayed_effect)
rmst(
  data = ex1_delayed_effect,
  var_label_tte = "month",
  var_label_event = "evntd",
  var_label_group = "trt",
  tau = 10,
  reference = "0"
)

# Formula interface
rmst(
  data = ex1_delayed_effect,
  formula = Surv(month, evntd) ~ trt,
  tau = 10,
  reference = "0"
)
 


FUNCTION: rpwexp_enroll 
TITLE: Generate piecewise exponential enrollment 
DESCRIPTION: 
With piecewise exponential enrollment rate generation any enrollment rate
distribution can be easily approximated.
\code{rpwexp_enroll()} is to support simulation of both the Lachin and Foulkes (1986)
sample size method for (fixed trial duration) as well as the
Kim and Tsiatis(1990) method (fixed enrollment rates and either
fixed enrollment duration or fixed minimum follow-up);
see \code{\link[gsDesign:nSurv]{gsDesign::nSurv()}}.
 
ARGUMENTS: n Number of observations.
Default of \code{NULL} yields random enrollment size.
enroll_rate A data frame containing period duration (\code{duration})
and enrollment rate (\code{rate}). for specified enrollment periods.
If necessary, last period will be extended to ensure enrollment
of specified \code{n}. 
EXAMPLE: 
# Example 1
# Piecewise uniform (piecewise exponential inter-arrival times) for 10k patients enrollment
# Enrollment rates of 5 for time 0-100, 15 for 100-300, and 30 thereafter
x <- rpwexp_enroll(
  n = 1e5,
  enroll_rate = data.frame(
    rate = c(5, 15, 30),
    duration = c(100, 200, 100)
  )
)
plot(x, 1:1e5,
  main = "Piecewise uniform enrollment simulation",
  xlab = "Time",
  ylab = "Enrollment"
)

# Example 2
# Exponential enrollment
x <- rpwexp_enroll(
  n = 1e5,
  enroll_rate = data.frame(rate = .03, duration = 1)
)
plot(x, 1:1e5,
  main = "Simulated exponential inter-arrival times",
  xlab = "Time",
  ylab = "Enrollment"
)
 


FUNCTION: rpwexp 
TITLE: The piecewise exponential distribution 
DESCRIPTION: 
The piecewise exponential distribution allows a simple method to specify
a distribution where the hazard rate changes over time.
It is likely to be useful for conditions where failure rates change,
but also for simulations where there may be a delayed treatment effect
or a treatment effect that that is otherwise changing
(for example, decreasing) over time.
\code{rpwexp()} is to support simulation of both the Lachin and Foulkes (1986)
sample size method for (fixed trial duration) as well as the
Kim and Tsiatis (1990) method (fixed enrollment rates and either
fixed enrollment duration or fixed minimum follow-up);
see \code{\link[gsDesign:nSurv]{gsDesign::nSurv()}}.
 
ARGUMENTS: n Number of observations to be generated.
fail_rate A data frame containing \code{duration} and \code{rate} variables.
\code{rate} specifies failure rates during the corresponding interval duration
specified in \code{duration}. The final interval is extended to be infinite
to ensure all observations are generated. 
EXAMPLE: 
# Example 1
# Exponential failure times
x <- rpwexp(
  n = 10000,
  fail_rate = data.frame(rate = 5, duration = 1)
)
plot(sort(x), (10000:1) / 10001,
  log = "y", main = "Exponential simulated survival curve",
  xlab = "Time", ylab = "P{Survival}"
)

# Example 2

# Get 10k piecewise exponential failure times.
# Failure rates are 1 for time 0 to 0.5, 3 for time 0.5 to 1, and 10 for > 1.
# Intervals specifies duration of each failure rate interval
# with the final interval running to infinity.
x <- rpwexp(
  n = 1e4,
  fail_rate = data.frame(rate = c(1, 3, 10), duration = c(.5, .5, 1))
)
plot(sort(x), (1e4:1) / 10001,
  log = "y", main = "PW Exponential simulated survival curve",
  xlab = "Time", ylab = "P{Survival}"
)
 


FUNCTION: sim_fixed_n 
TITLE: Simulation of fixed sample size design for time-to-event endpoint 
DESCRIPTION: 
\code{sim_fixed_n()} provides simulations of a single endpoint two-arm trial
where the enrollment, hazard ratio, and failure and dropout rates change
over time.
 
ARGUMENTS: n_sim Number of simulations to perform.
sample_size Total sample size per simulation.
target_event Targeted event count for analysis.
stratum A data frame with stratum specified in \code{stratum},
probability (incidence) of each stratum in \code{p}.
enroll_rate Piecewise constant enrollment rates by time period.
Note that these are overall population enrollment rates and the \code{stratum}
argument controls the random distribution between stratum.
fail_rate Piecewise constant control group failure rates, hazard ratio
for experimental vs. control, and dropout rates by stratum and time period.
total_duration Total follow-up from start of enrollment to data cutoff.
block As in \code{\link[=sim_pw_surv]{sim_pw_surv()}}. Vector of treatments to be included
in each block.
timing_type A numeric vector determining data cutoffs used;
see details. Default is to include all available cutoff methods.
rho_gamma A data frame with variables
\code{rho} and \code{gamma}, both greater than equal to zero,
to specify one Fleming-Harrington weighted logrank test per row. 
EXAMPLE: 
\dontshow{if (requireNamespace("dplyr", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)({ # examplesIf}
library(dplyr)
library(future)

# Example 1: logrank test ----
x <- sim_fixed_n(n_sim = 10, timing_type = 1, rho_gamma = data.frame(rho = 0, gamma = 0))
# Get power approximation
mean(x$z <= qnorm(.025))

# Example 2: WLR with FH(0,1) ----
sim_fixed_n(n_sim = 1, timing_type = 1, rho_gamma = data.frame(rho = 0, gamma = 1))
# Get power approximation
mean(x$z <= qnorm(.025))

\donttest{
# Example 3: MaxCombo, i.e., WLR-FH(0,0)+ WLR-FH(0,1)
# Power by test
# Only use cuts for events, events + min follow-up
x <- sim_fixed_n(
  n_sim = 10,
  timing_type = 2,
  rho_gamma = data.frame(rho = 0, gamma = c(0, 1))
)

# Get power approximation
x |>
  group_by(sim) |>
  filter(row_number() == 1) |>
  ungroup() |>
  summarize(power = mean(p_value < .025))

# Example 4
# Use two cores
set.seed(2023)
plan("multisession", workers = 2)
sim_fixed_n(n_sim = 10)
plan("sequential")
}
\dontshow{}) # examplesIf}
 


FUNCTION: sim_gs_n 
TITLE: Simulate group sequential designs with fixed sample size 
DESCRIPTION: 
This function uses the option "stop" for the error-handling behavior of the
foreach loop. This will cause the entire function to stop when errors are
encountered and return the first error encountered instead of returning
errors for each individual simulation.
 
ARGUMENTS: n_sim Number of simulations to perform.
sample_size Total sample size per simulation.
stratum A data frame with stratum specified in \code{stratum},
probability (incidence) of each stratum in \code{p}.
enroll_rate Piecewise constant enrollment rates by time period.
Note that these are overall population enrollment rates and the \code{stratum}
argument controls the random distribution between stratum.
fail_rate Piecewise constant control group failure rates, hazard ratio
for experimental vs. control, and dropout rates by stratum and time period.
block As in \code{\link[=sim_pw_surv]{sim_pw_surv()}}. Vector of treatments to be included
in each block.
test One or more test functions such as \code{\link[=wlr]{wlr()}}, \code{\link[=rmst]{rmst()}}, or
\code{\link[=milestone]{milestone()}} (\code{\link[=maxcombo]{maxcombo()}} can only be applied by itself). If a single
test function is provided, it will be applied at each cut. Alternatively a
list of functions created by \code{\link[=create_test]{create_test()}}. The list form is experimental
and currently limited. It only accepts one test per cutting (in the future
multiple tests may be accepted), and all the tests must consistently return
the same exact results (again this may be more flexible in the future).
Importantly, note that the simulated data set is always passed as the first
positional argument to each test function provided.
cut A list of cutting functions created by \code{\link[=create_cut]{create_cut()}}, see
examples.
original_design A design object from the gsDesign2 package, which is required when users
want to calculate updated bounds. The default is NULL leaving the updated bounds uncalculated.
ia_alpha_spending Spend alpha at interim analysis based on
\itemize{
\item \code{"min_planned_actual"}: the minimal of planned and actual alpha spending.
\item \code{"actual"}: the actual alpha spending.
}
fa_alpha_spending If targeted final event count is not achieved (under-running at final analysis),
specify how to do final spending. Generally, this should be specified in analysis plan.
\itemize{
\item \code{"info_frac"} = spend final alpha according to final information fraction
\item \code{"full_alpha"} = spend full alpha at final analysis.
}
... Arguments passed to the test function(s) provided by the argument
\code{test}. 
EXAMPLE: 
\dontshow{if (requireNamespace("gsDesign2", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)({ # examplesIf}
library(gsDesign2)

# Parameters for enrollment
enroll_rampup_duration <- 4 # Duration for enrollment ramp up
enroll_duration <- 16 # Total enrollment duration
enroll_rate <- define_enroll_rate(
  duration = c(
    enroll_rampup_duration,
    enroll_duration - enroll_rampup_duration
  ),
  rate = c(10, 30)
)

# Parameters for treatment effect
delay_effect_duration <- 3 # Delay treatment effect in months
median_ctrl <- 9 # Survival median of the control arm
median_exp <- c(9, 14) # Survival median of the experimental arm
dropout_rate <- 0.001
fail_rate <- define_fail_rate(
  duration = c(delay_effect_duration, 100),
  fail_rate = log(2) / median_ctrl,
  hr = median_ctrl / median_exp,
  dropout_rate = dropout_rate
)

# Other related parameters
alpha <- 0.025 # Type I error
beta <- 0.1 # Type II error
ratio <- 1 # Randomization ratio (experimental:control)

# Define cuttings of 2 IAs and 1 FA
# IA1
# The 1st interim analysis will occur at the later of the following 3 conditions:
# - At least 20 months have passed since the start of the study.
# - At least 100 events have occurred.
# - At least 20 months have elapsed after enrolling 200/400 subjects, with a
#   minimum of 20 months follow-up.
# However, if events accumulation is slow, we will wait for a maximum of 24 months.
ia1_cut <- create_cut(
  planned_calendar_time = 20,
  target_event_overall = 100,
  max_extension_for_target_event = 24,
  min_n_overall = 200,
  min_followup = 20
)

# IA2
# The 2nd interim analysis will occur at the later of the following 3 conditions:
# - At least 32 months have passed since the start of the study.
# - At least 250 events have occurred.
# - At least 10 months after IA1.
# However, if events accumulation is slow, we will wait for a maximum of 34 months.
ia2_cut <- create_cut(
  planned_calendar_time = 32,
  target_event_overall = 200,
  max_extension_for_target_event = 34,
  min_time_after_previous_analysis = 10
)

# FA
# The final analysis will occur at the later of the following 2 conditions:
# - At least 45 months have passed since the start of the study.
# - At least 300 events have occurred.
fa_cut <- create_cut(
  planned_calendar_time = 45,
  target_event_overall = 350
)

# Example 1: regular logrank test at all 3 analyses
sim_gs_n(
  n_sim = 3,
  sample_size = 400,
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  test = wlr,
  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
  weight = fh(rho = 0, gamma = 0)
)

\donttest{
# Example 2: weighted logrank test by FH(0, 0.5) at all 3 analyses
sim_gs_n(
  n_sim = 3,
  sample_size = 400,
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  test = wlr,
  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
  weight = fh(rho = 0, gamma = 0.5)
)

# Example 3: weighted logrank test by MB(3) at all 3 analyses
sim_gs_n(
  n_sim = 3,
  sample_size = 400,
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  test = wlr,
  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
  weight = mb(delay = 3)
)

# Example 4: weighted logrank test by early zero (6) at all 3 analyses
sim_gs_n(
  n_sim = 3,
  sample_size = 400,
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  test = wlr,
  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
  weight = early_zero(6)
)

# Example 5: RMST at all 3 analyses
sim_gs_n(
  n_sim = 3,
  sample_size = 400,
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  test = rmst,
  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
  tau = 20
)

# Example 6: Milestone at all 3 analyses
sim_gs_n(
  n_sim = 3,
  sample_size = 400,
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  test = milestone,
  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
  ms_time = 10
)
}

# Warning: this example will be executable when we add info info0 to the milestone test
# Example 7: WLR with fh(0, 0.5) test at IA1,
# WLR with mb(6, Inf) at IA2, and milestone test at FA
ia1_test <- create_test(wlr, weight = fh(rho = 0, gamma = 0.5))
ia2_test <- create_test(wlr, weight = mb(delay = 6, w_max = Inf))
fa_test <- create_test(milestone, ms_time = 10)
\dontrun{
sim_gs_n(
  n_sim = 3,
  sample_size = 400,
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  test = list(ia1 = ia1_test, ia2 = ia2_test, fa = fa_test),
  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut)
)
}

# WARNING: Multiple tests per cut will be enabled in a future version.
#          Currently does not work.
# Example 8: At IA1, we conduct 3 tests, LR, WLR with fh(0, 0.5), and RMST test.
# At IA2, we conduct 2 tests, LR and WLR with early zero (6).
# At FA, we conduct 2 tests, LR and milestone test.
ia1_test <- list(
  test1 = create_test(wlr, weight = fh(rho = 0, gamma = 0)),
  test2 = create_test(wlr, weight = fh(rho = 0, gamma = 0.5)),
  test3 = create_test(rmst, tau = 20)
)
ia2_test <- list(
  test1 = create_test(wlr, weight = fh(rho = 0, gamma = 0)),
  test2 = create_test(wlr, weight = early_zero(6))
)
fa_test <- list(
  test1 = create_test(wlr, weight = fh(rho = 0, gamma = 0)),
  test3 = create_test(milestone, ms_time = 20)
)
\dontrun{
sim_gs_n(
  n_sim = 3,
  sample_size = 400,
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  test = list(ia1 = ia1_test, ia2 = ia2_test, fa = fa_test),
  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut)
)
}

\donttest{
# Example 9: regular logrank test at all 3 analyses in parallel
plan("multisession", workers = 2)
sim_gs_n(
  n_sim = 3,
  sample_size = 400,
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  test = wlr,
  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
  weight = fh(rho = 0, gamma = 0)
)
plan("sequential")

# Example 10: group sequential design with updated bounds -- efficacy only
x <- gs_design_ahr(analysis_time = 1:3*12) |> to_integer()
sim_gs_n(
  n_sim = 1,
  sample_size = max(x$analysis$n),
  enroll_rate = x$enroll_rate,
  fail_rate = x$fail_rate,
  test = wlr,
  cut = list(ia1 = create_cut(planned_calendar_time = x$analysis$time[1]),
             ia2 = create_cut(planned_calendar_time = x$analysis$time[2]),
             fa = create_cut(planned_calendar_time = x$analysis$time[3])),
  weight = fh(rho = 0, gamma = 0),
  original_design = x
)

# Example 11: group sequential design with updated bounds -- efficacy & futility
x <- gs_design_ahr(
 alpha = 0.025, beta = 0.1, analysis_time = 1:3*12,
 upper = gs_spending_bound, upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
 lower = gs_spending_bound, lpar = list(sf = gsDesign::sfHSD, param = -4, total_spend = 0.01),
 test_upper = c(FALSE, TRUE, TRUE), test_lower = c(TRUE, FALSE, FALSE)) |> to_integer()
sim_gs_n(
  n_sim = 1,
  sample_size = max(x$analysis$n),
  enroll_rate = x$enroll_rate,
  fail_rate = x$fail_rate,
  test = wlr,
  cut = list(ia1 = create_cut(planned_calendar_time = x$analysis$time[1]),
             ia2 = create_cut(planned_calendar_time = x$analysis$time[2]),
             fa = create_cut(planned_calendar_time = x$analysis$time[3])),
  weight = fh(rho = 0, gamma = 0),
  original_design = x
)
}
\dontshow{}) # examplesIf}
 


FUNCTION: sim_pw_surv 
TITLE: Simulate a stratified time-to-event outcome randomized trial 
DESCRIPTION: 
\code{sim_pw_surv()} enables simulation of a clinical trial with
essentially arbitrary patterns of enrollment, failure rates and censoring.
The piecewise exponential distribution allows a simple method to specify
a distribution and enrollment pattern where the enrollment, failure,
and dropout rate changes over time.
While the main purpose may be to generate a trial that can be analyzed
at a single point in time or using group sequential methods,
the routine can also be used to simulate an adaptive trial design.
Enrollment, failure, and dropout rates are specified by treatment group,
stratum and time period.
Fixed block randomization is used; blocks must include treatments provided
in failure and dropout specification.
Default arguments are set up to allow very simple implementation of
a non-proportional hazards assumption for an unstratified design.
 
ARGUMENTS: n Number of observations.
If length(n) > 1, the length is taken to be the number required.
stratum A data frame with stratum specified in \code{stratum},
probability (incidence) of each stratum in \code{p}.
block Vector of treatments to be included in each block. Also used to
calculate the attribute "ratio" (for more details see the section Value
below).
enroll_rate Enrollment rates; see details and examples.
fail_rate Failure rates; see details and examples;
note that treatments need to be the same as input in block.
dropout_rate Dropout rates; see details and examples;
note that treatments need to be the same as input in block. 
EXAMPLE: 
\dontshow{if (requireNamespace("dplyr", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)({ # examplesIf}
library(dplyr)

# Example 1
sim_pw_surv(n = 20)

# Example 2
# 3:1 randomization
sim_pw_surv(
  n = 20,
  block = c(rep("experimental", 3), "control")
)

# Example 3
# Simulate 2 stratum; will use defaults for blocking and enrollRates
sim_pw_surv(
  n = 20,
  # 2 stratum,30% and 70% prevalence
  stratum = data.frame(stratum = c("Low", "High"), p = c(.3, .7)),
  fail_rate = data.frame(
    stratum = c(rep("Low", 4), rep("High", 4)),
    period = rep(1:2, 4),
    treatment = rep(c(
      rep("control", 2),
      rep("experimental", 2)
    ), 2),
    duration = rep(c(3, 1), 4),
    rate = c(.03, .05, .03, .03, .05, .08, .07, .04)
  ),
  dropout_rate = data.frame(
    stratum = c(rep("Low", 2), rep("High", 2)),
    period = rep(1, 4),
    treatment = rep(c("control", "experimental"), 2),
    duration = rep(1, 4),
    rate = rep(.001, 4)
  )
)
# Example 4
# If you want a more rectangular entry for a data.frame
fail_rate <- bind_rows(
  data.frame(stratum = "Low", period = 1, treatment = "control", duration = 3, rate = .03),
  data.frame(stratum = "Low", period = 1, treatment = "experimental", duration = 3, rate = .03),
  data.frame(stratum = "Low", period = 2, treatment = "experimental", duration = 3, rate = .02),
  data.frame(stratum = "High", period = 1, treatment = "control", duration = 3, rate = .05),
  data.frame(stratum = "High", period = 1, treatment = "experimental", duration = 3, rate = .06),
  data.frame(stratum = "High", period = 2, treatment = "experimental", duration = 3, rate = .03)
)

dropout_rate <- bind_rows(
  data.frame(stratum = "Low", period = 1, treatment = "control", duration = 3, rate = .001),
  data.frame(stratum = "Low", period = 1, treatment = "experimental", duration = 3, rate = .001),
  data.frame(stratum = "High", period = 1, treatment = "control", duration = 3, rate = .001),
  data.frame(stratum = "High", period = 1, treatment = "experimental", duration = 3, rate = .001)
)

sim_pw_surv(
  n = 12,
  stratum = data.frame(stratum = c("Low", "High"), p = c(.3, .7)),
  fail_rate = fail_rate,
  dropout_rate = dropout_rate
)
\dontshow{}) # examplesIf}
 


FUNCTION: simtrial-package 
TITLE: simtrial: Clinical Trial Simulation 
DESCRIPTION: 
\if{html\figure{logo.pngoptions: style='float: right' alt='logo' width='120'}}

Provides some basic routines for simulating a clinical trial. The primary intent is to provide some tools to generate trial simulations for trials with time to event outcomes. Piecewise exponential failure rates and piecewise constant enrollment rates are the underlying mechanism used to simulate a broad range of scenarios such as those presented in Lin et al. (2020) c("\\Sexpr[results=rd]{tools:::Rd_expr_doi(\"#1\")}", "10.1080/19466315.2019.1697738")\Sexpr{tools:::Rd_expr_doi("10.1080/19466315.2019.1697738")}. However, the basic generation of data is done using pipes to allow maximum flexibility for users to meet different needs.
 
ARGUMENTS:  
EXAMPLE:  


FUNCTION: summary 
TITLE: Summary of group sequential simulations. 
DESCRIPTION: 
Summary of group sequential simulations.
 
ARGUMENTS: object Simulation results generated by \code{\link[=sim_gs_n]{sim_gs_n()}}
design Asymptotic design generated by \code{\link[gsDesign2:gs_design_ahr]{gsDesign2::gs_design_ahr()}},
\code{\link[gsDesign2:gs_power_ahr]{gsDesign2::gs_power_ahr()}}, \code{\link[gsDesign2:gs_design_wlr]{gsDesign2::gs_design_wlr()}}, or \link[gsDesign2:gs_power_wlr]{gsDesign2::gs_power_wlr}.
bound The boundaries.
... Additional parameters (not used). 
EXAMPLE: 
library(gsDesign2)

# Parameters for enrollment
enroll_rampup_duration <- 4 # Duration for enrollment ramp up
enroll_duration <- 16 # Total enrollment duration
enroll_rate <- define_enroll_rate(
  duration = c(
    enroll_rampup_duration, enroll_duration - enroll_rampup_duration),
 rate = c(10, 30))

# Parameters for treatment effect
delay_effect_duration <- 3 # Delay treatment effect in months
median_ctrl <- 9 # Survival median of the control arm
median_exp <- c(9, 14) # Survival median of the experimental arm
dropout_rate <- 0.001
fail_rate <- define_fail_rate(
  duration = c(delay_effect_duration, 100),
  fail_rate = log(2) / median_ctrl,
  hr = median_ctrl / median_exp,
  dropout_rate = dropout_rate)

# Other related parameters
alpha <- 0.025 # Type I error
beta <- 0.1 # Type II error
ratio <- 1 # Randomization ratio (experimental:control)

# Build a one-sided group sequential design
design <- gs_design_ahr(
  enroll_rate = enroll_rate, fail_rate = fail_rate,
  ratio = ratio, alpha = alpha, beta = beta,
  analysis_time = c(12, 24, 36),
  upper = gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = alpha),
  lower = gs_b,
  lpar = rep(-Inf, 3))

# Define cuttings of 2 IAs and 1 FA
ia1_cut <- create_cut(target_event_overall = ceiling(design$analysis$event[1]))
ia2_cut <- create_cut(target_event_overall = ceiling(design$analysis$event[2]))
fa_cut <- create_cut(target_event_overall = ceiling(design$analysis$event[3]))

# Run simulations
simulation <- sim_gs_n(
  n_sim = 3,
  sample_size = ceiling(design$analysis$n[3]),
  enroll_rate = design$enroll_rate,
  fail_rate = design$fail_rate,
  test = wlr,
  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
  weight = fh(rho = 0, gamma = 0.5))

# Summarize simulations
bound <- gsDesign::gsDesign(k = 3, test.type = 1, sfu = gsDesign::sfLDOF)$upper$bound
simulation |> summary(bound = bound)

# Summarize simulation and compare with the planned design
simulation |> summary(design = design)
 


FUNCTION: to_sim_pw_surv 
TITLE: Convert enrollment and failure rates from \code{sim_fixed_n()} to
\code{sim_pw_surv()} format 
DESCRIPTION: 
\code{to_sim_pw_surv()} converts failure rates and dropout rates entered in
the simpler format for \code{\link[=sim_fixed_n]{sim_fixed_n()}} to that used for \code{\link[=sim_pw_surv]{sim_pw_surv()}}.
The \code{fail_rate} argument for \code{\link[=sim_fixed_n]{sim_fixed_n()}} requires enrollment rates,
failure rates hazard ratios and dropout rates by stratum for a 2-arm trial,
\code{\link[=sim_pw_surv]{sim_pw_surv()}} is in a more flexible but less obvious but more flexible
format. Since \code{\link[=sim_fixed_n]{sim_fixed_n()}} automatically analyzes data and \code{\link[=sim_pw_surv]{sim_pw_surv()}}
just produces a simulation dataset, the latter provides additional options
to analyze or otherwise evaluate individual simulations in ways that
\code{\link[=sim_fixed_n]{sim_fixed_n()}} does not.
 
ARGUMENTS: fail_rate Piecewise constant control group failure rates,
hazard ratio for experimental vs. control,
and dropout rates by stratum and time period. 
EXAMPLE: 
# Example 1
# Convert standard input
to_sim_pw_surv()

# Stratified example
fail_rate <- data.frame(
  stratum = c(rep("Low", 3), rep("High", 3)),
  duration = rep(c(4, 10, 100), 2),
  fail_rate = c(
    .04, .1, .06,
    .08, .16, .12
  ),
  hr = c(
    1.5, .5, 2 / 3,
    2, 10 / 16, 10 / 12
  ),
  dropout_rate = .01
)

x <- to_sim_pw_surv(fail_rate)

# Do a single simulation with the above rates
# Enroll 300 patients in ~12 months at constant rate
sim <- sim_pw_surv(
  n = 300,
  stratum = data.frame(stratum = c("Low", "High"), p = c(.6, .4)),
  enroll_rate = data.frame(duration = 12, rate = 300 / 12),
  fail_rate = x$fail_rate,
  dropout_rate = x$dropout_rate
)

# Cut after 200 events and do a stratified logrank test
sim |>
  cut_data_by_event(200) |> # Cut data
  wlr(weight = fh(rho = 0, gamma = 0)) # Stratified logrank
 


FUNCTION: wlr 
TITLE: Weighted logrank test 
DESCRIPTION: 
Weighted logrank test
 
ARGUMENTS: data Dataset (generated by \code{\link[=sim_pw_surv]{sim_pw_surv()}}) that has been cut by
\code{\link[=counting_process]{counting_process()}}, \code{\link[=cut_data_by_date]{cut_data_by_date()}}, or \code{\link[=cut_data_by_event]{cut_data_by_event()}}.
weight Weighting functions, such as \code{\link[=fh]{fh()}}, \code{\link[=mb]{mb()}}, and
\code{\link[=early_zero]{early_zero()}}.
return_variance A logical flag that, if \code{TRUE}, adds columns
estimated variance for weighted sum of observed minus expected;
see details; Default: \code{FALSE}.
ratio randomization ratio (experimental:control).
\itemize{
\item If the \code{data} is generated by simtrial, such as
\itemize{
\item \code{data = sim_pw_surv(...) |> cut_data_by_date(...)}
\item \code{data = sim_pw_surv(...) |> cut_data_by_event(...)}
\item \code{data = sim_pw_surv(...) |> cut_data_by_date(...) |> counting_process(...)}
\item \code{data = sim_pw_surv(...) |> cut_data_by_event(...) |> counting_process(...)}
there is no need to input the \code{ratio}, as simtrial gets the \code{ratio} via the
\code{block} arguments in \code{\link[=sim_pw_surv]{sim_pw_surv()}}.
}
\item If the \code{data} is a custom dataset (see Example 2) below,
\itemize{
\item Users are suggested to input the planned randomization ratio to \code{ratio};
\item If not, simtrial takes the empirical randomization ratio.
}
}
formula A formula to specify the columns that contain the
time-to-event, event, treatment, and stratum variables. Only used by the
default S3 method because the other classes aleady have the required column
names. For stratified designs, the formula should have the form \code{Surv(tte, event) ~ treatment + strata(stratum)}, where \code{tte}, \code{event}, \code{treatment},
and \code{stratum} are the column names from \code{data} with the time-to-event
measurement, event status, treatment group, and stratum, respectively. For
unstratified designs, the formula can omit the stratum column: \code{Surv(tte, event) ~ treatment}. 
EXAMPLE: 
# ---------------------- #
#      Example 1         #
#  Use dataset generated #
#     by simtrial        #
# ---------------------- #
x <- sim_pw_surv(n = 200) |> cut_data_by_event(100)

# Example 1A: WLR test with FH wights
x |> wlr(weight = fh(rho = 0, gamma = 0.5))
x |> wlr(weight = fh(rho = 0, gamma = 0.5), return_variance = TRUE)

# Example 1B: WLR test with MB wights
x |> wlr(weight = mb(delay = 4, w_max = 2))

# Example 1C: WLR test with early zero wights
x |> wlr(weight = early_zero(early_period = 4))

# Example 1D
# For increased computational speed when running many WLR tests, you can
# pre-compute the counting_process() step first, and then pass the result of
# counting_process() directly to wlr()
x <- x |> counting_process(arm = "experimental")
x |> wlr(weight = fh(rho = 0, gamma = 1))
x |> wlr(weight = mb(delay = 4, w_max = 2))
x |> wlr(weight = early_zero(early_period = 4))

# ---------------------- #
#      Example 2         #
#  Use cumsum dataset    #
# ---------------------- #
x <- data.frame(treatment = ifelse(ex1_delayed_effect$trt == 1, "experimental", "control"),
                stratum = rep("All", nrow(ex1_delayed_effect)),
                tte = ex1_delayed_effect$month,
                event = ex1_delayed_effect$evntd)

# Users can specify the randomization ratio to calculate the statistical information under H0
x |> wlr(weight = fh(rho = 0, gamma = 0.5), ratio = 2)

x |>
  counting_process(arm = "experimental") |>
  wlr(weight = fh(rho = 0, gamma = 0.5), ratio = 2)

# If users don't provide the randomization ratio, we will calculate the emperical ratio
x |> wlr(weight = fh(rho = 0, gamma = 0.5))

x |>
  counting_process(arm = "experimental") |>
  wlr(weight = fh(rho = 0, gamma = 0.5))

# ---------------------- #
#      Example 3         #
#  Use formula           #
# ---------------------- #
library("survival")

# Unstratified design
x <- sim_pw_surv(n = 200) |> cut_data_by_event(100) |> as.data.frame()
colnames(x) <- c("tte", "evnt", "strtm", "trtmnt")
wlr(x, weight = fh(0, 0.5), formula = Surv(tte, evnt) ~ trtmnt)

# Stratified design
x$strtm <- sample(c("s1", "s2"), size = nrow(x), replace = TRUE)
wlr(x, weight = fh(0, 0.5), formula = Surv(tte, evnt) ~ trtmnt + strata(strtm))

View detailed documentation

Overview

simtrial provides comprehensive simulation capabilities for clinical trials with time-to-event endpoints. It enables researchers to validate trial designs, assess operating characteristics, and optimize study parameters through Monte Carlo simulation.

Installation

# Install from CRAN
install.packages("simtrial")

# Install development version from GitHub
devtools::install_github("Merck/simtrial")

Key Features

Flexible Simulation Framework: Simulate various trial designs and scenarios
Time-to-Event Focus: Specialized for survival and time-to-event endpoints
Operating Characteristics: Comprehensive evaluation of design performance
Scenario Analysis: Test multiple assumptions and design variations
Integration: Works seamlessly with gsDesign and gsDesign2

Basic Usage

Simple Simulation

library(simtrial)

# Define trial parameters
enroll_rate <- define_enroll_rate(duration = 18, rate = 20)
fail_rate <- define_fail_rate(
  duration = c(4, 100),
  fail_rate = log(2) / 12,
  hr = c(1, 0.7),
  dropout_rate = 0.001
)

# Run simulation
sim_results <- sim_pw_surv(
  n_sim = 1000,
  study_duration = 36,
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  block = c("Experimental", "Control")
)

Group Sequential Simulation

# Simulate group sequential trial
gs_sim <- sim_gs_n(
  n_sim = 1000,
  alpha = 0.025,
  beta = 0.1,
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  timing = c(0.5, 0.75, 1),
  upper = list(bound = "OF"),
  lower = list(bound = "futility")
)

Advanced Features

Non-Proportional Hazards Simulation

# Simulate with non-proportional hazards
nph_sim <- sim_pw_surv(
  n_sim = 1000,
  study_duration = 36,
  enroll_rate = define_enroll_rate(duration = 18, rate = 20),
  fail_rate = define_fail_rate(
    duration = c(4, 8, 100),
    fail_rate = log(2) / 12,
    hr = c(1, 0.8, 0.6),
    dropout_rate = 0.001
  )
)

Adaptive Designs

# Simulate adaptive trial with interim modifications
adaptive_sim <- sim_fixed_n(
  n_sim = 1000,
  alpha = 0.025,
  power = 0.9,
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  study_duration = 36,
  interim_analysis = c(12, 24)
)

Simulation Outputs

Key Metrics

Power: Probability of rejecting null hypothesis
Type I Error: False positive rate under null hypothesis
Expected Study Duration: Average trial completion time
Expected Sample Size: Average number of subjects needed
Operating Characteristics: Comprehensive performance metrics

Result Analysis

# Analyze simulation results
summary(sim_results)
plot(sim_results)

# Calculate operating characteristics
oc <- operating_characteristics(sim_results)
print(oc)

Best Practices

Sufficient Simulations: Use adequate number of simulations (≥1000)
Scenario Testing: Test multiple scenarios and assumptions
Sensitivity Analysis: Assess robustness to assumption changes
Computational Efficiency: Use parallel processing for large simulations
Result Validation: Cross-check with analytical methods when possible

Common Applications

Power analysis for complex designs
Operating characteristics evaluation
Adaptive trial planning
Non-proportional hazards assessment
Regulatory submission support

Integration with Other Packages

simtrial works well with:

gsDesign2: For modern group sequential designs
gsDesign: For traditional group sequential methods
survival: For survival analysis
future: For parallel processing
ggplot2: For result visualization