| Title: | Clinical Trial Simulation |
|---|---|
| Description: | 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) <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. |
| Authors: | Keaven Anderson [aut], Yujie Zhao [ctb, cre], John Blischak [ctb], Nan Xiao [ctb], Yilong Zhang [aut], Jianxiao Yang [ctb], Lili Ling [ctb], Xintong Li [ctb], Ruixue Wang [ctb], Yi Cui [ctb], Ping Yang [ctb], Yalin Zhu [ctb], Heng Zhou [ctb], Amin Shirazi [ctb], Cole Manschot [ctb], Merck & Co., Inc., Rahway, NJ, USA and its affiliates [cph] |
| Maintainer: | Yujie Zhao <[email protected]> |
| License: | GPL-3 |
| Version: | 0.4.2 |
| Built: | 2024-11-18 15:27:28 UTC |
| Source: | https://github.com/merck/simtrial |
Convert summary table to a gt object
as_gt(x, ...) ## S3 method for class 'simtrial_gs_wlr' as_gt( x, title = "Summary of simulation results by WLR tests", subtitle = NULL, ... )as_gt(x, ...) ## S3 method for class 'simtrial_gs_wlr' as_gt( x, title = "Summary of simulation results by WLR tests", subtitle = NULL, ... )
x |
A object returned by |
... |
Additional parameters (not used). |
title |
Title of the gt table. |
subtitle |
Subtitle of the gt table. |
A gt table.
A gt table summarizing the simulation results.
# 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()# 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()
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.
counting_process(x, arm)counting_process(x, arm)
x |
A data frame with no missing values and contain variables:
|
arm |
Value in the input |
The function only considered two group situation.
The tie is handled by the Breslow's Method.
The output produced by counting_process() produces a
counting process dataset grouped by stratum and sorted within stratum
by increasing times where events occur. The object is assigned the class
"counting_process". It also has the attributes "n_ctrl" and "n_exp",
which are the totals of the control and experimental treatments,
respectively, from the input time-to-event data.
A data frame grouped by stratum and sorted within stratum by tte.
Remain rows with at least one event in the population, at least one subject
is at risk in both treatment group and control group.
Other variables in this represent the following within each stratum at
each time at which one or more events are observed:
event_total: Total number of events
event_trt: Total number of events at treatment group
n_risk_total: Number of subjects at risk
n_risk_trt: Number of subjects at risk in treatment group
s: Left-continuous Kaplan-Meier survival estimate
o_minus_e: In treatment group, observed number of events minus expected
number of events. The expected number of events is estimated by assuming
no treatment effect with hypergeometric distribution with parameters total
number of events, total number of events at treatment group and number of
events at a time. (Same assumption of log-rank test under the null
hypothesis)
var_o_minus_e: Variance of o_minus_e under the same assumption.
# 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))# 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))
Create a cutting function for use with sim_gs_n()
create_cut(...)create_cut(...)
... |
Arguments passed to |
A function that accepts a data frame of simulated trial data and returns a cut date
get_analysis_date(), sim_gs_n()
# 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)# 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)
Create a cutting test function for use with sim_gs_n()
create_test(test, ...)create_test(test, ...)
test |
A test function such as |
... |
Arguments passed to the cutting test function |
A function that accepts a data frame of simulated trial data and returns a test result
# 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)) ))# 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)) ))
Cut a dataset for analysis at a specified date
cut_data_by_date(x, cut_date)cut_data_by_date(x, cut_date)
x |
A time-to-event dataset, for example, generated by |
cut_date |
Date relative to start of randomization
( |
A dataset ready for survival analysis.
# 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)# 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)
Takes a time-to-event data set and cuts the data at which an event count is reached.
cut_data_by_event(x, event)cut_data_by_event(x, event)
x |
A time-to-event dataset, for example, generated by |
event |
Event count at which data cutoff is to be made. |
A data frame ready for survival analysis, including columns
time to event (tte), event, the stratum, and the treatment.
# 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)# 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)
Zero early weighting function
early_zero(early_period, fail_rate = NULL)early_zero(early_period, fail_rate = NULL)
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 |
A list of parameters of the zero early weighting function
Xu, Z., Zhen, B., Park, Y., & Zhu, B. (2017). "Designing therapeutic cancer vaccine trials with delayed treatment effect."
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))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))
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.
data(ex1_delayed_effect)data(ex1_delayed_effect)
Data frame with 4 variables:
id: Sequential numbering of unique identifiers.
month: Time-to-event.
event: 1 for event, 0 for censored.
trt: 1 for experimental, 0 for control.
Lin, Ray S., Ji Lin, Satrajit Roychoudhury, Keaven M. Anderson, Tianle Hu, Bo Huang, Larry F Leon, Jason J.Z. Liao, Rong Liu, Xiaodong Luo, Pralay Mukhopadhyay, Rui Qin, Kay Tatsuoka, Xuejing Wang, Yang Wang, Jian Zhu, Tai-Tsang Chen, Renee Iacona & Cross-Pharma Non-proportional Hazards Working Group. 2020. Alternative analysis methods for time to event endpoints under nonproportional hazards: A comparative analysis. Statistics in Biopharmaceutical Research 12(2): 187–198.
ex2_delayed_effect, ex3_cure_with_ph, ex4_belly, ex5_widening, ex6_crossing
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))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))
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.
data(ex2_delayed_effect)data(ex2_delayed_effect)
Data frame with 4 variables:
id: Sequential numbering of unique identifiers.
month: Time-to-event.
event: 1 for event, 0 for censored.
trt: 1 for experimental, 0 for control.
Lin, Ray S., Ji Lin, Satrajit Roychoudhury, Keaven M. Anderson, Tianle Hu, Bo Huang, Larry F Leon, Jason J.Z. Liao, Rong Liu, Xiaodong Luo, Pralay Mukhopadhyay, Rui Qin, Kay Tatsuoka, Xuejing Wang, Yang Wang, Jian Zhu, Tai-Tsang Chen, Renee Iacona & Cross-Pharma Non-proportional Hazards Working Group. 2020. Alternative analysis methods for time to event endpoints under nonproportional hazards: A comparative analysis. Statistics in Biopharmaceutical Research 12(2): 187–198.
ex1_delayed_effect, ex3_cure_with_ph, ex4_belly, ex5_widening, ex6_crossing
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))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))
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.
data(ex3_cure_with_ph)data(ex3_cure_with_ph)
Data frame with 4 variables:
id: Sequential numbering of unique identifiers.
month: Time-to-event.
event: 1 for event, 0 for censored.
trt: 1 for experimental, 0 for control.
Lin, Ray S., Ji Lin, Satrajit Roychoudhury, Keaven M. Anderson, Tianle Hu, Bo Huang, Larry F Leon, Jason J.Z. Liao, Rong Liu, Xiaodong Luo, Pralay Mukhopadhyay, Rui Qin, Kay Tatsuoka, Xuejing Wang, Yang Wang, Jian Zhu, Tai-Tsang Chen, Renee Iacona & Cross-Pharma Non-proportional Hazards Working Group. 2020. Alternative analysis methods for time to event endpoints under nonproportional hazards: A comparative analysis. Statistics in Biopharmaceutical Research 12(2): 187–198.
ex1_delayed_effect, ex2_delayed_effect, ex4_belly, ex5_widening, ex6_crossing
library(survival) data(ex3_cure_with_ph) km1 <- with(ex3_cure_with_ph, survfit(Surv(month, evntd) ~ trt)) km1 plot(km1)library(survival) data(ex3_cure_with_ph) km1 <- with(ex3_cure_with_ph, survfit(Surv(month, evntd) ~ trt)) km1 plot(km1)
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.
data(ex4_belly)data(ex4_belly)
Data frame with 4 variables:
id: Sequential numbering of unique identifiers.
month: Time-to-event.
event: 1 for event, 0 for censored.
trt: 1 for experimental, 0 for control.
Lin, Ray S., Ji Lin, Satrajit Roychoudhury, Keaven M. Anderson, Tianle Hu, Bo Huang, Larry F Leon, Jason J.Z. Liao, Rong Liu, Xiaodong Luo, Pralay Mukhopadhyay, Rui Qin, Kay Tatsuoka, Xuejing Wang, Yang Wang, Jian Zhu, Tai-Tsang Chen, Renee Iacona & Cross-Pharma Non-proportional Hazards Working Group. 2020. Alternative analysis methods for time to event endpoints under nonproportional hazards: A comparative analysis. Statistics in Biopharmaceutical Research 12(2): 187–198.
ex1_delayed_effect, ex2_delayed_effect, ex3_cure_with_ph, ex5_widening, ex6_crossing
library(survival) data(ex4_belly) km1 <- with(ex4_belly, survfit(Surv(month, evntd) ~ trt)) km1 plot(km1)library(survival) data(ex4_belly) km1 <- with(ex4_belly, survfit(Surv(month, evntd) ~ trt)) km1 plot(km1)
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.
data(ex5_widening)data(ex5_widening)
Data frame with 4 variables:
id: Sequential numbering of unique identifiers.
month: Time-to-event.
event: 1 for event, 0 for censored.
trt: 1 for experimental, 0 for control.
Lin, Ray S., Ji Lin, Satrajit Roychoudhury, Keaven M. Anderson, Tianle Hu, Bo Huang, Larry F Leon, Jason J.Z. Liao, Rong Liu, Xiaodong Luo, Pralay Mukhopadhyay, Rui Qin, Kay Tatsuoka, Xuejing Wang, Yang Wang, Jian Zhu, Tai-Tsang Chen, Renee Iacona & Cross-Pharma Non-proportional Hazards Working Group. 2020. Alternative analysis methods for time to event endpoints under nonproportional hazards: A comparative analysis. Statistics in Biopharmaceutical Research 12(2): 187–198.
ex1_delayed_effect, ex2_delayed_effect, ex3_cure_with_ph, ex4_belly, ex6_crossing
library(survival) data(ex5_widening) km1 <- with(ex5_widening, survfit(Surv(month, evntd) ~ trt)) km1 plot(km1)library(survival) data(ex5_widening) km1 <- with(ex5_widening, survfit(Surv(month, evntd) ~ trt)) km1 plot(km1)
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.
data(ex6_crossing)data(ex6_crossing)
Data frame with 4 variables:
id: Sequential numbering of unique identifiers.
month: Time-to-event.
event: 1 for event, 0 for censored.
trt: 1 for experimental, 0 for control.
Lin, Ray S., Ji Lin, Satrajit Roychoudhury, Keaven M. Anderson, Tianle Hu, Bo Huang, Larry F Leon, Jason J.Z. Liao, Rong Liu, Xiaodong Luo, Pralay Mukhopadhyay, Rui Qin, Kay Tatsuoka, Xuejing Wang, Yang Wang, Jian Zhu, Tai-Tsang Chen, Renee Iacona & Cross-Pharma Non-proportional Hazards Working Group. 2020. Alternative analysis methods for time to event endpoints under nonproportional hazards: A comparative analysis. Statistics in Biopharmaceutical Research 12(2): 187–198.
ex1_delayed_effect, ex2_delayed_effect, ex3_cure_with_ph, ex4_belly, ex5_widening
library(survival) data(ex6_crossing) km1 <- with(ex6_crossing, survfit(Surv(month, evntd) ~ trt)) km1 plot(km1)library(survival) data(ex6_crossing) km1 <- with(ex6_crossing, survfit(Surv(month, evntd) ~ trt)) km1 plot(km1)
Fleming-Harrington weighting function
fh(rho = 0, gamma = 0)fh(rho = 0, gamma = 0)
rho |
Non-negative number. |
gamma |
Non-negative number. |
A list of parameters of the Fleming-Harrington weighting function
sim_pw_surv(n = 200) |> cut_data_by_event(100) |> wlr(weight = fh(rho = 0, gamma = 1))sim_pw_surv(n = 200) |> cut_data_by_event(100) |> wlr(weight = fh(rho = 0, gamma = 1))
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.
fit_pwexp( srv = Surv(time = ex1_delayed_effect$month, event = ex1_delayed_effect$evntd), intervals = array(3, 3) )fit_pwexp( srv = Surv(time = ex1_delayed_effect$month, event = ex1_delayed_effect$evntd), intervals = array(3, 3) )
srv |
Input survival object (see |
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. |
A matrix with rows containing interval length, estimated rate, -2 * log-likelihood for each interval.
# 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)# 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)
Derive analysis date for interim/final analysis given multiple conditions
get_analysis_date( data, planned_calendar_time = NA, target_event_overall = NA, target_event_per_stratum = NA, max_extension_for_target_event = NA, previous_analysis_date = 0, min_time_after_previous_analysis = NA, min_n_overall = NA, min_n_per_stratum = NA, min_followup = NA )get_analysis_date( data, planned_calendar_time = NA, target_event_overall = NA, target_event_per_stratum = NA, max_extension_for_target_event = NA, previous_analysis_date = 0, min_time_after_previous_analysis = NA, min_n_overall = NA, min_n_per_stratum = NA, min_followup = NA )
data |
A simulated data generated by |
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
|
To obtain the analysis date, consider the following multiple conditions:
The planned calendar time for analysis.
The targeted events, encompassing both overall population and stratum-specific events.
The maximum time extension required to achieve the targeted events.
The planned minimum time interval after the previous analysis.
The minimum follow-up time needed to reach a certain number of patients in enrollments.
Users have the flexibility to employ all 5 conditions simultaneously or
selectively choose specific conditions to determine the analysis date.
Any unused conditions will default to NA and not affect the output.
Regardless of the number of conditions used, the analysis date is determined
by min(max(date1, date2, date4, date5, na.rm = TRUE), date3, na.rm = TRUE),
where date1, date2, date3, date4, date5 represent the analysis
dates determined solely by Condition 1, Condition 2, Condition 3,
Condition 4 and Condition 5, respectively.
A numerical value of the analysis date.
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 )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 )
Get date at which an event count is reached
get_cut_date_by_event(x, event)get_cut_date_by_event(x, event)
x |
A time-to-event dataset, for example, generated by |
event |
Event count at which dataset is to be cut off for analysis. |
A numeric value with the cte from the input dataset
at which the targeted event count is reached, or if the final event count
is never reached, the final cte at which an event occurs.
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)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)
WARNING: This experimental function is a work-in-progress. The function arguments will change as we add additional features.
maxcombo( data = cut_data_by_event(sim_pw_surv(n = 200), 150), rho = c(0, 0, 1), gamma = c(0, 1, 1), return_variance = FALSE, return_corr = FALSE )maxcombo( data = cut_data_by_event(sim_pw_surv(n = 200), 150), rho = c(0, 0, 1), gamma = c(0, 1, 1), return_variance = FALSE, return_corr = FALSE )
data |
A TTE dataset. |
rho |
Numeric vector. Must be greater
than or equal to zero. Must be the same length as |
gamma |
Numeric vector. Must be
greater than or equal to zero. Must be the same length as |
return_variance |
A logical flag that, if |
return_corr |
A logical flag that, if |
A list containing the test method (method),
parameters of this test method (parameter),
point estimate of the treatment effect (estimate),
standardized error of the treatment effect (se),
Z-score of each test of the MaxCombo (z),
p-values (p_value)
and the correlation matrix of each tests in MaxCombo (begin with v)
sim_pw_surv(n = 200) |> cut_data_by_event(150) |> maxcombo(rho = c(0, 0), gamma = c(0, 1), return_corr = TRUE)sim_pw_surv(n = 200) |> cut_data_by_event(150) |> maxcombo(rho = c(0, 0), gamma = c(0, 1), return_corr = TRUE)
Magirr and Burman weighting function
mb(delay = 4, w_max = Inf)mb(delay = 4, w_max = Inf)
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 |
Magirr and Burman (2019) proposed a weighted logrank test to have better power than the logrank test when the treatment effect is delayed, but to still maintain good power under a proportional hazards assumption. In Magirr (2021), (the equivalent of) a maximum weight was proposed as opposed to a fixed time duration over which weights would increase. The weights for some early interval specified by the user are the inverse of the combined treatment group empirical survival distribution; see details. After this initial period, weights are constant at the maximum of the previous weights. Another advantage of the test is that under strong null hypothesis that the underlying survival in the control group is greater than or equal to underlying survival in the experimental group, Type I error is controlled as the specified level.
We define to be the input variable delay.
This specifies an initial period during which weights increase.
We also set a maximum weight .
To define specific weights, we let denote the Kaplan-Meier
survival estimate at time for the combined data
(control plus experimental treatment groups).
The weight at time is then defined as
A list of parameters of the Magirr and Burman weighting function
Magirr, Dominic, and Carl‐Fredrik Burman. 2019. "Modestly weighted logrank tests." Statistics in Medicine 38 (20): 3782–3790.
Magirr, Dominic. 2021. "Non‐proportional hazards in immuno‐oncology: Is an old perspective needed?" Pharmaceutical Statistics 20 (3): 512–527.
sim_pw_surv(n = 200) |> cut_data_by_event(100) |> wlr(weight = mb(delay = 8, w_max = Inf))sim_pw_surv(n = 200) |> cut_data_by_event(100) |> wlr(weight = mb(delay = 8, w_max = Inf))
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 sim_pw_surv() function
under the above scenario.
mb_delayed_effectmb_delayed_effect
A data frame with 200 rows and 4 columns:
tte: Time to event.
Magirr, Dominic, and Carl‐Fredrik Burman. 2019. "Modestly weighted logrank tests." Statistics in Medicine 38 (20): 3782–3790.
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)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)
Milestone test for two survival curves
milestone(data, ms_time, test_type = c("log-log", "naive"))milestone(data, ms_time, test_type = c("log-log", "naive"))
data |
Data frame containing at least 3 columns:
|
ms_time |
Milestone analysis time. |
test_type |
Method to build the test statistics. There are 2 options:
|
A list frame containing:
method - The method, always "milestone".
parameter - Milestone time point.
estimate - Survival difference between the experimental and control arm.
se - Standard error of the control and experimental arm.
z - Test statistics.
Klein, J. P., Logan, B., Harhoff, M., & Andersen, P. K. (2007). "Analyzing survival curves at a fixed point in time." Statistics in Medicine, 26(24), 4505–4519.
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")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")
WARNING: This experimental function is a work-in-progress. The function arguments and/or returned output format may change as we add additional features.
multitest(data, ...)multitest(data, ...)
data |
Trial data cut by |
... |
One or more test functions. Use |
A list of test results, one per test. If the test functions are named
in the call to multitest(), the returned list uses the same names.
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 )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 )
Fixed block randomization. The 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 n is not an exact multiple
of the block length.
randomize_by_fixed_block(n = 10, block = c(0, 0, 1, 1))randomize_by_fixed_block(n = 10, block = c(0, 0, 1, 1))
n |
Sample size to be randomized. |
block |
Vector of treatments to be included in each block. |
A treatment group sequence (vector) of length n with
treatments from block permuted within each block having
block size equal to the length of block.
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())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())
RMST difference of 2 arms
rmst( data, tau = 10, var_label_tte = "tte", var_label_event = "event", var_label_group = "treatment", formula = NULL, reference = "control", alpha = 0.05 )rmst( data, tau = 10, var_label_tte = "tte", var_label_event = "event", var_label_group = "treatment", formula = NULL, reference = "control", alpha = 0.05 )
data |
A time-to-event dataset with a column |
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: |
reference |
A group label indicating the reference group. |
alpha |
Type I error. |
The argument formula is provided as a convenience to easily specify the
TTE, event, and grouping variables using the syntax Surv(tte, event) ~ group). Surv() is from the {survival} package (survival::Surv()).
You can also explicitly
name the arguments passed to Surv(), for example the following is
equivalent Surv(event = event, time = tte) ~ group). Note however that the
function Surv() is never actually executed. Similarly, any other functions
applied in the formula are also ignored, thus you shouldn't apply any
transformation functions such as log() since these will have no effect.
The z statistics.
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" )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" )
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.
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 gsDesign::nSurv().
rpwexp(n = 100, fail_rate = data.frame(duration = c(1, 1), rate = c(10, 20)))rpwexp(n = 100, fail_rate = data.frame(duration = c(1, 1), rate = c(10, 20)))
n |
Number of observations to be generated. |
fail_rate |
A data frame containing |
Using the cumulative = TRUE option, enrollment times that piecewise
constant over time can be generated.
The generated random numbers.
# 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}" )# 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}" )
With piecewise exponential enrollment rate generation any enrollment rate
distribution can be easily approximated.
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 gsDesign::nSurv().
rpwexp_enroll( n = NULL, enroll_rate = data.frame(duration = c(1, 2), rate = c(2, 5)) )rpwexp_enroll( n = NULL, enroll_rate = data.frame(duration = c(1, 2), rate = c(2, 5)) )
n |
Number of observations.
Default of |
enroll_rate |
A data frame containing period duration ( |
A vector of random enrollment times.
# 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" )# 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" )
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.
sim_fixed_n( n_sim = 1000, sample_size = 500, target_event = 350, stratum = data.frame(stratum = "All", p = 1), enroll_rate = data.frame(duration = c(2, 2, 10), rate = c(3, 6, 9)), fail_rate = data.frame(stratum = "All", duration = c(3, 100), fail_rate = log(2)/c(9, 18), hr = c(0.9, 0.6), dropout_rate = rep(0.001, 2)), total_duration = 30, block = rep(c("experimental", "control"), 2), timing_type = 1:5, rho_gamma = data.frame(rho = 0, gamma = 0) )sim_fixed_n( n_sim = 1000, sample_size = 500, target_event = 350, stratum = data.frame(stratum = "All", p = 1), enroll_rate = data.frame(duration = c(2, 2, 10), rate = c(3, 6, 9)), fail_rate = data.frame(stratum = "All", duration = c(3, 100), fail_rate = log(2)/c(9, 18), hr = c(0.9, 0.6), dropout_rate = rep(0.001, 2)), total_duration = 30, block = rep(c("experimental", "control"), 2), timing_type = 1:5, rho_gamma = data.frame(rho = 0, gamma = 0) )
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 |
enroll_rate |
Piecewise constant enrollment rates by time period.
Note that these are overall population enrollment rates and the |
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 |
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
|
timing_type has up to 5 elements indicating different options
for data cutoff:
1: Uses the planned study duration.
2: The time the targeted event count is achieved.
3: The planned minimum follow-up after enrollment is complete.
4: The maximum of planned study duration and targeted event count cuts
(1 and 2).
5: The maximum of targeted event count and minimum follow-up cuts
(2 and 3).
A data frame including columns:
event: Event count.
ln_hr: Log-hazard ratio.
z: Normal test statistic; < 0 favors experimental.
cut: Text describing cutoff used.
duration: Duration of trial at cutoff for analysis.
sim: Sequential simulation ID.
One row per simulated dataset per cutoff specified in timing_type,
per test statistic specified.
If multiple Fleming-Harrington tests are specified in rho_gamma,
then columns rho and gamma are also included.
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)) # 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")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)) # 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")
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.
sim_gs_n( n_sim = 1000, sample_size = 500, stratum = data.frame(stratum = "All", p = 1), enroll_rate = data.frame(duration = c(2, 2, 10), rate = c(3, 6, 9)), fail_rate = data.frame(stratum = "All", duration = c(3, 100), fail_rate = log(2)/c(9, 18), hr = c(0.9, 0.6), dropout_rate = rep(0.001, 2)), block = rep(c("experimental", "control"), 2), test = wlr, cut = NULL, ... )sim_gs_n( n_sim = 1000, sample_size = 500, stratum = data.frame(stratum = "All", p = 1), enroll_rate = data.frame(duration = c(2, 2, 10), rate = c(3, 6, 9)), fail_rate = data.frame(stratum = "All", duration = c(3, 100), fail_rate = log(2)/c(9, 18), hr = c(0.9, 0.6), dropout_rate = rep(0.001, 2)), block = rep(c("experimental", "control"), 2), test = wlr, cut = NULL, ... )
n_sim |
Number of simulations to perform. |
sample_size |
Total sample size per simulation. |
stratum |
A data frame with stratum specified in |
enroll_rate |
Piecewise constant enrollment rates by time period.
Note that these are overall population enrollment rates and the |
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 |
test |
One or more test functions such as |
cut |
A list of cutting functions created by |
... |
Arguments passed to the test function(s) provided by the argument
|
WARNING: This experimental function is a work-in-progress. The function arguments will change as we add additional features.
A data frame summarizing the simulation ID, analysis date, z statistics or p-values.
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) ) # 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) ## Not run: 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) ) ## End(Not run) # 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) ) ## Not run: 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) ) ## End(Not run) # 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")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) ) # 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) ## Not run: 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) ) ## End(Not run) # 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) ) ## Not run: 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) ) ## End(Not run) # 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")
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.
sim_pw_surv( n = 100, stratum = data.frame(stratum = "All", p = 1), block = c(rep("control", 2), rep("experimental", 2)), enroll_rate = data.frame(rate = 9, duration = 1), fail_rate = data.frame(stratum = rep("All", 4), period = rep(1:2, 2), treatment = c(rep("control", 2), rep("experimental", 2)), duration = rep(c(3, 1), 2), rate = log(2)/c(9, 9, 9, 18)), dropout_rate = data.frame(stratum = rep("All", 2), period = rep(1, 2), treatment = c("control", "experimental"), duration = rep(100, 2), rate = rep(0.001, 2)) )sim_pw_surv( n = 100, stratum = data.frame(stratum = "All", p = 1), block = c(rep("control", 2), rep("experimental", 2)), enroll_rate = data.frame(rate = 9, duration = 1), fail_rate = data.frame(stratum = rep("All", 4), period = rep(1:2, 2), treatment = c(rep("control", 2), rep("experimental", 2)), duration = rep(c(3, 1), 2), rate = log(2)/c(9, 9, 9, 18)), dropout_rate = data.frame(stratum = rep("All", 2), period = rep(1, 2), treatment = c("control", "experimental"), duration = rep(100, 2), rate = rep(0.001, 2)) )
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 |
block |
Vector of treatments to be included in each block. |
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. |
A data frame with the following variables for each observation:
stratum.
enroll_time: Enrollment time for the observation.
Treatment: Treatment group; this will be one of the values
in the input block.
fail_time: Failure time generated using rpwexp().
dropout_time: Dropout time generated using rpwexp().
cte: Calendar time of enrollment plot the minimum of
failure time and dropout time.
fail: Indicator that cte was set using failure time;
i.e., 1 is a failure, 0 is a dropout.
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 )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 )
Summary of group sequential simulations.
## S3 method for class 'simtrial_gs_wlr' summary(object, design = NULL, bound = NULL, ...)## S3 method for class 'simtrial_gs_wlr' summary(object, design = NULL, bound = NULL, ...)
object |
Simulation results generated by |
design |
Asymptotic design generated by |
bound |
The boundaries. |
... |
Additional parameters (not used). |
A data frame
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)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)
sim_fixed_n() to
sim_pw_surv() formatto_sim_pw_surv() converts failure rates and dropout rates entered in
the simpler format for sim_fixed_n() to that used for sim_pw_surv().
The fail_rate argument for sim_fixed_n() requires enrollment rates,
failure rates hazard ratios and dropout rates by stratum for a 2-arm trial,
sim_pw_surv() is in a more flexible but less obvious but more flexible
format. Since sim_fixed_n() automatically analyzes data and sim_pw_surv()
just produces a simulation dataset, the latter provides additional options
to analyze or otherwise evaluate individual simulations in ways that
sim_fixed_n() does not.
to_sim_pw_surv( fail_rate = data.frame(stratum = "All", duration = c(3, 100), fail_rate = log(2)/c(9, 18), hr = c(0.9, 0.6), dropout_rate = rep(0.001, 2)) )to_sim_pw_surv( fail_rate = data.frame(stratum = "All", duration = c(3, 100), fail_rate = log(2)/c(9, 18), hr = c(0.9, 0.6), dropout_rate = rep(0.001, 2)) )
fail_rate |
Piecewise constant control group failure rates, hazard ratio for experimental vs. control, and dropout rates by stratum and time period. |
A list of two data frame components formatted for
sim_pw_surv(): fail_rate and dropout_rate.
# 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# 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
Weighted logrank test
wlr(data, weight, return_variance = FALSE, ratio = NULL) ## S3 method for class 'tte_data' wlr(data, weight, return_variance = FALSE, ratio = NULL) ## S3 method for class 'counting_process' wlr(data, weight, return_variance = FALSE, ratio = NULL)wlr(data, weight, return_variance = FALSE, ratio = NULL) ## S3 method for class 'tte_data' wlr(data, weight, return_variance = FALSE, ratio = NULL) ## S3 method for class 'counting_process' wlr(data, weight, return_variance = FALSE, ratio = NULL)
data |
Dataset (generated by |
weight |
Weighting functions, such as |
return_variance |
A logical flag that, if |
ratio |
randomization ratio (experimental:control).
|
- Standardized normal Fleming-Harrington weighted logrank test.
- Stratum index.
- Number of distinct times at which events occurred in
stratum .
- Ordered times at which events in stratum
, were observed;
for each observation, represents the time post study entry.
- Total number of events in stratum that occurred
at time .
- Total number of events in stratum in the
experimental treatment group that occurred at time .
- Total number of study subjects in stratum
who were followed for at least duration.
- Expected observations in experimental treatment group
given random selection of from those in
stratum at risk at time .
- Hypergeometric variance for as
produced in Var from counting_process().
- Total number of study subjects in
stratum in the experimental treatment group
who were followed for at least duration .
- Expected observations in experimental group in
stratum at time conditioning on the overall number
of events and at risk populations at that time and sampling at risk
observations without replacement:
- Kaplan-Meier estimate of survival in combined
treatment groups immediately prior to time .
- Real parameters for Fleming-Harrington test.
- Numerator for signed logrank test in stratum
- Variance used in denominator for Fleming-Harrington
weighted logrank tests
The stratified Fleming-Harrington weighted logrank test is then computed as:
A list containing the test method (method),
parameters of this test method (parameter),
point estimate of the treatment effect (estimate),
standardized error of the treatment effect (se),
Z-score (z), p-values (p_value).
# ---------------------- # # 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) class(x) <- c("tte_data", class(x)) # 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 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) class(x) <- c("tte_data", class(x)) # 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))