graphicalMCP-gsDesign2

GSD with Graphical Multiplicity Control

Guide users through group sequential design with graphical multiplicity control using the graphicalMCP and gsDesign2 R packages. Use this skill whenever the user asks about: group sequential designs with multiple hypotheses, graphical multiplicity testing, sequential p-values with gsDesign2, combining graphicalMCP with gsDesign2, clinical trial designs with multiple endpoints and populations, Maurer-Bretz procedures, alpha-spending with multiplicity graphs, or adapting the gMCPLite vignette template. Also trigger when users mention spending time, information fraction, or sequential p-values in the context of group sequential or graphical testing.

Group Sequential Design with graphicalMCP + gsDesign2

This skill helps users design and analyze clinical trials that combine graphical multiplicity control (graphicalMCP) with group sequential designs (gsDesign2). The workflow follows the Maurer-Bretz (2013) framework.

When to use this skill

  • Setting up a multiplicity graph for multiple hypotheses (endpoints x populations)
  • Designing group sequential bounds for each hypothesis using gsDesign2
  • Computing sequential p-values from observed data
  • Testing hypotheses using graphicalMCP with sequential p-values
  • Verifying rejection decisions with updated group sequential bounds

Required packages

library(dplyr)
library(tibble)
library(gsDesign)
library(gsDesign2)
library(graphicalMCP)

gsDesign2 must export sequential_pval(). Install from GitHub if needed:

remotes::install_github("Merck/gsDesign2")

Workflow overview

The workflow has 4 phases:

Phase 1: Design specification

  1. Define hypotheses — typically endpoints (OS, PFS, ORR) crossed with populations (subgroup, overall).
  2. Build the multiplicity graph — assign initial alpha weights and transition matrix using graphicalMCP::graph_create().
  3. Choose a sample-size-driving hypothesis (typically OS in the subgroup). Design it with gsDesign2::gs_design_ahr() targeting the desired power (e.g., 90%). Use info_frac = NULL and specify analysis_time as calendar months; gsDesign2 derives the information fraction from the enrollment/failure rate assumptions and analysis timing.
  4. Derive enrollment rates from the driving hypothesis. The subgroup enrollment rate comes directly from the design output. The complement enrollment rate is scaled by prevalence: rate_complement = rate_sub * (1 - prevalence) / prevalence. Build stratified enrollment for overall population designs using define_enroll_rate() with stratum columns.
  5. Compute power for remaining hypotheses:
    • Time-to-event hypotheses (OS, PFS) in the subgroup or overall: use gsDesign2::gs_power_ahr() with the derived enrollment rates. Pass event = NULL so analysis_time drives the design. For overall population, use stratified fail_rate with different HRs per stratum.
    • Binary endpoints (ORR): use gsDesign2::fixed_design_rd() with sample sizes derived from the driving hypothesis. These get NULL in the design list.
  6. Specify analysis timing rules — document when each analysis is triggered (minimum follow-up after FPE, event count thresholds, maximum extensions).
  7. Store designs in an ordered list matching the hypothesis order in the graph. Use NULL for non-GSD hypotheses.

Phase 2: Results entry

  1. Record event counts at each analysis for each hypothesis.
  2. Record nominal one-sided p-values for each analysis of each hypothesis.
  3. Compute spending times — typically events / max(events) using the subgroup information fraction. The spending time must reach 1 at the final analysis of each hypothesis.

Phase 3: Hypothesis testing

  1. Compute sequential p-values using gsDesign2::sequential_pval() for each group sequential hypothesis. For non-GSD hypotheses, the nominal p-value is the sequential p-value.
  2. Test with graphicalMCP using graphicalMCP::graph_test_shortcut() with the sequential p-values and total FWER alpha.

Phase 4: Verification

  1. Extract the graph update sequence using graphicalMCP::graph_update() to see the multiplicity graph at each rejection step.
  2. Update group sequential bounds at the maximum alpha allocated to each hypothesis using gsDesign2::gs_update_ahr().
  3. Compare nominal p-values to updated bounds to confirm rejection decisions.

Key code patterns

For detailed code templates covering each phase, read references/code_patterns.md.

