Simulated power inflated for multi-stratum designs: enrollment pooling in get_dist_tite/get_info_tite drops per-stratum periods

[lhunt13]

Package: appendMCP 0.3.0 (main @ 435200e). Affected functions: get_dist_tite() (R/utils-sim.R) and get_info_tite() (R/utils-process.R), feeding exec_sims → Tables 6a/6b. Related to but distinct from #2 (threshold alignment); the numbers below are measured on top of the #2 fix.

Summary

When get_dist_tite()/get_info_tite() pool arm-level enrollment down to stratum totals for ahr_dd()/gsDesign2::ahr(), they group by (rate, duration) only:

enroll_rate |>
  dplyr::group_by(.data$rate, .data$duration) |>
  dplyr::mutate(rate = sum(.data$arm_rate)) |>
  dplyr::slice_head(n = 1)

Enrollment periods with identical (rate, duration) values in different strata collapse into a single row that keeps only the first stratum's label. In multi-stratum designs where strata share period values — most commonly the zero-rate lead-in periods of a sub-study that opens late — all but the first stratum lose those periods and are treated as enrolling from time 0.

The affected strata then have far more follow-up by each analysis time, so the simulation's expected events, statistical information, and z-statistic means (EZ = -log(ahr)*sqrt(info)) are all inflated. Simulated power (Tables 6a/6b) for TTE hypotheses can substantially exceed the full-alpha local power (Table 5), which is an upper bound.

get_boundaries() pools the same data correctly (group_by(stratum, rate, duration, index)), which is why the design boundaries and local power are unaffected — and why the discrepancy shows up as sim-vs-Table-5 disagreement.

Demonstration in the shipped platform_study example

H8 (OS, sub-study 2, opens at month 14 via two zero-rate lead-in periods shared by strata S2_Bpos/S2_Bneg, tested at analyses 2 & 3):

• After the buggy pooling, S2_Bpos keeps its 14-month delay but S2_Bneg enrolls from month 0.
• ahr_dd() then expects 142.2 events at H8's first look instead of the correct 104.8, information 26.2 instead of 20.1, EZ 2.44 instead of 2.14.
• Feeding ahr_dd() the correctly pooled per-stratum enrollment reproduces gsDesign2::gs_power_ahr's reference values (events 104.77/148.46, info 20.103/27.403, EZ 2.1434/2.5024) to machine precision — for H8 and for every other TTE hypothesis in the example.
• Simulated power for H8 at its first look: 0.55 before the fix vs 0.43 after (full-alpha local power: 0.35 — see "remaining question" below for the residual).

Every hypothesis in a late-opening sub-study (H3–H11) is affected; H1/H2 (sub-study 1, opens at time 0) and H12 are bit-for-bit unchanged by the fix, matching the mechanism.

Fix

Pool per stratum and enrollment period, exactly mirroring get_boundaries():

enroll_rate |>
  dplyr::group_by(.data$stratum, .data$rate, .data$duration, .data$index) |>
  dplyr::summarise(rate = sum(.data$arm_rate), .groups = "drop") |>
  dplyr::arrange(.data$index)

PR with this fix and regression tests (anchored against gsDesign2::ahr() with explicitly correct stratified enrollment) to follow.

Remaining question: information scale for simulated z-statistics

After this fix, a smaller systematic gap remains between simulated power and Table 5 local power, e.g. in platform_study (first look, full alpha):

Hyp	Table 5 local	simulated (both fixes)
H1 (OS)	0.723	0.744
H8 (OS)	0.350	0.433
H12 (OS, single look)	0.565	0.611

This residual is fully explained by a modeling-convention difference, not by sampling error or either bug: get_dist_tite() simulates Z ~ N(-log(ahr)*sqrt(info), 1) (information under the alternative), whereas gs_power_ahr's crossing probabilities imply a smaller drift, close to -log(ahr)*sqrt(info0) (null/Schoenfeld scale). For H12 (single look, so the comparison is a one-dimensional marginal): implied drift 2.1224 vs theta*sqrt(info0) = 2.1126 vs the sim's theta*sqrt(info) = 2.2483. An MVN computation using the sim's own EZ/CovZ reproduces the simulated powers to three decimals, so the drift convention is the entire remainder.

Since info >= info0 in these designs, the current choice makes the simulation systematically optimistic relative to the package's own design engine. Should get_dist_tite() use info0 (making the simulation consistent with gs_power_ahr), or is the H1-information scale intentional? Happy to submit a follow-up PR either way.