Merck / simtrial

Clinical trial simulation for time-to-event endpoints
https://merck.github.io/simtrial/
GNU General Public License v3.0
16 stars 8 forks source link

Add the info to the output of `sim_gs_n()` #259

Closed LittleBeannie closed 4 weeks ago

LittleBeannie commented 2 months ago

The PR will fix the 2 issues below.

  1. The weighting function (such as fh_weight, mb_weight and early_zero_weight) now all gives a column of "weight".
  2. The info and info0 are added as output of sim_gs_n.
LittleBeannie commented 1 month ago

To test the calculation of the simulated info and info0, I compared the numbers with the asymptotic results (gsDesign2 branch https://github.com/Merck/gsDesign2/tree/441-info-frac-driven-design-does-not-work-for-gs_design_wlr).

Comparision 1: Regular logrank test image

# asymptotic info
x <- gs_design_wlr(enroll_rate = enroll_rate,
                   fail_rate = fail_rate,
                   ratio = 1,
                   alpha = 0.025, beta = 0.2,
                   weight = function(x, arm0, arm1) {wlr_weight_fh(x, arm0, arm1, rho = 0, gamma = 0)},
                   upper = gs_spending_bound,
                   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
                   lower = gs_b,
                   lpar = rep(-Inf, 1),
                   analysis_time = 36) |> to_integer()

y <- gs_design_ahr(enroll_rate = enroll_rate,
                   fail_rate = fail_rate,
                   ratio = 1,
                   alpha = 0.025, beta = 0.2,
                   upper = gs_spending_bound,
                   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
                   lower = gsDesign2::gs_b,
                   lpar = rep(-Inf, 1),
                   analysis_time = 36) |> to_integer()

# simulated info
registerDoFuture()
registerDoRNG()
plan("multisession", workers = 24)
set.seed(2024)
sim_res <- foreach(
  sim_id = seq_len(1e4),
  .combine = "rbind",
  .errorhandling = "pass"
) %dorng% {
  do.call(cbind,
          sim_pw_surv(n = round(x$analysis$n, 0),
                      stratum = data.frame(stratum = "All", p = 1),
                      block = c(rep("control", 2), rep("experimental", 2)),
                      enroll_rate = x$enroll_rate,
                      fail_rate = to_sim_pw_surv(fail_rate)$fail_rate,
                      dropout_rate = to_sim_pw_surv(fail_rate)$dropout_rate) |>
            cut_data_by_event(round(x$analysis$event, 0)) |>
            wlr(weight = fh(rho = 0, gamma = 0))
          ) |>
    as.data.frame() |>
    dplyr::mutate(sim_id = sim_id)
}
plan("sequential")

tibble(From = c("gs_design_ahr", "gs_design_wlr", "simulation"),
       Info = c(y$analysis$info, x$analysis$info, mean(as.numeric(sim_res$info))),
       Info0 = c(y$analysis$info0, x$analysis$info0, mean(as.numeric(sim_res$info0)))
       ) %>% 
  gt() %>% 
  tab_header("Statistcial information of the logrank test")
LittleBeannie commented 1 month ago

Comparision 2: WLR-FH(0, 0.5) image

# asymptotic info
x <- gs_design_wlr(enroll_rate = enroll_rate,
                   fail_rate = fail_rate,
                   ratio = 1, alpha = 0.025, beta = 0.2,
                   weight = function(x, arm0, arm1) {
                     wlr_weight_fh(x, arm0, arm1, rho = 0, gamma = 0.5)
                     },
                   upper = gs_spending_bound,
                   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
                   lower = gsDesign2::gs_b,
                   lpar = rep(-Inf, 1),
                   analysis_time = 36) |> to_integer()

# simulated info
registerDoFuture()
registerDoRNG()
plan("multisession", workers = 24)
set.seed(2024)
sim_res <- foreach(
  sim_id = seq_len(1e4),
  .combine = "rbind",
  .errorhandling = "pass"
) %dorng% {
  do.call(cbind,
          sim_pw_surv(n = round(x$analysis$n, 0),
                      stratum = data.frame(stratum = "All", p = 1),
                      block = c(rep("control", 2), rep("experimental", 2)),
                      enroll_rate = x$enroll_rate,
                      fail_rate = to_sim_pw_surv(fail_rate)$fail_rate,
                      dropout_rate = to_sim_pw_surv(fail_rate)$dropout_rate) |>
            cut_data_by_event(round(x$analysis$event, 0)) |>
            wlr(weight = fh(rho = 0, gamma = 0.5))
    ) |>
    as.data.frame() |>
    mutate(sim_id = sim_id)
  }
plan("sequential")

tibble(From = c("gs_design_wlr", "simulation"),
       Info = c(x$analysis$info, mean(as.numeric(sim_res$info))),
       Info0 = c(x$analysis$info0, mean(as.numeric(sim_res$info0)))
       ) %>% 
  gt() %>% 
  tab_header("Statistcial information of the WLR-FH(0, 0.5) test")
LittleBeannie commented 1 month ago

Comparision 3: WLR-MB(12, 2) image

# asymptotic info
x <- gs_design_wlr(enroll_rate = enroll_rate,
                   fail_rate = fail_rate,
                   ratio = 1, alpha = 0.025, beta = 0.2,
                   weight = function(x, arm0, arm1) {
                     gsDesign2::wlr_weight_mb(x, arm0, arm1, tau = 12, wmax = 2)
                     },
                   upper = gsDesign2::gs_spending_bound,
                   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
                   lower = gsDesign2::gs_b,
                   lpar = rep(-Inf, 1),
                   analysis_time = 36) |> gsDesign2::to_integer()

# simulated info
registerDoFuture()
registerDoRNG()
plan("multisession", workers = 24)
set.seed(2024)
sim_res <- foreach(
  sim_id = seq_len(1e4),
  .combine = "rbind",
  .errorhandling = "pass"
) %dorng% {
  do.call(cbind,
          sim_pw_surv(n = round(x$analysis$n, 0),
                      stratum = data.frame(stratum = "All", p = 1),
                      block = c(rep("control", 2), rep("experimental", 2)),
                      enroll_rate = x$enroll_rate,
                      fail_rate = to_sim_pw_surv(fail_rate)$fail_rate,
                      dropout_rate = to_sim_pw_surv(fail_rate)$dropout_rate) |>
            cut_data_by_event(round(x$analysis$event, 0)) |>
            wlr(weight = mb(delay = 12, w_max = 2))
          ) |>
    as.data.frame() |>
    mutate(sim_id = sim_id)
}
plan("sequential")

tibble(From = c("gs_design_wlr", "simulation"),
       Info = c(x$analysis$info, mean(as.numeric(sim_res$info))),
       Info0 = c(x$analysis$info0, mean(as.numeric(sim_res$info0)))
       ) %>% 
  gt() %>% 
  tab_header("Statistcial information of the WLR-MB(12, 2) test")
LittleBeannie commented 1 month ago

Comparision 4: WLR-Xu(6) image

# asymptotic info
x <- gs_design_wlr(enroll_rate = enroll_rate,
                   fail_rate = fail_rate,
                   ratio = 1, alpha = 0.025, beta = 0.2,
                   weight = function(x, arm0, arm1) {x > 6},
                   upper = gs_spending_bound,
                   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
                   lower = gs_b,
                   lpar = rep(-Inf, 1),
                   analysis_time = 36) |> gsDesign2::to_integer()
# simulated info
registerDoFuture()
registerDoRNG()
plan("multisession", workers = 24)
set.seed(2024)
sim_res <- foreach(
  sim_id = seq_len(1e4),
  .combine = "rbind",
  .errorhandling = "pass"
) %dorng% {
  do.call(cbind,
          sim_pw_surv(n = round(x$analysis$n, 0),
                      stratum = data.frame(stratum = "All", p = 1),
                      block = c(rep("control", 2), rep("experimental", 2)),
                      enroll_rate = x$enroll_rate,
                      fail_rate = to_sim_pw_surv(fail_rate)$fail_rate,
                      dropout_rate = to_sim_pw_surv(fail_rate)$dropout_rate) |>
            cut_data_by_event(round(x$analysis$event, 0)) |>
            wlr(weight = early_zero(6))
          ) |>
    as.data.frame() |>
    mutate(sim_id = sim_id)
}
plan("sequential")

tibble(From = c("gs_design_wlr", "simulation"),
       Info = c(x$analysis$info, mean(as.numeric(sim_res$info))),
       Info0 = c(x$analysis$info0, mean(as.numeric(sim_res$info0)))
       ) %>% 
  gt() %>% 
  tab_header("Statistcial information of the WLR-Xu(6) test")