Validation of Group Sequential Functions Against gsDesign and rpact

alpha <- 0.025
info_frac_list <- list(
  equally_spaced_3 = c(1/3, 2/3, 1),
  equally_spaced_2 = c(0.5, 1),
  unequal = c(0.2, 0.6, 1),
  early_look = c(0.1, 0.5, 1)
)

tol <- 1e-10
all_pass <- TRUE

for (name in names(info_frac_list)) {
  t <- info_frac_list[[name]]

  # O'Brien-Fleming
  graphicalMCP_of <- spending_of(alpha, t)
  gsd_of <- gsDesign::sfLDOF(alpha, t, param = NULL)$spend
  if (!isTRUE(all.equal(graphicalMCP_of, gsd_of, tolerance = tol))) {
    cat("FAIL: OBF spending,", name, "\n")
    all_pass <- FALSE
  }

  # Pocock
  graphicalMCP_poc <- spending_pocock(alpha, t)
  gsd_poc <- gsDesign::sfLDPocock(alpha, t, param = NULL)$spend
  if (!isTRUE(all.equal(graphicalMCP_poc, gsd_poc, tolerance = tol))) {
    cat("FAIL: Pocock spending,", name, "\n")
    all_pass <- FALSE
  }

  # HSD with gamma = -4
  graphicalMCP_hsd <- spending_hsd(alpha, t, gamma = -4)
  gsd_hsd <- gsDesign::sfHSD(alpha, t, param = -4)$spend
  if (!isTRUE(all.equal(graphicalMCP_hsd, gsd_hsd, tolerance = tol))) {
    cat("FAIL: HSD spending,", name, "\n")
    all_pass <- FALSE
  }

  # Linear (spending_linear is alpha * t by definition)
  graphicalMCP_lin <- spending_linear(alpha, t)
  expected_lin <- alpha * t
  if (!isTRUE(all.equal(graphicalMCP_lin, expected_lin, tolerance = tol))) {
    cat("FAIL: Linear spending,", name, "\n")
    all_pass <- FALSE
  }
}

if (all_pass) cat("All spending function comparisons with gsDesign PASS\n")
#> All spending function comparisons with gsDesign PASS

all_pass <- TRUE

for (name in names(info_frac_list)) {
  t <- info_frac_list[[name]]

  # O'Brien-Fleming
  graphicalMCP_of <- spending_of(alpha, t)
  rpact_of <- getDesignGroupSequential(
    typeOfDesign = "OF",
    informationRates = t,
    alpha = alpha
  )$alphaSpent
  if (!isTRUE(all.equal(graphicalMCP_of, rpact_of, tolerance = 1e-6))) {
    cat("FAIL: OBF spending,", name, "\n")
    all_pass <- FALSE
  }

  # Pocock
  graphicalMCP_poc <- spending_pocock(alpha, t)
  rpact_poc <- getDesignGroupSequential(
    typeOfDesign = "P",
    informationRates = t,
    alpha = alpha
  )$alphaSpent
  if (!isTRUE(all.equal(graphicalMCP_poc, rpact_poc, tolerance = 1e-6))) {
    cat("FAIL: Pocock spending,", name, "\n")
    all_pass <- FALSE
  }
}
#> FAIL: OBF spending, equally_spaced_3 
#> FAIL: Pocock spending, equally_spaced_3 
#> FAIL: OBF spending, equally_spaced_2 
#> FAIL: Pocock spending, equally_spaced_2 
#> FAIL: OBF spending, unequal 
#> FAIL: Pocock spending, unequal 
#> FAIL: OBF spending, early_look 
#> FAIL: Pocock spending, early_look

if (all_pass) cat("All spending function comparisons with rpact PASS\n")

t <- c(1/3, 2/3, 1)

spending_detail <- data.frame(
  `Info Fraction` = t,
  `graphicalMCP (OBF)` = spending_of(alpha, t),
  `gsDesign (sfLDOF)` = gsDesign::sfLDOF(alpha, t, param = NULL)$spend,
  `graphicalMCP (Pocock)` = spending_pocock(alpha, t),
  `gsDesign (sfLDPocock)` = gsDesign::sfLDPocock(alpha, t, param = NULL)$spend,
  check.names = FALSE
)
knitr::kable(spending_detail, digits = 10,
             caption = "Spending function comparison (alpha = 0.025)")