Important design considerations

  • H1 drives sample size: One hypothesis (typically OS in the subgroup) is designed with gs_design_ahr() to determine enrollment rates and sample sizes. All other hypotheses derive their enrollment from H1 using gs_power_ahr() (time-to-event) or fixed_design_rd() (binary).
  • Stratified overall designs: Overall population designs use stratified define_enroll_rate() and define_fail_rate() with separate stratum rows (e.g., “BM+” and “BM-”) and stratum-specific HRs or rates.
  • gs_power_ahr() API: Does not accept info_frac. Use event = NULL with analysis_time to let timing drive the design. If event is not set to NULL, the default c(30, 40, 50) may cause length mismatches for designs with fewer analyses.
  • fixed_design_rd() output: Returns a fixed_design object. Wrap with summary() before piping to gt() or kable().
  • Zero initial alpha: If a hypothesis starts with alpha=0 (receives alpha only through reallocation), use another hypothesis’s alpha for the bounds structure in the design. The actual testing uses the reallocated alpha from the graph.
  • Spending time vs information fraction: Spending time determines how alpha is allocated across analyses. Information fraction drives the correlation structure. Both are needed for bound computation.
  • Non-binding futility: Use binding = FALSE so efficacy bounds are computed ignoring the futility bound, preserving Type I error control even if the trial continues past a futility crossing. Theoretical basis: Liu & Anderson (2008) Theorem 1.
  • Hung-Wang-O’Neill warning: Naively testing a secondary endpoint at level α whenever the primary is significant does NOT control FWER in a group sequential trial (Hung, Wang & O’Neill, 2007). Must use sequential p-values with a proper closed testing or graphical procedure.
  • Well-ordered spending functions: For the Maurer-Bretz graphical procedure to be consonant (sequentially rejective), the nominal significance levels α*_t(γ) must be non-decreasing in γ. Qualified families: power (αt^ρ, all ρ > 0), Pocock-type, OBF-type (for γ < 0.318, covering all practical significance levels). See Maurer & Bretz (2013).
  • Time travel: If OS hypotheses are rejected and alpha passes to previously-tested PFS hypotheses, those PFS tests can be re-evaluated with updated bounds. This controls Type I error per Liu & Anderson (2008).
  • One-sided testing: Maurer-Bretz designs assume one-sided testing or non-binding futility bounds.
  • Alpha spending function: Lan-DeMets spending approximating O’Brien-Fleming (sfLDOF) is a common default.

Key references

  • Maurer W, Bretz F. Multiple testing in group sequential trials using graphical approaches. Stat Biopharm Res 2013; 5:311–320.
  • Liu Q, Anderson KM. On adaptive extensions of group sequential trials for clinical investigations. JASA 2008; 103:1621–1630.
  • Hung HMJ, Wang SJ, O’Neill R. Statistical considerations for testing multiple endpoints in group sequential or adaptive clinical trials. J Biopharm Stat 2007; 17:1201–1210.
  • Spending time rules: Use pmin(planned_IF, actual_IF) at interim analyses to protect against over-spending. Align spending time across populations (H2 uses H1’s spending time). At final analysis, spending time = 1.
  • Stratified ORR: Use gs_power_rd() with weight = "invar" for stratified overall population ORR, not fixed_design_rd().
  • Repeated p-values: In verification, compute both sequential p-values (cumulative evidence) and repeated p-values (single-analysis evidence) to understand each analysis’s contribution.

Code Patterns

Code Patterns for graphicalMCP + gsDesign2

Table of Contents

  1. Multiplicity graph setup
  2. Sample-size-driving design (H1)
  3. Deriving enrollment from H1
  4. Power for remaining hypotheses
  5. Results entry template
  6. Sequential p-value computation
  7. Hypothesis testing with graphicalMCP
  8. Verification with updated bounds
  9. Spending time rules
  10. Stratified ORR power with gs_power_rd()
  11. Displaying design bounds
  12. Repeated p-values for verification
  13. Simulation with illness-death model
  14. P-value computation, theoretical curves, and KM plots
  15. WPGSD analysis

Multiplicity graph setup

Define hypotheses, alpha allocation, and transition matrix:

# Hypothesis names (endpoints x populations)
nameHypotheses <- c(
  "H1: OS\n Subgroup",
  "H2: OS\n All subjects",
  "H3: PFS\n Subgroup",
  "H4: PFS\n All subjects",
  "H5: ORR\n Subgroup",
  "H6: ORR\n All subjects"
)
nHypotheses <- length(nameHypotheses)

# Transition matrix (row i -> col j reallocation weights)
m <- matrix(c(
  0, 1, 0, 0, 0, 0,
  0, 0, .5, .5, 0, 0,
  0, 0, 0, 1, 0, 0,
  0, 0, 0, 0, .5, .5,
  0, 0, 0, 0, 0, 1,
  .5, .5, 0, 0, 0, 0
), nrow = 6, byrow = TRUE)

# Initial alpha allocation (one-sided)
alphaHypotheses <- c(.01, .01, .004, 0.000, 0.0005, .0005)
fwer <- sum(alphaHypotheses)

# Create graph (weights must sum to 1, so divide by fwer)
g0 <- graphicalMCP::graph_create(
  hypotheses = stats::setNames(alphaHypotheses / fwer, nameHypotheses),
  transitions = m,
  hyp_names = nameHypotheses
)

