Custom Fixed Design Simulations: A Tutorial on Writing Code from the Ground Up

cut_date_a <- get_analysis_date(data = uncut_data,
                                planned_calendar_time = 24,
                                target_event_overall = 300)

cut_date_b <- get_analysis_date(data = uncut_data,
                                min_followup = 12,
                                target_event_overall = 300)

cut_date_c <- get_analysis_date(data = uncut_data,
                                max_extension_for_target_event = 12,
                                target_event_overall = 300)

cut_date_d <- get_analysis_date(data = uncut_data,
                                min_n_overall = 100 * 0.8,
                                min_followup = 12)

cut_date <- cut_date_d
cat("The cutoff date is ", round(cut_date, 2))

cut_data <- uncut_data |> cut_data_by_date(cut_date)
cut_data |> head() |> gt() |> tab_header(paste0("An Overview of TTE data Cut at ", round(cut_date, 2), "Months"))

# Logrank test
sim_res_lr <- cut_data |> wlr(weight = fh(rho = 0, gamma = 0))

# weighted logrank test by Fleming-Harrington weights
sim_res_fh <- cut_data |> wlr(weight = fh(rho = 0, gamma = 0.5))

# Modestly weighted logrank test
sim_res_mb <- cut_data |> wlr(weight = mb(delay = Inf, w_max = 2))

# Weighted logrank test by Xu 2017's early zero weights
sim_res_xu <- cut_data |> wlr(weight = early_zero(early_period = 3))

# RMST test
sim_res_rmst <- cut_data |> rmst(tau = 10)

# Milestone test
sim_res_ms <- cut_data |> milestone(ms_time = 10)

# Maxcombo tests comboing multiple weighted logrank test with Fleming-Harrington weights
sim_res_mc <- cut_data |> maxcombo(rho = c(0, 0), gamma = c(0, 0.5))

sim_res <- tribble(
  ~Method, ~Parameter, ~Z, ~Estimate, ~SE, ~`P value`,
  sim_res_lr$method, sim_res_lr$parameter, sim_res_lr$z, sim_res_lr$estimate, sim_res_lr$se, pnorm(-sim_res_lr$z),
  sim_res_fh$method, sim_res_fh$parameter, sim_res_fh$z, sim_res_fh$estimate, sim_res_fh$se, pnorm(-sim_res_fh$z),
  sim_res_mb$method, sim_res_mb$parameter, sim_res_mb$z, sim_res_mb$estimate, sim_res_mb$se, pnorm(-sim_res_mb$z),
  sim_res_xu$method, sim_res_xu$parameter, sim_res_xu$z, sim_res_xu$estimate, sim_res_xu$se, pnorm(-sim_res_xu$z),
  sim_res_rmst$method, sim_res_rmst$parameter|> as.character(), sim_res_rmst$z, sim_res_rmst$estimate, sim_res_rmst$se, pnorm(-sim_res_rmst$z),
  sim_res_ms$method, sim_res_ms$parameter |> as.character(), sim_res_ms$z, sim_res_ms$estimate, sim_res_ms$se, pnorm(-sim_res_ms$z),
  sim_res_mc$method, sim_res_mc$parameter, NA, NA, NA, sim_res_mc$p_value
  ) 

sim_res |> gt() |> tab_header("One Simulation Results")