all_pass <- TRUE

test_cases <- list(
  list(alpha = 0.025, t = c(1/3, 2/3, 1), sfu = gsDesign::sfLDOF,
       graphicalMCP_fn = spending_of, label = "OBF, 3 analyses"),
  list(alpha = 0.025, t = c(0.5, 1), sfu = gsDesign::sfLDOF,
       graphicalMCP_fn = spending_of, label = "OBF, 2 analyses"),
  list(alpha = 0.025, t = c(1/3, 2/3, 1), sfu = gsDesign::sfLDPocock,
       graphicalMCP_fn = spending_pocock, label = "Pocock, 3 analyses"),
  list(alpha = 0.01, t = c(0.2, 0.6, 1), sfu = gsDesign::sfLDOF,
       graphicalMCP_fn = spending_of, label = "OBF, unequal, alpha=0.01")
)

for (tc in test_cases) {
  # gsDesign boundaries
  K <- length(tc$t)
  gsd_design <- gsDesign::gsDesign(
    k = K,
    test.type = 1,
    alpha = tc$alpha,
    sfu = tc$sfu,
    timing = tc$t[-K]
  )
  gsd_z <- gsd_design$upper$bound
  gsd_nom <- pnorm(gsd_z, lower.tail = FALSE)

  # graphicalMCP boundaries
  graphicalMCP_bounds <- graphicalMCP:::gs_boundaries(tc$alpha, tc$t, tc$graphicalMCP_fn)

  if (!isTRUE(all.equal(graphicalMCP_bounds$bounds_z, gsd_z, tolerance = 1e-4))) {
    cat("FAIL Z:", tc$label, "\n")
    all_pass <- FALSE
  }
  if (!isTRUE(all.equal(graphicalMCP_bounds$bounds_nominal, gsd_nom, tolerance = 1e-4))) {
    cat("FAIL nominal:", tc$label, "\n")
    all_pass <- FALSE
  }
}
#> FAIL nominal: Pocock, 3 analyses

if (all_pass) cat("All boundary comparisons with gsDesign PASS\n")

all_pass <- TRUE

rpact_cases <- list(
  list(alpha = 0.025, t = c(1/3, 2/3, 1), type = "OF",
       graphicalMCP_fn = spending_of, label = "OBF, 3 analyses"),
  list(alpha = 0.025, t = c(0.5, 1), type = "OF",
       graphicalMCP_fn = spending_of, label = "OBF, 2 analyses"),
  list(alpha = 0.025, t = c(1/3, 2/3, 1), type = "P",
       graphicalMCP_fn = spending_pocock, label = "Pocock, 3 analyses")
)

for (tc in rpact_cases) {
  rpact_design <- getDesignGroupSequential(
    typeOfDesign = tc$type,
    informationRates = tc$t,
    alpha = tc$alpha
  )
  rpact_z <- rpact_design$criticalValues
  rpact_nom <- rpact_design$stageLevels

  graphicalMCP_bounds <- graphicalMCP:::gs_boundaries(tc$alpha, tc$t, tc$graphicalMCP_fn)

  if (!isTRUE(all.equal(graphicalMCP_bounds$bounds_z, rpact_z, tolerance = 1e-4))) {
    cat("FAIL Z:", tc$label, "\n")
    all_pass <- FALSE
  }
  if (!isTRUE(all.equal(graphicalMCP_bounds$bounds_nominal, rpact_nom, tolerance = 1e-4))) {
    cat("FAIL nominal:", tc$label, "\n")
    all_pass <- FALSE
  }
}
#> FAIL Z: OBF, 3 analyses 
#> FAIL nominal: OBF, 3 analyses 
#> FAIL Z: OBF, 2 analyses 
#> FAIL nominal: OBF, 2 analyses 
#> FAIL Z: Pocock, 3 analyses 
#> FAIL nominal: Pocock, 3 analyses

if (all_pass) cat("All boundary comparisons with rpact PASS\n")

t <- c(1/3, 2/3, 1)