# Plot with alpha levels on vertices (using α character)
alphaLabels <- sprintf("\u03b1 = %s", format(alphaHypotheses, scientific = FALSE))
vertexLabels <- paste(nameHypotheses, alphaLabels, sep = "\n")

plot(g0,
  layout = layout6,
  vertex.size = 45,
  vertex.label = vertexLabels,
  vertex.label.cex = 0.9,
  vertex.color = vertex_colors,
  margin = 0.25
)

Sample-size-driving design (H1)

One hypothesis drives sample size — typically OS in the subgroup, designed with gs_design_ahr() targeting the desired power. Use info_frac = NULL and specify analysis_time as calendar months from enrollment start. gsDesign2 derives the information fraction from the enrollment/failure rate assumptions and analysis timing.

fail_rate_os_sub <- gsDesign2::define_fail_rate(
  duration = Inf,
  fail_rate = log(2) / osmedian,
  hr = 0.65,
  dropout_rate = 0.001
)

# Initial enrollment rate with ramp-up (gs_design_ahr will scale)
enroll_rate_sub_init <- gsDesign2::define_enroll_rate(
  duration = c(2, 2, 2, 8),
  rate = c(0.25, 0.50, 0.75, 1.00)  # Relative rates
)

ossub <- gsDesign2::gs_design_ahr(
  enroll_rate = enroll_rate_sub_init,
  fail_rate = fail_rate_os_sub,
  alpha = alphaHypotheses[1],
  beta = 0.1,
  binding = FALSE,
  analysis_time = c(20, 28, 38),   # Calendar months from enrollment start
  info_frac = NULL,                  # Derived from analysis_time + enroll/fail rates
  info_scale = "h0_info",
  upper = gsDesign2::gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = alphaHypotheses[1]),
  test_lower = FALSE
)

Analysis timing rules

Analysis timing should be pre-specified with rules relative to final patient enrolled (FPE):

Analysis Timing rule (after FPE) Max extension Endpoints assessed
IA1 6 months after FPE None ORR (final), PFS and OS (interim)
IA2 14 months after FPE AND targeted PFS events in subgroup +3 months PFS (final), OS (interim)
Final 24 months after FPE AND targeted OS events in subgroup +6 months OS (final)

These rules ensure adequate follow-up for each endpoint while providing flexibility for event-driven timing. The analysis_time values in gs_design_ahr() should reflect enrollment duration + follow-up (e.g., 14 months enrollment + 6 months = 20 months).

Deriving enrollment from H1

The driving hypothesis determines enrollment rates and sample sizes for all other designs.

# Extract subgroup enrollment rate from H1 design
enroll_rate_sub <- ossub$enroll_rate
n_sub <- sum(enroll_rate_sub$rate * enroll_rate_sub$duration)
n_complement <- n_sub * (1 - prevalence) / prevalence
n_total <- n_sub + n_complement

# Build stratified enrollment for overall population designs
enroll_rate_overall <- gsDesign2::define_enroll_rate(
  stratum = rep(c("BM+", "BM-"), each = nrow(enroll_rate_sub)),
  duration = rep(enroll_rate_sub$duration, 2),
  rate = c(
    enroll_rate_sub$rate,                                    # Subgroup rates from H1
    enroll_rate_sub$rate * (1 - prevalence) / prevalence     # Complement rates
  )
)

Power for remaining hypotheses

Time-to-event in subgroup — gs_power_ahr()

pfssub <- gsDesign2::gs_power_ahr(
  enroll_rate = enroll_rate_sub,
  fail_rate = fail_rate_pfs_sub,
  ratio = 1,
  event = NULL,                    # IMPORTANT: must be NULL to use analysis_time
  analysis_time = c(20, 28),       # 2 analyses for PFS
  info_scale = "h0_info",
  upper = gsDesign2::gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = alphaHypotheses[3]),
  test_lower = FALSE,
  binding = FALSE
)

Time-to-event in overall population — stratified gs_power_ahr()

Use stratified define_fail_rate() with different HRs per stratum:

fail_rate_os_overall <- gsDesign2::define_fail_rate(
  stratum = c("BM+", "BM-"),
  duration = c(Inf, Inf),
  fail_rate = log(2) / osmedian,
  hr = c(0.65, 0.85),
  dropout_rate = 0.001
)

os <- gsDesign2::gs_power_ahr(
  enroll_rate = enroll_rate_overall,
  fail_rate = fail_rate_os_overall,
  ratio = 1,
  event = NULL,
  analysis_time = c(20, 28, 38),
  info_scale = "h0_info",
  upper = gsDesign2::gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = alphaHypotheses[2]),
  test_lower = FALSE,
  binding = FALSE
)

Zero initial alpha workaround

If a hypothesis starts with alpha=0 (e.g., H4: PFS Overall), the spending function will error. Use another hypothesis’s alpha for the bounds structure:

