Implement sim_gs_n()

DESCRIPTION

@@ -44,6 +44,7 @@ Imports:
     future,
     methods,
     mvtnorm,
+    rlang,
     stats,
     survival,
     utils

NAMESPACE

@@ -3,20 +3,26 @@
 export(counting_process)
 export(cut_data_by_date)
 export(cut_data_by_event)
+export(early_zero)
 export(early_zero_weight)
+export(fh)
 export(fh_weight)
 export(fit_pwexp)
 export(get_analysis_date)
 export(get_cut_date_by_event)
+export(maxcombo)
+export(mb)
 export(mb_weight)
 export(pvalue_maxcombo)
 export(randomize_by_fixed_block)
 export(rmst)
 export(rpwexp)
 export(rpwexp_enroll)
 export(sim_fixed_n)
+export(sim_gs_n)
 export(sim_pw_surv)
 export(to_sim_pw_surv)
+export(wlr)
 importFrom(Rcpp,sourceCpp)
 importFrom(data.table,":=")
 importFrom(data.table,.N)

R/maxcombo.R

@@ -0,0 +1,57 @@
+#  Copyright (c) 2023 Merck & Co., Inc., Rahway, NJ, USA and its affiliates.
+#  All rights reserved.
+#
+#  This file is part of the simtrial program.
+#
+#  simtrial is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#' Maxcombo test
+#'
+#' @param data a tte dataset
+#' @param test1 maxcombo test1
+#' @param test2 maxcombo test2
+#' @param ... additional tests
+#'
+#' @return pvalues
+#' @export
+#'
+#' @examples
+#' sim_pw_surv(n = 200) |>
+#'   cut_data_by_event(150) |>
+#'   maxcombo(test1 = wlr(data, rho = 0, gamma = 0) |> quote(),
+#'            test2 = wlr(data, rho = 0, gamma = 0.5) |> quote())
+maxcombo <- function(data, test1, test2, ...){
+  all_args <- match.call(expand.dots = FALSE)
+  args <- all_args[-1]  # Exclude the first element (function name)
+
+  n_test <- length(args) - 1
+  rho_vector <- NULL
+  gamma_vector <- NULL
+
+  for (i in 1:n_test) {
+    test_i <- get(paste0("test", i))
+    rho_vector <- c(rho_vector, test_i$rho)
+    gamma_vector <- c(gamma_vector, test_i$gamma)
+  }
+
+  ans <- data |>
+    counting_process(arm = "experimental") |>
+    fh_weight(
+      rho_gamma = data.frame(rho = rho_vector, gamma = gamma_vector),
+      return_corr = TRUE)
+
+  ans <-  data.frame(p_value = pvalue_maxcombo(ans))
+
+  return(ans)
+}

R/sim_gs_n.R