one_sim <- function(sim_id = 1, 
                    # arguments from Step 1: design characteristic
                    n, stratum, enroll_rate, fail_rate, dropout_rate, block, 
                    # arguments from Step 2； cutting method
                    min_n_overall, min_followup,
                    # arguments from Step 3； testing method
                    fh, mb, xu, rmst, ms, mc
                    ) {
    # Step 1: simulate a time-to-event data
    uncut_data <- sim_pw_surv(
      n = n,
      stratum = stratum,
      block = block,
      enroll_rate = enroll_rate,
      fail_rate = fail_rate,
      dropout_rate = dropout_rate) 
    
    ## Step 2: Cut data
    cut_date <- get_analysis_date(min_n_overall = min_n_overall, min_followup = min_followup, data = uncut_data)
    cut_data <- uncut_data |> cut_data_by_date(cut_date)
    
    # Step 3: Run tests
    sim_res_lr <- cut_data |> wlr(weight = fh(rho = 0, gamma = 0))
    sim_res_fh <- cut_data |> wlr(weight = fh(rho = fh$rho, gamma = fh$gamma))
    sim_res_mb <- cut_data |> wlr(weight = mb(delay = mb$delay, w_max = mb$w_max))
    sim_res_xu <- cut_data |> wlr(weight = early_zero(early_period = xu$early_period))
    sim_res_rmst <- cut_data |> rmst(tau = rmst$tau)
    sim_res_ms <- cut_data |> milestone(ms_time = ms$ms_time)
    sim_res_mc <- cut_data |> maxcombo(rho = mc$rho, gamma = mc$gamma)
    
    sim_res <- tribble(
      ~`Sim ID`, ~Method, ~Parameter, ~Z, ~Estimate, ~SE, ~`P value`,
      sim_id, sim_res_lr$method, sim_res_lr$parameter, sim_res_lr$z, sim_res_lr$estimate, sim_res_lr$se, pnorm(-sim_res_lr$z),
      sim_id, sim_res_fh$method, sim_res_fh$parameter, sim_res_fh$z, sim_res_fh$estimate, sim_res_fh$se, pnorm(-sim_res_fh$z),
      sim_id, sim_res_mb$method, sim_res_mb$parameter, sim_res_mb$z, sim_res_mb$estimate, sim_res_mb$se, pnorm(-sim_res_mb$z),
      sim_id, sim_res_xu$method, sim_res_xu$parameter, sim_res_xu$z, sim_res_xu$estimate, sim_res_xu$se, pnorm(-sim_res_xu$z),
      sim_id, sim_res_rmst$method, sim_res_rmst$parameter|> as.character(), sim_res_rmst$z, sim_res_rmst$estimate, sim_res_rmst$se, pnorm(-sim_res_rmst$z),
      sim_id, sim_res_ms$method, sim_res_ms$parameter |> as.character(), sim_res_ms$z, sim_res_ms$estimate, sim_res_ms$se, pnorm(-sim_res_ms$z),
      sim_id, sim_res_mc$method, sim_res_mc$parameter, NA, NA, NA, sim_res_mc$p_value
  ) 
      
    return(sim_res)
}

set.seed(2025)

plan("multisession", workers = 2)
ans <- foreach(
  sim_id = seq_len(n_sim),
  .errorhandling = "stop",
  .options.future = list(seed = TRUE)
  ) %dofuture% {
    ans_new <- one_sim(
      sim_id = sim_id, 
      # arguments from Step 1: design characteristic
      n = n, 
      stratum = stratum, 
      enroll_rate = enroll_rate, 
      fail_rate = to_sim_pw_surv(fail_rate)$fail_rate, 
      dropout_rate = differential_dropout_rate, 
      block = block, 
      # arguments from Step 2； cutting method
      min_n_overall = 500 * 0.8,
      min_followup = 12,
      # arguments from Step 3； testing method
      fh = list(rho = 0, gamma = 0.5), 
      mb = list(delay = Inf, w_max = 2), 
      xu = list(early_period = 3), 
      rmst = list(tau = 10), 
      ms = list(ms_time = 10), 
      mc = list(rho = c(0, 0), gamma = c(0, 0.5))
      )
                              
    ans_new
  }

ans <- data.table::rbindlist(ans)

plan("sequential")

ans |> head() |> gt() |> tab_header("Overview Each Simulation results")

ans_non_mc <- ans |>
  filter(Method != "MaxCombo") |>
  group_by(Method, Parameter) %>% 
  summarise(Power = mean(Z > -qnorm(0.025))) |>
  ungroup()

ans_mc <- ans |>
  filter(Method == "MaxCombo") |>
  summarize(Power = mean(`P value` < 0.025), Method = "MaxCombo", Parameter = "FH(0, 0) + FH(0, 0.5)") 

ans_non_mc |>
  union(ans_mc) |>
  gt() |>
  tab_header("Summary from 100 simulations")