# H4 starts with alpha=0; use H3's alpha for bounds structure
pfs <- gsDesign2::gs_power_ahr(
  enroll_rate = enroll_rate_overall,
  fail_rate = fail_rate_pfs_overall,
  ratio = 1,
  event = NULL,
  analysis_time = c(20, 28),
  info_scale = "h0_info",
  upper = gsDesign2::gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = alphaHypotheses[3]),  # H3's alpha
  test_lower = FALSE,
  binding = FALSE
)

Binary endpoint — fixed_design_rd()

For rate difference tests (e.g., ORR), use fixed_design_rd() with sample sizes from H1. Wrap with summary() before piping to gt().

# Subgroup
orr_sub <- gsDesign2::fixed_design_rd(
  alpha = alphaHypotheses[5],
  power = NULL,
  p_c = 0.30,
  p_e = 0.45,
  rd0 = 0,
  n = ceiling(n_sub) * 2  # Total N (both arms)
)
summary(orr_sub) %>%
  gt::gt() %>%
  gt::fmt_number(columns = "Bound", decimals = 2) %>%
  gt::fmt_number(columns = "Power", decimals = 3)

# Overall (weighted average rates across strata)
orr_ctrl_overall <- orr_ctrl_sub * prevalence + orr_ctrl_complement * (1 - prevalence)
orr_exp_overall <- orr_exp_sub * prevalence + orr_exp_complement * (1 - prevalence)

orr_all <- gsDesign2::fixed_design_rd(
  alpha = alphaHypotheses[6],
  power = NULL,
  p_c = orr_ctrl_overall,
  p_e = orr_exp_overall,
  rd0 = 0,
  n = ceiling(n_total) * 2
)
summary(orr_all) %>%
  gt::gt() %>%
  gt::fmt_number(columns = "Bound", decimals = 2) %>%
  gt::fmt_number(columns = "Power", decimals = 3)

Design list (ordered to match graph hypotheses)

# NULL for non-GSD hypotheses (e.g., ORR tested at a single analysis)
gsD2list <- list(ossub, os, pfssub, pfs, NULL, NULL)

Results entry template

For calendar-time-based interim timing (common interim dates across endpoints), see calendar-time analysis timing.

# Event counts per hypothesis per analysis
events_pfs_all <- c(675, 750)
events_pfs_sub <- c(265, 310)
events_os_all <- c(529, 700, 800)
events_os_sub <- c(185, 245, 295)

inputResults <- tibble(
  H = c(rep(1, 3), rep(2, 3), rep(3, 2), rep(4, 2), 5, 6),
  Pop = c(rep("Subgroup", 3), rep("All", 3),
          rep("Subgroup", 2), rep("All", 2),
          "Subgroup", "All"),
  Endpoint = c(rep("OS", 6), rep("PFS", 4), rep("ORR", 2)),
  nominalP = c(
    .03, .0001, .000001,   # H1: OS Subgroup (3 analyses)
    .2, .15, .1,            # H2: OS All (3 analyses)
    .2, .001,               # H3: PFS Subgroup (2 analyses)
    .3, .2,                 # H4: PFS All (2 analyses)
    .00001,                 # H5: ORR Subgroup (1 analysis)
    .1                      # H6: ORR All (1 analysis)
  ),
  Analysis = c(1:3, 1:3, 1:2, 1:2, 1, 1),
  events = c(events_os_sub, events_os_all,
             events_pfs_sub, events_pfs_all, NA, NA),
  # Spending time: subgroup info fraction for all hypotheses
  spendingTime = c(
    events_os_sub / max(events_os_sub),
    events_os_sub / max(events_os_sub),
    events_pfs_sub / max(events_pfs_sub),
    events_pfs_sub / max(events_pfs_sub),
    NA, NA
  )
)

Sequential p-value computation

EOCtab <- inputResults %>%
  group_by(H) %>%
  slice(1) %>%
  ungroup() %>%
  select("H", "Pop", "Endpoint", "nominalP")
EOCtab$seqp <- .9999

for (EOCtabline in 1:nHypotheses) {
  EOCtab$seqp[EOCtabline] <-
    ifelse(is.null(gsD2list[[EOCtabline]]),
      EOCtab$nominalP[EOCtabline],
      {
        tem <- filter(inputResults, H == EOCtabline)
        gsDesign2::sequential_pval(
          gs_design = gsD2list[[EOCtabline]],
          event = tem$events,
          z = -stats::qnorm(tem$nominalP),
          ustime = tem$spendingTime,
          interval = c(1e-05, 0.9999)
        )
      }
    )
}
EOCtab <- EOCtab %>% select(-"nominalP")

Hypothesis testing with graphicalMCP