@@ -0,0 +1,222 @@
+#  Copyright (c) 2023 Merck & Co., Inc., Rahway, NJ, USA and its affiliates.
+#  All rights reserved.
+#
+#  This file is part of the simtrial program.
+#
+#  simtrial is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#' Simulate group sequantial designs with fixed sample size
+#' @inheritParams sim_fixed_n
+#' @param test a functional call of the test such as \code{wlr()} or \code{maxcombo()}
+#' @param cutting a functional call of the cutting for IA(s) and FA, see examples
+#' @param seed random seed
+#'
+#' @return a data frame summaring the simulation ID, analysis date, z statistics or p-values
+#' @export
+#'
+#' @examples
+#' 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_col <- 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_col,
+#'                               hr = median_col / 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 (exp:col)
+#'
+#' # 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 <- get_analysis_date(data,
+#'                          planned_calendar_time = 20,
+#'                          target_event_overall = 100,
+#'                          max_extension_for_target_event = 24,
+#'                          min_n_overall = 200,
+#'                          min_followup = 20) |> quote()
+#'
+#' # 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 <- get_analysis_date(data,
+#'                          planned_calendar_time = 32,
+#'                          target_event_overall = 200,
+#'                          max_extension_for_target_event = 34,
+#'                          min_time_after_previous_analysis = 10) |> quote()
+#'
+#' # 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 <- get_analysis_date(data,
+#'                        planned_calendar_time = 45,
+#'                        target_event_overall = 350) |> quote()
+#'
+#' # Test 1: regular logrank test
+#' sim_gs_n(
+#'   n_sim = 3,
+#'   sample_size = 400,
+#'   enroll_rate = enroll_rate,
+#'   fail_rate = fail_rate,
+#'   test = wlr(data, weight = fh(rho = 0, gamma = 0)) |> quote(),
+#'   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
+#'   seed = 2024)
+#'
+#' # Test 2: weighted logrank test by FH(0, 0.5)
+#' sim_gs_n(
+#'   n_sim = 3,
+#'   sample_size = 400,
+#'   enroll_rate = enroll_rate,
+#'   fail_rate = fail_rate,
+#'   test = wlr(data, weight = fh(rho = 0, gamma = 0.5)) |> quote(),
+#'   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
+#'   seed = 2024)
+#'
+#'
+#' # Test 3: weighted logrank test by MB(6)
+#' sim_gs_n(
+#'   n_sim = 3,
+#'   sample_size = 400,
+#'   enroll_rate = enroll_rate,
+#'   fail_rate = fail_rate,
+#'   test = wlr(data, weight = mb(delay = 3)) |> quote(),
+#'   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
+#'   seed = 2024)
+#'
+#' # Test 4: weighted logrank test by early zero (6)
+#' sim_gs_n(
+#'   n_sim = 3,
+#'   sample_size = 400,
+#'   enroll_rate = enroll_rate,
+#'   fail_rate = fail_rate,
+#'   test = wlr(data, weight = early_zero(6)) |> quote(),
+#'   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
+#'   seed = 2024)
+#'
+#' # Test 5: RMST
+#' sim_gs_n(
+#'   n_sim = 3,
+#'   sample_size = 400,
+#'   enroll_rate = enroll_rate,
+#'   fail_rate = fail_rate,
+#'   test = rmst(data, tau = 20) |> quote(),
+#'   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
+#'   seed = 2024)
+#'
+#' # Test 6: maxcombo (FH(0,0) + FH(0, 0.5))
+#' sim_gs_n(
+#'   n_sim = 3,
+#'   sample_size = 400,
+#'   enroll_rate = enroll_rate,
+#'   fail_rate = fail_rate,
+#'   test = maxcombo(data,
+#'                   test1 = wlr(data, rho = 0, gamma = 0) |> quote(),
+#'                   test2 = wlr(data, rho = 0, gamma = 0.5) |> quote()) |> quote(),
+#'   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
+#'   seed = 2024)
+sim_gs_n <- function(
+    # number of simulations
+  n_sim = 1000,
+  # sample size
+  sample_size = 500,
+  # multinomial probability distribution for stratum enrollment
+  stratum = data.frame(stratum = "All", p = 1),
+  # enrollment rates
+  enroll_rate = data.frame(duration = c(2, 2, 10), rate = c(3, 6, 9)),
+  # failure rates
+  fail_rate = data.frame(
+    stratum = "All",
+    duration = c(3, 100),
+    fail_rate = log(2) / c(9, 18),
+    hr = c(.9, .6),
+    dropout_rate = rep(.001, 2)
+  ),
+  # fixed block randomization specification
+  block = rep(c("experimental", "control"), 2),
+  # default is to to logrank testing
+  # but alternative tests (such as rmst, maxcombo) can be specified
+  test = wlr(weight = fh(rho = 0, gamma = 0)) |> quote(),
+  # cutting for IA(s) and FA
+  cutting = NULL,
+  # random seed
+  seed = 2024
+){
+  # input checking
+  # TODO
+
+  # simulate for n_sim times
+  ans <- NULL
+  for (sim_id in 1:n_sim) {
+    set.seed(seed + sim_id)
+    # generate data
+    simu_data <- sim_pw_surv(
+      n = sample_size,
+      stratum = stratum,
+      block = block,
+      enroll_rate = enroll_rate,
+      fail_rate = to_sim_pw_surv(fail_rate)$fail_rate,
+      dropout_rate = to_sim_pw_surv(fail_rate)$dropout_rate)
+
+    # initialize the cut date of IA(s) and FA
+    n_analysis <- length(cutting)
+    cut_date <- rep(-100, n_analysis)
+    ans_1sim <- NULL
+
+    for (i_analysis in 1:n_analysis) {
+
+      # get cut date
+      if (i_analysis < n_analysis) {
+        cut_date[i_analysis] <- cutting[[paste0("ia", i_analysis)]] |> eval(envir = rlang::env(data = simu_data))
+      } else {
+        cut_date[i_analysis] <- cutting[["fa"]] |> eval(envir = rlang::env(data = simu_data))
+      }
+
+      # cut the data
+      simu_data_cut <- simu_data |> cut_data_by_date(cut_date[i_analysis])
+
+      # test
+      ans_1sim_new <- eval(test, envir = rlang::env(data = simu_data_cut))
+      ans_1sim_new$analysis <- i_analysis
+      ans_1sim_new$cut_date <- cut_date[i_analysis]
+      ans_1sim_new$sim_id <- sim_id
+
+      # rbind simulation results for all IA(s) and FA in 1 simulation
+      ans_1sim <- rbind(ans_1sim, ans_1sim_new)
+    }
+
+    ans <- rbind(ans, ans_1sim)
+  }
+  return(ans)
+}

