Inflated Type 1 Error in MaxCombo

[fb-elong]

Fixed design is using a fixed boundary for the calculation that needs to be calculated adaptively.

https://github.com/Merck/gsDesign2/blob/main/R/fixed_design_maxcombo.R#L105

upper = gs_b, upar = qnorm(1 - alpha),

[fb-elong]

The issue can be fixed in

      upper = gs_spending_combo,
      upar = list(sf = gsDesign::sfLDOF, total_spend = alpha),

Example Code:

devtools::load_all()
library(dplyr)

survival_at_24_months <- 0.35
hr <- log(.35)/log(.25)
control_median <- 12
control_rate <- c(log(2) / control_median, (log(.25) - log(.2)) / 12)

scenarios <- tribble(
  ~Scenario, ~Name,           ~Period, ~duration, ~Survival,
  0,         "Control",       0,       0,         1,
  0,         "Control",       1,       24,        .25,
  0,         "Control",       2,       12,        .2,
  1,         "PH",            0,       0,         1,
  1,         "PH",            1,       24,        .35,
  1,         "PH",            2,       12,        .2^hr,
  2,         "3-month delay", 0,       0,         1,
  2,         "3-month delay", 1,       3,         exp(-3 * control_rate[1]),
  2,         "3-month delay", 2,       21,        .35,
  2,         "3-month delay", 3,       12,        .2^hr,
  3,         "6 month delay", 0,       0,         1,
  3,         "6-month delay", 1,       6,         exp(-6 * control_rate[1]),
  3,         "6-month delay", 2,       18,        .35,
  3,         "6-month delay", 3,       12,        .2^hr,
  4,         "Crossing",      0,       0,         1,
  4,         "Crossing",      1,       3,         exp(-3 * control_rate[1] * 1.3),
  4,         "Crossing",      2,       21,        .35,
  4,         "Crossing",      3,       12,        .2^hr,
  5,         "Weak null",     0,       0,         1,
  5,         "Weak null",     1,       24,        .25,
  5,         "Weak null",     2,       12,        .2,
  6,         "Strong null",   0,       0,         1,         
  6,         "Strong null",   1,       3,         exp(-3 * control_rate[1] * 1.5),
  6,         "Strong null",   2,       3,         exp(-6 * control_rate[1]),
  6,         "Strong null",   3,       18,        .25,
  6,         "Strong null",   4,       12,        .2,
)
scenarios <-
  scenarios |> group_by(Scenario) |> 
  mutate(Month = cumsum(duration)) 

fr <- 
  scenarios |>
  group_by(Scenario) |>
  #  filter(Scenario == 2) |>
  mutate(x_rate = -(log(Survival) - log(lag(Survival, default = 1))) / 
           duration, 
         rate = ifelse(Month > 24, control_rate[2], control_rate[1]),
         hr = x_rate / rate) |>
  select(-x_rate) |>
  filter(Period > 0, Scenario > 0) |> ungroup()
fr <- fr |> mutate(fail_rate = rate, dropout_rate =0.001, stratum = "All")


mwlr <- fixed_design_mb(
  tau = 12,
  enroll_rate = define_enroll_rate(duration = 12, rate = 1),
  fail_rate = fr |> filter(Scenario == 2),
  alpha = 0.025, power = .85, ratio = 1,
  study_duration = 36
) |> to_integer()


er <- mwlr$enroll_rate
ss <- NULL

scen <- 5
# MaxCombo
maxcombo <- fixed_design_maxcombo(
  rho = c(0, 0), gamma = c(0, 0.5), tau = rep(-1, 2),
  enroll_rate = er,
  fail_rate = fr |> filter(Scenario == scen),
  alpha = 0.025, power = NULL, ratio = 1,
  study_duration = 36
)
ss <- rbind(ss, maxcombo$analysis |> select(-c("bound", "alpha")) |>
              mutate(Scenario = scen, design="MaxCombo"))
ss

[LittleBeannie]

Hi @fb-elong , will you work on this issue?

[elong0527]

sure, a PR has been submitted. Please consider to rerun simulation and cross check the results.

[LittleBeannie]

Hi @elong0527 , thanks for getting this issue fixed. After running the following codes, the asymptotic type I error is well controlled at 0.025. I further explored 1 million simulations and the simulated type I error is also around 0.025.

The only slight difference is the correlation (yellow-highlighted in the screenshot below), but I guess it should be fine. So, I will get this PR merged. Thank you so much!

Looping @keaven for information.

library(dplyr)
library(tibble)
library(gsDesign2)
library(simtrial)
library(doFuture)
library(foreach)

# ------------------------------- #
#   define enroll and fail rate   #
# ------------------------------- #
# sample size from the MWLR test
enroll_rate <- define_enroll_rate(duration = 12, rate = 698/12)
# week null
fail_rate <- define_fail_rate(duration = c(24, 12),
                              fail_rate = c(log(2) / 12, (log(.25) - log(.2)) / 12),
                              dropout_rate = 0.001,
                              hr = 1)

# ------------------------------- #
#   design by maxcombo            #
# ------------------------------- #
# asymptotic design
x <- fixed_design_maxcombo(
  rho = c(0, 0), gamma = c(0, 0.5), tau = rep(-1, 2),
  enroll_rate = define_enroll_rate(duration = 12, rate = 698/12),
  fail_rate = fail_rate,
  alpha = 0.025, power = NULL, ratio = 1,
  study_duration = 36)

x |> summary() |> as_gt()

# calculate the asymptotic correlation of maxcombo
gsDesign2:::gs_utility_combo(
  enroll_rate = define_enroll_rate(duration = 12, rate = 698/12),
  fail_rate = fail_rate,
  fh_test = data.frame(rho = c(0, 0), gamma = c(0, 0.5), tau = c(-1, -1), test = 1:2,
                       analysis = rep(1, 2), analysis_time = c(36, 36)),
  ratio = 1,
  algorithm = GenzBretz(maxpts = 1e5, abseps = 1e-8))$corr[1, 2]

# ------------------------------- #
#   run simulations               #
# ------------------------------- #
# build a utility function for a single simulation
sim_fn <- function(enroll_rate, fail_rate){
  temp <- to_sim_pw_surv(fail_rate)
  
  # Generate a dataset
  dat <- sim_pw_surv(n = 698, enroll_rate = enroll_rate,
                     fail_rate = temp$fail_rate, dropout_rate = temp$dropout_rate)
  
  analysis_data <- cut_data_by_date(dat, 36)
  # Do the maxcombo test
  maxcombopz <- (analysis_data |> maxcombo(rho = c(0,0), gamma = c(0,0.5), return_corr = TRUE))
  maxcombop <- maxcombopz$p_value
  maxcombocorr <- maxcombopz$corr[1, 2]
  
  return(data.frame(p_value = maxcombop, 
                    corr = maxcombocorr))
}

# run parallel simulation 
plan(multisession(workers = 24))
set.seed(3219)
tictoc::tic()
results <- foreach(sim = 1:1e6, .combine = 'rbind', .options.future = list(seed = TRUE)) %dofuture%{
  sim_fn(enroll_rate = enroll_rate, fail_rate = fail_rate) %>% mutate(sim_id = sim)
}
tictoc::toc()

results |>
  summarize(power = mean(p_value <= .025), correlation = mean(corr))