result <- graphicalMCP::graph_test_shortcut(
  graph = g0,
  p = EOCtab$seqp,
  alpha = fwer,
  verbose = TRUE
)

adj_p <- result$outputs$adjusted_p
rej <- result$outputs$rejected

# Convert to logical if needed
if (is.numeric(rej)) {
  rej_logical <- rep(FALSE, nHypotheses)
  rej_logical[rej] <- TRUE
} else {
  rej_logical <- as.logical(rej)
}

EOCtab$Rejected <- rej_logical
EOCtab$adjPValues <- adj_p

Verification with updated bounds

Extract graph update sequence

rejected_order <- which(EOCtab$Rejected)[order(EOCtab$adjPValues[EOCtab$Rejected])]

graphs <- list(g0)
if (length(rejected_order) > 0) {
  gu <- graphicalMCP::graph_update(g0, delete = rejected_order)
  graphs <- gu$intermediate_graphs
}

# Get max alpha allocated to each hypothesis
lastWeights <- as.numeric(graphs[[length(graphs)]]$hypotheses)
for (j in seq_along(rejected_order)) {
  h <- rejected_order[j]
  lastWeights[h] <- as.numeric(graphs[[j]]$hypotheses[h])
}
EOCtab$lastAlpha <- fwer * lastWeights

Plot graph sequence with alpha labels

When plotting the graph update sequence, show actual alpha levels (not weights) on each vertex:

for (i in seq_along(graphs)) {
  gi_alpha <- fwer * as.numeric(graphs[[i]]$hypotheses)
  gi_labels <- paste(
    names(graphs[[i]]$hypotheses),
    sprintf("\u03b1 = %s", format(gi_alpha, scientific = FALSE)),
    sep = "\n"
  )
  plot(graphs[[i]],
    layout = layout6,
    vertex.size = 45,
    vertex.label = gi_labels,
    vertex.label.cex = 0.9,
    vertex.color = vertex_colors,
    margin = 0.25
  )
}

Update bounds for verification

for (i in 1:nHypotheses) {
  hresults <- inputResults %>% filter(H == i)
  d2 <- gsD2list[[i]]

  if (!is.null(d2) && EOCtab$lastAlpha[i] > 0) {
    d2_upd <- gsDesign2::gs_update_ahr(
      x = d2,
      alpha = EOCtab$lastAlpha[i],
      ustime = hresults$spendingTime,
      event_tbl = data.frame(
        analysis = hresults$Analysis,
        event = hresults$events
      )
    )
    # Compare nominal p-values to updated bounds
    # Rejected: at least one nominal p <= bound nominal p
    # Not rejected: all nominal p > bound nominal p
  }
}

Spending time rules

Spending time controls how alpha is allocated across analyses. It differs from the information fraction, which drives the correlation structure.

The min() rule

At interim analyses, use the minimum of the planned information fraction and the actual information fraction to protect against over-spending when events accumulate faster than planned:

# Planned info fractions from H1 design
os_bm_info_frac <- ossub$analysis$info_frac
pfs_bm_info_frac <- pfssub$analysis$info_frac

# Actual events observed at each analysis
events_os_bm <- c(185, 245, 295)    # From simulated/observed data
events_pfs_bm <- c(265, 310)

# Spending time: min(planned, actual) at interims; 1 at final
spendingTime_H1 <- pmin(os_bm_info_frac, events_os_bm / max(events_os_bm))
spendingTime_H3 <- pmin(pfs_bm_info_frac, events_pfs_bm / max(events_pfs_bm))

Alignment across populations

The overall population (H2, H4) uses the same spending time as its corresponding subgroup hypothesis (H1, H3). This ensures a consistent, pre-specified spending schedule:

inputResults <- tibble(
  H = c(rep(1, 3), rep(2, 3), rep(3, 2), rep(4, 2), 5, 6),
  spendingTime = c(
    pmin(os_bm_info_frac, events_os_bm / max(events_os_bm)),    # H1: OS BM+
    pmin(os_bm_info_frac, events_os_bm / max(events_os_bm)),    # H2: OS All (same as H1)
    pmin(pfs_bm_info_frac, events_pfs_bm / max(events_pfs_bm)), # H3: PFS BM+
    pmin(pfs_bm_info_frac, events_pfs_bm / max(events_pfs_bm)), # H4: PFS All (same as H3)
    NA, NA                                                        # H5, H6: ORR (not GSD)
  ),
  # ... other columns
)

Stratified ORR power with gs_power_rd()

For the overall population ORR, use gs_power_rd() with stratified inputs and INVAR (inverse-variance) weighting instead of fixed_design_rd():

