From f60c5bbbf6478c11458a0de5ad72f8f8d86f0410 Mon Sep 17 00:00:00 2001 From: Lamar Hunt III Date: Mon, 13 Jul 2026 14:31:00 -0500 Subject: [PATCH] Fix group-sequential threshold misalignment for hypotheses first tested after analysis 1 In update_p_thresholds_cpp, per-look rejection boundaries were written into the threshold matrix starting at analysis column 1, while get_maurer_bretz_z_raw_cpp places observed p-values at the hypothesis's analyses_analysed columns. For any multi-look hypothesis whose first tested analysis is not analysis 1 (e.g. OS tested at analyses 2 & 3), the stringent interim boundary landed on an analysis where the hypothesis is never tested and was silently discarded, while the real interim look was tested against the final-look boundary. This inflated simulated unconditional power (Tables 6a/6b) above the hypothesis's full-alpha local power and over-spent alpha; design boundaries and local power (Tables 1/4/5) were unaffected. Thresholds are now indexed by analyses_analysed, mirroring the p-value placement, and forward-filled across later analyses to match the existing forward-fill of observed p-values (so a hypothesis retested after weight propagation uses its most recent look's boundary). Columns before a hypothesis's first look keep the -1 sentinel, matching NA p-values there. Adds regression tests: exact boundary placement for a staggered (analyses_analysed = 2:3) hypothesis, and the invariant that simulated unconditional power cannot exceed full-alpha local power. Co-Authored-By: Claude Fable 5 --- NEWS.md | 6 ++ src/sim_functions.cpp | 27 +++++++-- tests/testthat/test-sim-thresholds.R | 85 ++++++++++++++++++++++++++++ 3 files changed, 112 insertions(+), 6 deletions(-) create mode 100644 tests/testthat/test-sim-thresholds.R diff --git a/NEWS.md b/NEWS.md index 0c30a2e..c558905 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +# appendMCP (development version) + +## Bug Fixes + +- Fixed misalignment of group-sequential rejection boundaries in `update_p_thresholds_cpp`: thresholds are now placed at each hypothesis's `analyses_analysed` columns (matching the observed p-value placement) instead of starting at analysis 1. Previously, any multi-look hypothesis whose first tested analysis was not analysis 1 (e.g. OS tested at analyses 2 & 3) had its interim boundary silently discarded and its interim tested against the final-look boundary, inflating simulated unconditional power (Tables 6a/6b) above the full-alpha local power and over-spending alpha. Design boundaries and local power (Tables 1/4/5) were unaffected. + # appendMCP 0.3.0 ## New Features diff --git a/src/sim_functions.cpp b/src/sim_functions.cpp index 07c1218..6311885 100644 --- a/src/sim_functions.cpp +++ b/src/sim_functions.cpp @@ -922,6 +922,7 @@ NumericMatrix update_p_thresholds_cpp(DataFrame analyses, DataFrame hypotheses, const IntegerVector hyp_index = hypotheses["index"]; const NumericVector possible_weight = hypotheses["possible_weight"]; const List p_thresholds = hypotheses["p_thresholds"]; + const List analyses_analysed = hypotheses["analyses_analysed"]; for (int j = 0; j < n_hyp; j++) { if (weights[j] > 0) { @@ -940,14 +941,28 @@ NumericMatrix update_p_thresholds_cpp(DataFrame analyses, DataFrame hypotheses, } } if (best_h >= 0) { + // Place each per-look threshold at the analysis it belongs to + // (analyses_analysed), mirroring how get_maurer_bretz_z_raw_cpp + // places observed p-values. const NumericVector p_thr_j = p_thresholds[best_h]; - const int len = std::min(static_cast(p_thr_j.size()), n_analyses); - for (int a = 0; a < len; a++) { - p_thr(j, a) = p_thr_j[a]; + const IntegerVector aa = analyses_analysed[best_h]; + const int len = std::min(static_cast(p_thr_j.size()), + static_cast(aa.size())); + for (int i = 0; i < len; i++) { + const int col = aa[i] - 1; // convert to 0-based analysis column + if (col >= 0 && col < n_analyses) { + p_thr(j, col) = p_thr_j[i]; + } } - if (len > 0) { - const double last_val = p_thr_j[len - 1]; - for (int a = len; a < n_analyses; a++) { + // Forward-fill between/after looks: p_obs is forward-filled the same + // way, so a hypothesis retested at a later analysis (after weight + // propagation) uses its most recent look's boundary. Columns before + // the first look stay -1; p_obs is NA there. + double last_val = -1.0; + for (int a = 0; a < n_analyses; a++) { + if (p_thr(j, a) >= 0.0) { + last_val = p_thr(j, a); + } else if (last_val >= 0.0) { p_thr(j, a) = last_val; } } diff --git a/tests/testthat/test-sim-thresholds.R b/tests/testthat/test-sim-thresholds.R new file mode 100644 index 0000000..800c604 --- /dev/null +++ b/tests/testthat/test-sim-thresholds.R @@ -0,0 +1,85 @@ +library(appendMCP) + +# Regression tests for threshold/analysis alignment in the simulated operating +# characteristics: per-look group-sequential boundaries must be placed at each +# hypothesis's analyses_analysed columns, matching where +# get_maurer_bretz_z_raw_cpp places the observed p-values. A hypothesis whose +# first look is not analysis 1 (e.g. OS tested at analyses 2 & 3) previously +# had its interim boundary written to analysis 1 and its final boundary applied +# at the interim, inflating simulated power above the full-alpha local power. + +staggered_config <- list( + study_name = "staggered_looks", alpha = 0.025, sims = 4000L, + analyses = tibble::tribble( + ~endpoint, ~strata, ~treatments, ~sample_size, ~events, + "PFS", c("All"), c("Control", "Treatment"), NA, 300, + "PFS", c("All"), c("Control", "Treatment"), NA, 420, + "OS", c("All"), c("Control", "Treatment"), NA, 300), + hypotheses = tibble::tribble( + ~type, ~endpoint, ~strata, ~control, ~test, ~analyses_analysed, ~sf, ~sfpar, ~nominal, ~test_method, + "Primary", "PFS", c("All"), "Control", "Treatment", 1:2, "asOF", NULL, NULL, "logrank", + "Primary", "OS", c("All"), "Control", "Treatment", 2:3, "asOF", NULL, NULL, "logrank"), + enroll_rate = tibble::tribble( + ~stratum, ~treatments, ~rate, ~duration, ~ratio, + "All", c("Control", "Treatment"), 60, 18, c(1, 1)), + distribution_tte = tibble::tribble( + ~endpoint, ~stratum, ~treatment, ~duration, ~fail_rate, ~dropout_rate, + "PFS", "All", "Control", Inf, log(2)/10, 0, + "PFS", "All", "Treatment", Inf, log(2)/15, 0, + "OS", "All", "Control", Inf, log(2)/23, 0, + "OS", "All", "Treatment", Inf, log(2)/31, 0), + graph = list(g = rbind(c(0, 1), c(0, 0)), w = c(1, 0))) + +set.seed(1) +staggered_result <- process_config(staggered_config) + +test_that("thresholds are placed at each hypothesis's analyses_analysed columns", { + hyp <- process_rejection_rules( + dplyr::mutate(staggered_result$hypotheses, information = information_factor), + alpha = staggered_config$alpha + ) + + # H2 (OS, looks at analyses 2 & 3) holding full alpha + M2 <- update_p_thresholds_cpp(staggered_result$analyses, hyp, + staggered_config$graph$g, c(0, 1)) + h2_bounds <- hyp$p_thresholds[hyp$index == 2 & hyp$possible_weight == 1][[1]] + expect_length(h2_bounds, 2) + # Interim boundary at analysis 2, final boundary at analysis 3 + expect_equal(M2[2, 2], h2_bounds[1]) + expect_equal(M2[2, 3], h2_bounds[2]) + # H2 has no look at analysis 1: threshold must stay at the -1 sentinel + expect_equal(M2[2, 1], -1) + # Sanity: alpha-spending interim boundary is stricter than the final one + expect_lt(M2[2, 2], M2[2, 3]) + + # H1 (PFS, looks at analyses 1 & 2) holding full alpha: placement unchanged, + # with the final boundary carried forward to analysis 3 to mirror the + # forward-fill of observed p-values + M1 <- update_p_thresholds_cpp(staggered_result$analyses, hyp, + staggered_config$graph$g, c(1, 0)) + h1_bounds <- hyp$p_thresholds[hyp$index == 1 & hyp$possible_weight == 1][[1]] + expect_equal(M1[1, 1], h1_bounds[1]) + expect_equal(M1[1, 2], h1_bounds[2]) + expect_equal(M1[1, 3], h1_bounds[2]) +}) + +test_that("simulated power of a late-starting hypothesis respects its local power", { + # Unconditional (simulated) power can never exceed the hypothesis's local + # power at full alpha; before the alignment fix, H2's interim was tested + # against its final-look boundary and exceeded it by ~10 percentage points. + table5 <- as.data.frame(staggered_result$tables$table5) + table6a <- as.data.frame(staggered_result$tables$table6a) + + local_h2 <- table5[grepl("H2", table5$Hypothesis), ] + local_h2 <- local_h2[order(local_h2$Analysis), "Local power"] + + uncond_h2 <- table6a[table6a$Metric == "Power" & + table6a$`Hypothesis subset` == "H2", ] + uncond_h2 <- uncond_h2[order(uncond_h2$Analysis), ] + + expect_equal(uncond_h2$Value[uncond_h2$Analysis == 1], 0) + + mc_tol <- 0.025 # ~3 binomial SDs at 4000 sims; the pre-fix violation was ~0.10 + expect_lte(uncond_h2$Value[uncond_h2$Analysis == 2], local_h2[1] + mc_tol) + expect_lte(uncond_h2$Value[uncond_h2$Analysis == 3], local_h2[2] + mc_tol) +})