R/wlr.R

@@ -0,0 +1,66 @@
+#  Copyright (c) 2023 Merck & Co., Inc., Rahway, NJ, USA and its affiliates.
+#  All rights reserved.
+#
+#  This file is part of the simtrial program.
+#
+#  simtrial is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#' Weighted logrank test
+#'
+#' @param data cutted dataset generated by sim_pw_surv
+#' @param weight weighting functions, such as \code{fh_weight}, \code{mb_weight}, and \code{early_zero_weight}.
+#'
+#' @return test results
+#' @export
+#' @examples
+#' sim_pw_surv(n = 200) |>
+#'   cut_data_by_event(150) |>
+#'   wlr(weight = fh(rho = 0, gamma = 0))
+#'
+#' sim_pw_surv(n = 200) |>
+#'   cut_data_by_event(150) |>
+#'   wlr(weight = mb(delay = 4, w_max = 2))
+#'
+#' sim_pw_surv(n = 200) |>
+#'   cut_data_by_event(150) |>
+#'   wlr(weight = early_zero(early_period = 4))
+#'
+wlr <- function(data, weight){
+
+  if ("fh" %in% class(weight)) {
+    ans <- data |>
+      counting_process(arm = "experimental") |>
+      fh_weight(rho_gamma = data.frame(rho = weight$rho, gamma = weight$gamma))
+
+  } else if ("mb" %in% class(weight)) {
+    ans <- data |>
+      counting_process(arm = "experimental") |>
+      mb_weight(delay = weight$delay, w_max = weight$w_max) |>
+      dplyr::summarize(
+        s = sum(o_minus_e * mb_weight),
+        v = sum(var_o_minus_e * mb_weight^2),
+        z = s / sqrt(v)) |>
+      dplyr::select(z)
+  } else if ("early_period" %in% class(weight)){
+    ans <- data |>
+      counting_process(arm = "experimental") |>
+      early_zero_weight(early_period = weight$early_period) |>
+      dplyr::summarize(
+        s = sum(o_minus_e * weight),
+        v = sum(var_o_minus_e * weight^2),
+        z = s / sqrt(v)) |>
+      dplyr::select(z)
+  }
+  return(ans)
+}