orr_all <- gsDesign2::gs_power_rd(
  p_c = tibble(stratum = c("BM+", "BM-"),
               rate = c(orr_ctrl_bm_pos, orr_ctrl_bm_neg)),
  p_e = tibble(stratum = c("BM+", "BM-"),
               rate = c(orr_exp_bm_pos, orr_exp_bm_neg)),
  n = tibble(stratum = c("BM+", "BM-"),
             n = c(n_bm, n_bm_neg),
             analysis = c(1, 1)),
  rd0 = 0,
  ratio = 1,
  weight = "invar",
  upper = gsDesign2::gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = alphaHypotheses[6]),
  lower = gsDesign2::gs_b,
  lpar = -Inf,
  test_lower = FALSE,
  info_scale = "h0_h1_info",
  binding = FALSE
)

summary(orr_all) %>%
  gt::gt() %>%
  gt::tab_header(title = "Power for ORR in Overall population (stratified)")

Displaying design bounds

Use gs_bound_summary() to display group sequential bounds in a table:

gsDesign2::gs_bound_summary(ossub) %>%
  gt::gt() %>%
  gt::tab_header(title = "Design for OS in the BM+ population")

Repeated p-values for verification

In the verification step, compute both sequential p-values (cumulative evidence through analysis k) and repeated p-values (evidence at analysis k alone):

for (aa in seq_len(n_analyses)) {
  # Sequential p-value: observed Z at analyses 1..aa, padded with 0.9999 after
  z_seq <- c(all_z[1:aa], rep(-qnorm(0.9999), n_analyses - aa))
  seq_p_by_analysis[aa] <- gsDesign2::sequential_pval(
    gs_design = d2,
    event = n.I,
    z = z_seq,
    ustime = usTime,
    interval = c(1e-05, 0.9999)
  )

  # Repeated p-value: observed Z at analysis aa ONLY, 0.9999 elsewhere
  z_rep <- rep(-qnorm(0.9999), n_analyses)
  z_rep[aa] <- all_z[aa]
  rep_p_by_analysis[aa] <- gsDesign2::sequential_pval(
    gs_design = d2,
    event = n.I,
    z = z_rep,
    ustime = usTime,
    interval = c(1e-05, 0.9999)
  )
}

The sequential p-value is non-decreasing: it uses all evidence through analysis k. The repeated p-value isolates the contribution of a single analysis, analogous to the repeated confidence interval.

Simulation with illness-death model

The vignette template uses an illness-death model for realistic simulation of correlated OS, PFS, and ORR endpoints. The model has four states: 0 (alive, no response/progression), 1 (responded), 2 (progressed), 3 (dead).

See the illness-death model skill for full details. Brief usage:

source("inst/simulation/sim_illness_death.R")
source("inst/simulation/cut_illness_death.R")

# Build transition rates from clinical assumptions
transition_rate <- build_transition_rates(
  strata = c("BM+", "BM-"),
  treatments = c("control", "experimental"),
  median_pfs = c("BM+" = 5, "BM-" = 5),
  median_os = c("BM+" = 12, "BM-" = 12),
  orr = list(
    "BM+" = c(control = 0.15, experimental = 0.30),
    "BM-" = c(control = 0.15, experimental = 0.12)
  ),
  hr_pfs = c("BM+" = 0.65, "BM-" = 1.2),
  hr_os = c("BM+" = 0.70, "BM-" = 1.1)
)

# Simulate one trial
sim_data <- sim_illness_death(
  n = 500,
  stratum = data.frame(stratum = c("BM+", "BM-"), p = c(0.5, 0.5)),
  block = c("control", "control", "experimental", "experimental"),
  enroll_rate = data.frame(rate = c(6.25, 12.5, 18.75, 25), duration = c(2, 2, 2, 12)),
  transition_rate = transition_rate
)

# Determine analysis cut dates
analyses <- list(
  list(min_followup = 6, endpoint = NULL, event_target = NULL,
       target_stratum = NULL, max_followup = NULL),
  list(min_followup = 14, endpoint = "PFS", event_target = 310,
       target_stratum = "BM+", max_followup = 17),
  list(min_followup = 24, endpoint = "OS", event_target = 295,
       target_stratum = "BM+", max_followup = 30)
)
cut_dates <- get_analysis_dates(sim_data, analyses)

# Cut data to ADTTE format
adtte <- lapply(seq_along(cut_dates), function(i) {
  d <- cut_illness_death(sim_data, cut_dates[i])
  d$ANALYSIS <- i
  d
})

P-value computation, theoretical curves, and KM plots

These topics are covered in detail in the illness-death model skill: - logrank_pval() — 1-sided logrank test with stratification pitfall - rd_pval() — Stratified risk difference test (ORR) - Theoretical survival curves from .pfs_cdf() / .os_cdf() - Prevalence-weighted overall population curves - Computing p-values across analyses for the 6-hypothesis template

WPGSD analysis

The wpgsd package accounts for group sequential and population-induced correlations.