graphicalMCP_b <- graphicalMCP:::gs_boundaries(0.025, t, spending_of)
gsd_d <- gsDesign::gsDesign(k = 3, test.type = 1, alpha = 0.025,
                             sfu = gsDesign::sfLDOF, timing = c(1/3, 2/3))
rpact_d <- getDesignGroupSequential(
  typeOfDesign = "OF", informationRates = t, alpha = 0.025
)

boundary_detail <- data.frame(
  Analysis = 1:3,
  `graphicalMCP (Z)` = graphicalMCP_b$bounds_z,
  `gsDesign (Z)` = gsd_d$upper$bound,
  `rpact (Z)` = rpact_d$criticalValues,
  `graphicalMCP (nom. p)` = graphicalMCP_b$bounds_nominal,
  `gsDesign (nom. p)` = pnorm(gsd_d$upper$bound, lower.tail = FALSE),
  `rpact (nom. p)` = rpact_d$stageLevels,
  check.names = FALSE
)
knitr::kable(boundary_detail, digits = 6,
             caption = "Boundary comparison: OBF with 3 equally spaced analyses")

p_obs <- c(0.024, 0.01)
t <- c(0.5, 1)

# Compute repeated p-values at each analysis
graphicalMCP_rep <- c(
  repeated_p(p_obs[1], t[1], spending_of),
  repeated_p(p_obs, t, spending_of)
)

# For each repeated p-value, verify boundary matches observed p
verify <- data.frame(Analysis = integer(), `Observed p` = numeric(),
                     `Repeated p (alpha)` = numeric(),
                     `graphicalMCP boundary` = numeric(),
                     `rpact boundary` = numeric(),
                     check.names = FALSE)

for (k in 1:2) {
  alpha_k <- graphicalMCP_rep[k]
  if (alpha_k >= 1 - 1e-6 || alpha_k <= 1e-6) next

  graphicalMCP_b <- graphicalMCP:::gs_boundaries(alpha_k, t[1:k], spending_of)

  # rpact requires last info rate = 1; only compare at final analysis
  rpact_nom <- NA
  if (k == 2) {
    rpact_d <- suppressMessages(getDesignGroupSequential(
      typeOfDesign = "asOF", informationRates = t, alpha = alpha_k
    ))
    rpact_nom <- rpact_d$stageLevels[k]
  }

  verify <- rbind(verify, data.frame(
    Analysis = k,
    `Observed p` = p_obs[k],
    `Repeated p (alpha)` = alpha_k,
    `graphicalMCP boundary` = graphicalMCP_b$bounds_nominal[k],
    `rpact boundary` = rpact_nom,
    check.names = FALSE
  ))
}

knitr::kable(verify, digits = 6,
             caption = "Boundary verification: observed p should equal boundary at repeated p alpha")

all_pass <- TRUE

for (sc in seq_scenarios) {
  K <- length(sc$p)

  # Compute repeated p-values at each analysis
  rep_p_vals <- numeric(K)
  for (k in 1:K) {
    rep_p_vals[k] <- repeated_p(sc$p[1:k], sc$t[1:k], sc$fn)
  }

  # Compute sequential p-values at each analysis
  seq_p_vals <- numeric(K)
  for (k in 1:K) {
    seq_p_vals[k] <- sequential_p(sc$p[1:k], sc$t[1:k], sc$fn)
  }

  # Sequential p-values should equal cummin of repeated p-values
  cummin_rep <- cummin(rep_p_vals)

  if (!isTRUE(all.equal(seq_p_vals, cummin_rep, tolerance = 1e-4))) {
    cat("FAIL:", sc$label, "\n")
    cat("  seq_p:      ", round(seq_p_vals, 6), "\n")
    cat("  cummin(rep): ", round(cummin_rep, 6), "\n")
    all_pass <- FALSE
  }
}

if (all_pass) cat("All sequential_p == cummin(repeated_p) checks PASS\n")
#> All sequential_p == cummin(repeated_p) checks PASS

p_obs <- c(0.05, 0.02, 0.01)
t <- c(1/3, 2/3, 1)

# graphicalMCP sequential p-values at each analysis
graphicalMCP_seq <- numeric(3)
for (k in 1:3) {
  graphicalMCP_seq[k] <- sequential_p(p_obs[1:k], t[1:k], spending_of)
}

