## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup, message = FALSE, warning = FALSE----------------------------------
library(graphicalMCP)
library(gsDesign)
library(rpact)

## ----spending-gsdesign--------------------------------------------------------
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")

## ----spending-rpact-----------------------------------------------------------
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
  }
}

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

## ----spending-detail----------------------------------------------------------
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)")

## ----boundaries-gsdesign------------------------------------------------------
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
  }
}

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

## ----boundaries-rpact---------------------------------------------------------
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
  }
}

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

## ----boundaries-detail--------------------------------------------------------
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")

## ----corr-matrix--------------------------------------------------------------
all_pass <- TRUE

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

  # graphicalMCP
  graphicalMCP_corr <- graphicalMCP:::gs_corr(t)

  # Manual formula
  manual_corr <- outer(t, t, function(ti, tj) sqrt(pmin(ti, tj) / pmax(ti, tj)))

  if (!isTRUE(all.equal(graphicalMCP_corr, manual_corr, tolerance = 1e-12))) {
    cat("FAIL:", name, "\n")
    all_pass <- FALSE
  }

  # Verify properties
  stopifnot(
    all(diag(graphicalMCP_corr) == 1),
    isSymmetric(graphicalMCP_corr)
  )
}

if (all_pass) cat("All correlation matrix comparisons PASS\n")

## ----corr-detail--------------------------------------------------------------
t <- c(1/3, 2/3, 1)
knitr::kable(graphicalMCP:::gs_corr(t), digits = 6,
             caption = "Correlation matrix for info fractions (1/3, 2/3, 1)")

## ----repeated-p-boundary------------------------------------------------------
scenarios <- list(
  list(
    p = c(0.024, 0.01),
    t = c(0.5, 1),
    type = "asOF",
    fn = spending_of,
    label = "OBF, 2 analyses"
  ),
  list(
    p = c(0.05, 0.02, 0.01),
    t = c(1/3, 2/3, 1),
    type = "asOF",
    fn = spending_of,
    label = "OBF, 3 analyses"
  ),
  list(
    p = c(0.024, 0.01),
    t = c(0.5, 1),
    type = "asP",
    fn = spending_pocock,
    label = "Pocock, 2 analyses"
  ),
  list(
    p = c(0.1, 0.05, 0.01),
    t = c(0.2, 0.6, 1),
    type = "asOF",
    fn = spending_of,
    label = "OBF, unequal spacing"
  )
)

all_pass <- TRUE

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

  for (k in 1:K) {
    # Compute repeated p-value at analysis k
    alpha_rep <- repeated_p(
      p = sc$p[1:k],
      info_frac = sc$t[1:k],
      spending_fn = sc$fn
    )

    # Skip if repeated p-value is at the boundary (1 or near 0)
    if (alpha_rep >= 1 - 1e-6 || alpha_rep <= 1e-6) next

    # graphicalMCP boundary at this alpha
    graphicalMCP_bounds <- graphicalMCP:::gs_boundaries(
      alpha_rep, sc$t[1:k], sc$fn
    )
    graphicalMCP_nom_k <- graphicalMCP_bounds$bounds_nominal[k]

    # The observed p-value should equal the boundary at analysis k
    if (!isTRUE(all.equal(sc$p[k], graphicalMCP_nom_k, tolerance = 1e-4))) {
      cat("FAIL (graphicalMCP boundary):", sc$label, "analysis", k, "\n")
      cat("  observed p:", sc$p[k], "  boundary:", graphicalMCP_nom_k, "\n")
      all_pass <- FALSE
    }

    # rpact boundary at this alpha (only when info_frac ends at 1,
    # since rpact requires the last information rate to be 1)
    if (k == K) {
      rpact_design <- suppressMessages(getDesignGroupSequential(
        typeOfDesign = sc$type,
        informationRates = sc$t,
        alpha = alpha_rep
      ))
      rpact_nom_k <- rpact_design$stageLevels[k]
      if (!isTRUE(all.equal(sc$p[k], rpact_nom_k, tolerance = 1e-3))) {
        cat("FAIL (rpact boundary):", sc$label, "analysis", k, "\n")
        cat("  observed p:", sc$p[k], "  boundary:", rpact_nom_k, "\n")
        all_pass <- FALSE
      }
    }
  }
}

if (all_pass) cat("All repeated p-value boundary checks PASS\n")

## ----repeated-p-detail--------------------------------------------------------
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")

## ----sequential-p-gsdesign----------------------------------------------------
all_pass <- TRUE

seq_scenarios <- list(
  list(
    p = c(0.024, 0.01),
    t = c(0.5, 1),
    sfu = gsDesign::sfLDOF,
    fn = spending_of,
    label = "OBF, 2 analyses"
  ),
  list(
    p = c(0.05, 0.02, 0.01),
    t = c(1/3, 2/3, 1),
    sfu = gsDesign::sfLDOF,
    fn = spending_of,
    label = "OBF, 3 analyses"
  ),
  list(
    p = c(0.024, 0.01),
    t = c(0.5, 1),
    sfu = gsDesign::sfLDPocock,
    fn = spending_pocock,
    label = "Pocock, 2 analyses"
  ),
  list(
    p = c(0.1, 0.05, 0.01),
    t = c(0.2, 0.6, 1),
    sfu = gsDesign::sfLDOF,
    fn = spending_of,
    label = "OBF, unequal spacing"
  )
)

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

  # graphicalMCP
  graphicalMCP_seq <- sequential_p(sc$p, sc$t, sc$fn)

  # gsDesign
  gsd_design <- gsDesign::gsDesign(
    k = K,
    test.type = 1,
    alpha = 0.025,
    sfu = sc$sfu,
    timing = sc$t[-K]
  )
  gsd_seq <- gsDesign::sequentialPValue(
    gsD = gsd_design,
    n.I = gsd_design$n.I,
    Z = z_obs
  )

  if (!isTRUE(all.equal(graphicalMCP_seq, gsd_seq, tolerance = 1e-4))) {
    cat("FAIL:", sc$label, "\n")
    cat("  graphicalMCP:", graphicalMCP_seq, "\n")
    cat("  gsDesign:    ", gsd_seq, "\n")
    all_pass <- FALSE
  }
}

if (all_pass) cat("All sequential p-value comparisons with gsDesign PASS\n")

## ----seq-equals-cummin-rep----------------------------------------------------
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")

## ----sequential-p-detail------------------------------------------------------
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]
))

## ----wt-rpact-----------------------------------------------------------------
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)"))

## ----wt-shape-----------------------------------------------------------------
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"))

## ----summary-table------------------------------------------------------------
summary_table <- data.frame(
  Function = c(
    "spending_of / spending_pocock / spending_hsd / spending_linear",
    "gs_boundaries (Z and nominal p)",
    "gs_corr",
    "repeated_p",
    "sequential_p",
    "sequential_p == cummin(repeated_p)",
    "spending_wt (Wang-Tsiatis)"
  ),
  `Compared Against` = c(
    "gsDesign (sfLDOF, sfLDPocock, sfHSD, sfLinear), rpact (alphaSpent)",
    "gsDesign (gsDesign), rpact (getDesignGroupSequential)",
    "Analytical formula",
    "Boundary inversion: gsDesign + rpact boundaries",
    "gsDesign (sequentialPValue)",
    "Internal consistency",
    "rpact (getDesignGroupSequential with typeOfDesign OF/P/WT)"
  ),
  check.names = FALSE
)
knitr::kable(summary_table, caption = "Summary of validation comparisons")