Correlation matrix construction

For nested populations (BM+ ⊂ Overall), the event-count matrix D has entries: - D[Hi_As, Hi_At] = events_i at min(s,t) (within-hypothesis, across-analysis) - D[Hi_As, Hj_At] = intersection events at min(s,t) (cross-hypothesis) - Intersection for nested populations = subgroup events at min(s,t)

Bug in generate_corr() for k > 2: The function incorrectly computes within-hypothesis cross-analysis entries for non-adjacent analyses. Use generate_corr() only for k = 2. For k > 2, build D manually and compute corr = diag(1/sqrt(diag(D))) %*% D %*% diag(1/sqrt(diag(D))).

generate_bounds and closed_test

# Sub-graph: 2 hypotheses with full reallocation
m_sub <- matrix(c(0, 1, 1, 0), nrow = 2, byrow = TRUE)
w <- c(0.5, 0.5)  # weights must be > 0 (gsDesign requires alpha > 0)

bound_wpgsd <- wpgsd::generate_bounds(
  type = 3, k = n_analyses, w = w, m = m_sub, corr = corr,
  alpha = total_alpha_for_endpoint,
  sf = list(gsDesign::sfLDOF, gsDesign::sfLDOF),
  sfparm = list(0, 0),
  t = list(info_frac_h1, info_frac_h2)
)

ct <- wpgsd::closed_test(bound_wpgsd, p_obs)

Theoretical basis: Maurer-Bretz framework

Reference: Maurer W, Bretz F. Multiple testing in group sequential trials using graphical approaches. Stat Biopharm Res 2013; 5:311–320.

Algorithm 1: Sequentially rejective graphical procedure for GSD

The Maurer-Bretz algorithm extends the graphical approach of Bretz et al. (2009) to group sequential trials. At each analysis, it operates exactly like the non-sequential graphical procedure but uses group sequential nominal significance levels instead of fixed levels.

  1. Initialize: Set \(t = 1\), \(I = \{1, \ldots, h\}\) (active hypotheses)
  2. At analysis \(t\): Compute nominal p-values \(p_{i,t}\) and group sequential nominal significance levels \(\alpha^*_{i,t}(w_i(I) \cdot \alpha)\) for each \(i \in I\)
  3. If \(p_{j,t} \leq \alpha^*_{j,t}(w_j(I) \cdot \alpha)\) for some \(j\): Reject \(H_j\), update the graph (weights and transitions), repeat step 2
  4. If no rejection and \(t < k\): Continue to next analysis \(t \to t + 1\)

Key insight: The nominal level \(\alpha^*_{i,t}(\gamma)\) is the group sequential boundary at analysis \(t\) for a design with total alpha \(= \gamma\). This is NOT the same as \(w_i(I) \cdot \alpha^*_{i,t}(\alpha)\) — i.e., one cannot simply multiply a fixed-alpha boundary by the weight.

Equivalence with sequential p-values

Algorithm 1 can equivalently be performed using sequential p-values (Liu & Anderson, 2008):

  1. At each analysis \(t\), compute sequential p-values \(p^s_{i,t}\) for each \(i\)
  2. Reject \(H_j\) if \(p^s_{j,t} \leq \alpha \cdot w_j(I)\)
  3. Update graph and repeat

This is exactly what graphicalMCP::graph_test_shortcut(g, p, alpha) does when called with sequential p-values at each analysis.

Well-ordered spending functions

For the procedure to be consonant (sequentially rejective), the spending function family must produce nominal significance levels \(\alpha^*_t(\gamma)\) that are non-decreasing in \(\gamma\). This is called the well-ordering property (Liu & Anderson, 2008).

Sufficient condition: The spent levels \(\alpha_t(\gamma) = a(\gamma, y_t) - a(\gamma, y_{t-1})\) are nondecreasing in \(\gamma\) for all analyses \(t\). For differentiable spending functions, \(\partial^2 a(\gamma, y) / \partial\gamma \partial y \geq 0\) suffices.

Qualified families:

Spending function Well-ordered? Notes
Power: \(a(\gamma,y) = \gamma y^\rho\) Yes (all \(\rho > 0\)) Includes \(\rho=3\) (OBF-like) and \(\rho=1\) (Pocock-like)
Pocock-type: \(a(\gamma,y) = \gamma \ln(1+(e-1)y)\) Yes
OBF-type: \(a(\gamma,y) = 2(1 - \Phi(\Phi^{-1}(1-\gamma/2)/\sqrt{y}))\) Yes for \(\gamma < 0.318\) Covers \(\alpha = 0.025\) and all practical levels
sfLDOF in gsDesign Yes for \(\gamma < 0.318\) Same as OBF-type
sfHSD with \(\gamma < 0\) Verify case by case Generally well-ordered for practical levels