# gsDesign sequential p-value (at the final analysis)
gsd_d <- gsDesign::gsDesign(k = 3, test.type = 1, alpha = 0.025,
                             sfu = gsDesign::sfLDOF, timing = c(1/3, 2/3))
gsd_final <- gsDesign::sequentialPValue(
  gsD = gsd_d, n.I = gsd_d$n.I, Z = qnorm(1 - p_obs)
)

# Repeated p-values for comparison
graphicalMCP_rep <- numeric(3)
for (k in 1:3) {
  graphicalMCP_rep[k] <- repeated_p(p_obs[1:k], t[1:k], spending_of)
}

seq_detail <- data.frame(
  Analysis = 1:3,
  `Observed p` = p_obs,
  `Repeated p` = graphicalMCP_rep,
  `Sequential p` = graphicalMCP_seq,
  `cummin(rep)` = cummin(graphicalMCP_rep),
  check.names = FALSE
)
knitr::kable(seq_detail, digits = 6,
             caption = "Sequential vs. repeated p-values (OBF, 3 analyses)")

cat(sprintf(
  "\ngsDesign sequentialPValue (final): %.6f\ngraphicalMCP sequential_p (final): %.6f\n",
  gsd_final, graphicalMCP_seq[3]
))
#> 
#> gsDesign sequentialPValue (final): 0.010554
#> graphicalMCP sequential_p (final): 0.010554

alpha <- 0.025
info_frac <- c(1/3, 2/3, 1)

wt_comparison <- data.frame()

for (delta in c(0, 0.1, 0.25, 0.4, 0.5)) {
  # graphicalMCP: compute boundaries via spending_wt
  b_graphicalMCP <- gs_boundaries(
    alpha, info_frac,
    function(a, t) spending_wt(a, t, delta = delta)
  )

  # rpact: Wang-Tsiatis boundaries
  rpact_type <- if (delta == 0) {
    "OF"
  } else if (delta == 0.5) {
    "P"
  } else {
    "WT"
  }
  rpact_args <- list(
    sided = 1, alpha = alpha,
    informationRates = info_frac,
    typeOfDesign = rpact_type
  )
  if (rpact_type == "WT") rpact_args$deltaWT <- delta
  rpact_gsd <- do.call(getDesignGroupSequential, rpact_args)

  max_diff <- max(abs(b_graphicalMCP$bounds_z - rpact_gsd$criticalValues))

  wt_comparison <- rbind(wt_comparison, data.frame(
    Delta = delta,
    Type = ifelse(delta == 0, "O'Brien-Fleming",
           ifelse(delta == 0.5, "Pocock",
                  paste0("WT (", delta, ")"))),
    Z1_graphicalMCP = b_graphicalMCP$bounds_z[1],
    Z1_rpact = rpact_gsd$criticalValues[1],
    Z2_graphicalMCP = b_graphicalMCP$bounds_z[2],
    Z2_rpact = rpact_gsd$criticalValues[2],
    Z3_graphicalMCP = b_graphicalMCP$bounds_z[3],
    Z3_rpact = rpact_gsd$criticalValues[3],
    Max.Diff = max_diff,
    check.names = FALSE
  ))
}

knitr::kable(wt_comparison, digits = 6, row.names = FALSE,
             caption = paste("Wang-Tsiatis Z-scale boundaries:",
                             "graphicalMCP vs rpact (alpha = 0.025)"))

shape_check <- data.frame()

for (delta in c(0, 0.25, 0.5)) {
  b <- gs_boundaries(
    alpha, info_frac,
    function(a, t) spending_wt(a, t, delta = delta)
  )
  C_values <- b$bounds_z * info_frac^(0.5 - delta)

  shape_check <- rbind(shape_check, data.frame(
    Delta = delta,
    C1 = C_values[1],
    C2 = C_values[2],
    C3 = C_values[3],
    Max.Variation = max(C_values) - min(C_values)
  ))
}

knitr::kable(shape_check, digits = 6, row.names = FALSE,
             caption = paste("Verification that C = Z * t^(0.5 - delta)",
                             "is constant across analyses"))