Why naive hierarchical testing fails in GSD

Hung, Wang & O’Neill (2007) showed that the following naive strategy inflates the FWER: “test the secondary endpoint at level \(\alpha\) whenever the primary endpoint is significant at an interim or final analysis.”

The problem: in a GSD with non-binding futility, the primary can be significant at an interim but not at the final analysis (if the trial continues and evidence weakens). Testing the secondary at a later analysis with full \(\alpha\) double-counts the alpha already spent by the primary.

Solution: Use sequential p-values for each hypothesis and apply the graphical procedure (Algorithm 1) at each analysis. This is exactly the approach implemented in the graphicalMCP-gsDesign2 workflow.

Calendar-Time Analysis Timing for OS/PFS/ORR

This note explains how to define analysis timing on the calendar and convert it into spending time for sequential p-value calculations.

Planning vs execution

  • Planning: set a calendar schedule for interim analyses. Use common calendar dates across endpoints where required.
  • Execution: use actual data-cut dates (or a defined cutting routine) to compute realized elapsed times. Keep the same common cut dates for endpoints tied to the same interim.

Design-time timing rules

  • Interim 1 uses a common calendar date for ORR, PFS, and OS.
  • Interim 2 uses a common calendar date for PFS and OS.
  • Final analyses can be endpoint-specific (e.g., OS final later than PFS final).

Example timing table (calendar months from FPFV)

Assume FPFV = month 0.

Endpoint Analysis PlannedMonth Notes
ORR Interim 1 12 Common interim 1 date
PFS Interim 1 12 Common interim 1 date
OS Interim 1 12 Common interim 1 date
PFS Interim 2 18 Common interim 2 date
OS Interim 2 18 Common interim 2 date
PFS Final 24 Endpoint-specific final
OS Final 30 Endpoint-specific final

Customize planned months for your study

Replace the example months with your planned calendar schedule. Use common interim 1 and interim 2 dates across endpoints, then endpoint-specific finals.

t_interim1 <- 12
t_interim2 <- 18
t_final_pfs <- 24
t_final_os <- 30

analysis_timing <- tibble(
  Endpoint = c("ORR", "PFS", "OS", "PFS", "OS", "PFS", "OS"),
  Analysis = c(1, 1, 1, 2, 2, 3, 3),
  PlannedMonth = c(t_interim1, t_interim1, t_interim1,
                   t_interim2, t_interim2, t_final_pfs, t_final_os)
)

Convert calendar time to spendingTime

Define elapsed time as planned months (or actual months at execution). Compute spending time within each endpoint as elapsed time divided by the endpoint-specific final elapsed time. Ensure the final analysis equals 1.

analysis_timing <- tibble(
  Endpoint = c("ORR", "PFS", "OS", "PFS", "OS", "PFS", "OS"),
  Analysis = c(1, 1, 1, 2, 2, 3, 3),
  PlannedMonth = c(12, 12, 12, 18, 18, 24, 30)
)

analysis_timing <- analysis_timing %>%
  group_by(Endpoint) %>%
  mutate(spendingTime = PlannedMonth / max(PlannedMonth)) %>%
  ungroup()

# Join to your results entry for sequential p-values
inputResults <- inputResults %>%
  left_join(analysis_timing, by = c("Endpoint", "Analysis")) %>%
  mutate(spendingTime = ifelse(Endpoint == "ORR", NA, spendingTime))

Execution-time variant using data-cut dates

Use actual cut dates (or a defined cutting routine) to compute elapsed time and spendingTime. Keep the common interim 1 date across ORR, PFS, OS and the common interim 2 date across PFS, OS.

fpfv_date <- as.Date("2025-01-15")

analysis_timing_exec <- tibble(
  Endpoint = c("ORR", "PFS", "OS", "PFS", "OS", "PFS", "OS"),
  Analysis = c(1, 1, 1, 2, 2, 3, 3),
  CutDate = as.Date(c("2026-01-15", "2026-01-15", "2026-01-15",
                      "2026-07-15", "2026-07-15", "2027-01-15", "2027-07-15"))
)

analysis_timing_exec <- analysis_timing_exec %>%
  mutate(ElapsedMonths = as.numeric(difftime(CutDate, fpfv_date, units = "days")) / 30.4375) %>%
  group_by(Endpoint) %>%
  mutate(spendingTime = ElapsedMonths / max(ElapsedMonths)) %>%
  ungroup()

Notes

  • Use endpoint-specific denominators so PFS and OS can share interim dates but still reach 1 at their own finals.
  • If you follow a cutting routine (e.g., simtrial-style calendar cuts), replace PlannedMonth with observed elapsed time at each data cut and recompute spendingTime.
  • The common-date rule at interims ensures synchronized decision points across endpoints at design time.