--- title: "MRCT with complex survival analysis design" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{MRCT with complex survival analysis design} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{r load, message = FALSE} library(FastSurvival) ``` ## What this vignette shows A multiregional clinical trial (MRCT) is run under a single protocol across more than one region to support simultaneous approvals. Beyond the overall treatment effect, a regional regulator often asks whether the effect in its own region is consistent with the overall result. The Japanese Ministry of Health, Labour and Welfare set out two criteria for this purpose, usually called Method 1 and Method 2, and the operating characteristic of interest is the *regional consistency probability*: the chance that the regional criterion is met, often evaluated conditional on the overall test being significant. Closed-form formulas for these probabilities exist for continuous endpoints and have been extended to other endpoint types under a normal approximation of the treatment effect. Those formulas are convenient but they assume a fixed design, a single final analysis, and balanced, simultaneous accrual across regions. Real oncology MRCTs routinely break those assumptions: the design is group-sequential with interim looks, and one region (here, Japan) starts enrolling later than the others, so its information accrues on a different calendar. When the assumptions behind the closed-form formulas no longer hold, simulation is the natural way to obtain the consistency probabilities. FastSurvival does not provide a regional-consistency function, and it is not meant to: three-or-more-region designs and consistency criteria are outside its scope. What it does provide are fast, validated building blocks. This vignette shows that the existing simulation trio, `simdata_fast()`, `analysis_fast()`, and `simsummary_fast()`, can be composed to evaluate the conditional regional consistency probabilities of Method 1 and Method 2 at every interim and final analysis, under a group-sequential design with a late-enrolling region. The only piece written by hand is a short post-processing layer in base R that applies the consistency criteria to the per-region estimates the trio returns. ## The two criteria Write the estimated treatment effect on the hazard-reduction scale, that is `1 - HR` for the entire trial population and `1 - HR_s` for region `s`, where a positive value indicates benefit (this is the survival-endpoint formulation of Teng et al., 2018). The region of interest is region 1, and `pi` is the effect-retention fraction, conventionally `0.5`. Method 1 asks that region 1 retain at least a fraction `pi` of the overall effect, `(1 - HR_1) > pi * (1 - HR)`. Method 2 asks that the effect trend in the benefit direction in every region, `1 - HR_s > 0` for all `s`, equivalently that every regional hazard ratio fall below 1. The *conditional* regional consistency probability is the probability of meeting the criterion given a statistically significant overall result. In a group-sequential design the natural counterpart is, at each look, the probability of meeting the criterion among the trials that stop for efficacy at that look. ## The scenario A three-region phase III oncology superiority trial compares a treatment group against control on overall survival, with control median 4.3 months and treatment median 5.811 months (hazard ratio about 0.74). Region 1 is Japan, a small region (10% of the patients) whose accrual starts three months after the other regions. Enrollment runs over a 12.5-month window and survival is exponential. ```{r scenario} NSIM <- 10000 # number of simulated trials SEED <- 1 mst.T <- 5.811 # treatment median survival mst.C <- 4.3 # control median survival N.T <- c(25, 112, 113) # treatment sample size per region (region 1 = Japan) N.C <- c(25, 112, 113) # control sample size per region a <- 12.5 # accrual duration a_delay <- 3 # accrual delay for region 1 (Japan starts later than others) PI <- 0.5 # effect-retention fraction for Method 1 E <- c(142, 248, 354) # target total events at the three looks ``` The design is group-sequential with three looks at roughly 40%, 70%, and 100% of the final event count. The first look is futility-only, the second is efficacy-only, and the final look tests both. Boundaries are on the overall log-rank `Z` scale, where treatment benefit is a negative `Z`; in practice they would come from `gsDesign` or `rpact`. A boundary set to `NA` omits that rule at that look. ```{r boundaries} eff.bound <- c(NA, -2.437, -2) # lower efficacy boundary on logrank.z fut.bound <- c(0.381, NA, -2) # upper (non-binding) futility boundary ``` ## Step 1: simulate the trial with a region-specific accrual delay The accrual specification in `simdata_fast()` applies to the whole call, so a region-specific delay is expressed by generating each region with its own accrual window and stacking the results. Region 1 accrues uniformly over `[a_delay, a]`; the other regions over `[0, a]`. Each region is generated from an independent stream (`seed + k`), which is the right model for independent regional subpopulations. A `subgroup` column tags the region, and the shared `sim` index links the regional blocks into one trial per simulation. ```{r simulate} K <- length(N.T) blocks <- lapply(seq_len(K), function(k) { a.time_k <- if (k == 1L) c(a_delay, a) else c(0, a) dat <- simdata_fast( nsim = NSIM, n = c(N.C[k], N.T[k]), a.time = a.time_k, a.prop = 1, e.median = list(mst.C, mst.T), seed = SEED + k ) dat$subgroup <- k dat }) data_all <- do.call(rbind, blocks) head(data_all) ``` ## Step 2: one event-driven, by-subgroup analysis pass `analysis_fast()` with `event.looks` cuts each simulated trial at the calendar time of the target event over the whole trial population, exactly the timing a group-sequential MRCT uses. With `by.subgroup = TRUE` the same cutoff is reused for the overall analysis and for each region, so the regional estimates are taken at the identical data cut as the overall test. We request the log-rank statistic for the overall efficacy decision and the Cox hazard ratio for the consistency criteria. ```{r analyze} res <- analysis_fast( data = data_all, control = 1, event.looks = E, stat = c("logrank", "coxph"), side = 1, by.subgroup = TRUE ) head(res[, c("sim", "look", "population", "n.event", "logrank.z", "cox.hr")]) ``` The output is in long form: for each `(sim, look)` there is one `overall` row and one row per region (`subgroup_1`, `subgroup_2`, `subgroup_3`). ## Step 3: assemble the per-look effect estimates and the two indicators A few lines of base R reshape the long output into one row per `(sim, look)` carrying the overall log-rank `Z`, the overall hazard ratio, the region-1 hazard ratio, and the Method 1 and Method 2 indicators. Where a region has no events yet at an early cut (a real possibility for the late-starting region), its hazard ratio is `NA`; such an undemonstrated criterion is treated as not met. ```{r indicators} ov <- res[res$population == "overall", ] key <- function(d) paste(d$sim, d$look) # Region-1 hazard ratio aligned to the overall rows r1 <- res[res$population == "subgroup_1", ] ov$HR.1 <- r1$cox.hr[match(key(ov), key(r1))] # Method 2: all regional hazard ratios below 1 at this (sim, look) reg <- res[res$population != "overall", ] all_below_1 <- tapply(reg$cox.hr, key(reg), function(h) all(h < 1)) ov$M2 <- as.integer(all_below_1[key(ov)]) # Method 1: region 1 retains at least a fraction PI of the overall effect ov$M1 <- as.integer((1 - ov$HR.1) > PI * (1 - ov$cox.hr)) # Undemonstrated criteria (NA hazard ratio) count as not met ov$M1[is.na(ov$M1)] <- 0L ov$M2[is.na(ov$M2)] <- 0L ``` ## Step 4: apply the group-sequential boundaries The looks are arranged into a simulation-by-look matrix and the boundaries are applied in sequence. Each trial stops at its first boundary crossing; when both boundaries are crossed at the same look, efficacy takes precedence. ```{r gsd} sims <- sort(unique(ov$sim)); nsim <- length(sims) looks <- sort(unique(ov$look)); nlk <- length(looks) cell <- cbind(match(ov$sim, sims), match(ov$look, looks)) mk <- function(v) { m <- matrix(NA_real_, nsim, nlk); m[cell] <- v; m } Zm <- mk(ov$logrank.z) M1m <- mk(ov$M1) M2m <- mk(ov$M2) Cutm <- mk(ov$cutoff) effB <- matrix(eff.bound, nsim, nlk, byrow = TRUE) futB <- matrix(fut.bound, nsim, nlk, byrow = TRUE) eff_cross <- !is.na(effB) & Zm <= effB fut_cross <- !is.na(futB) & Zm >= futB stop_any <- eff_cross | fut_cross first_stop <- apply(stop_any, 1L, function(r) { w <- which(r); if (length(w)) w[1L] else nlk }) reject <- eff_cross[cbind(seq_len(nsim), first_stop)] ``` ## Step 5: conditional regional consistency probabilities For each look, restrict to the trials that stop for efficacy there and take the fraction that also meet each criterion. `Efficacy` is the marginal probability of an efficacy stop at the look, `CON_M1` and `CON_M2` are the conditional consistency probabilities given an efficacy stop, and `JOI_M1` and `JOI_M2` are the joint probabilities of an efficacy stop together with consistency. ```{r consistency} con1 <- con2 <- joi1 <- joi2 <- eff <- atime <- rep(NA_real_, nlk) for (k in seq_len(nlk)) { atime[k] <- mean(Cutm[, k], na.rm = TRUE) # calendar time of look k (all trials) sel <- reject & first_stop == k # trials stopping for efficacy at look k eff[k] <- mean(sel) joi1[k] <- sum(M1m[sel, k]) / nsim # joint prob (0 when no efficacy stop) joi2[k] <- sum(M2m[sel, k]) / nsim if (any(sel)) { # conditional prob (undefined otherwise) con1[k] <- mean(M1m[sel, k]) con2[k] <- mean(M2m[sel, k]) } } consistency <- data.frame( look = looks, Info.Fraction = E / max(E), Events = E, Analysis_time = atime, Efficacy = eff, cum.power = cumsum(eff), CON_M1 = con1, JOI_M1 = joi1, CON_M2 = con2, JOI_M2 = joi2 ) consistency ``` The first look is futility-only (`eff.bound[1] = NA`), so no trial can stop for efficacy there. The efficacy-conditional quantities are therefore degenerate at look 1: the joint probabilities `JOI_M1` and `JOI_M2` are zero, while the conditional probabilities `CON_M1` and `CON_M2` are undefined (`NA`, a zero-over-zero ratio). `Analysis_time` is the mean calendar time of the look itself, taken over all trials, and so is reported at every look. At the later looks the table gives, per analysis, the conditional regional consistency probability for each method. Method 2 is the looser criterion and so its conditional probability runs higher than Method 1's at both efficacy looks. For each method the conditional consistency is comparable at the two efficacy looks, if anything slightly higher at the interim, because the trials that cross the stringent interim efficacy boundary are enriched for a strong overall effect, which also makes the regional criteria more likely to be met. This calendar-driven behavior, together with the late-starting region 1 whose events accrue on a different schedule, is exactly what the closed-form formulas, which assume a single final analysis and simultaneous accrual, cannot capture. ## Tying back to the trio: standard operating characteristics The third trio function, `simsummary_fast()`, summarizes the same boundaries on the overall population to give the usual group-sequential operating characteristics, here as an independent check on the efficacy column above. ```{r summary} gs <- simsummary_fast( res[res$population == "overall", ], eff.col = "logrank.z", efficacy = eff.bound, fut.col = "logrank.z", futility = fut.bound, direction = "lower" ) gs ``` The per-look `prob.stop.efficacy` from `simsummary_fast()` matches the `Efficacy` column computed by hand, and its `overall` row gives the cumulative power and the expected information and calendar time at stopping. ## Notes and extensions The hazard ratio used for the consistency criteria is the Cox estimate from `coxph_fast()`. To use the log-scale criterion of Teng et al. (2018) instead, replace `cox.hr` with `cox.coef` and test `cox.coef_1 < PI * cox.coef_overall` for Method 1 and `cox.coef_s < 0` for all regions for Method 2. The conditional probability here is defined on the first efficacy-stopping look, which is the standard group-sequential interpretation; the joint and marginal columns are provided alongside so other conditioning conventions are easy to derive. The same composition extends naturally. More regions are added by lengthening the sample-size vectors. A different consistency endpoint, for example a restricted-mean or milestone-survival contrast, is obtained by requesting that statistic in `analysis_fast()` and reading its per-region estimate in place of the hazard ratio. Non-proportional hazards are handled by switching the overall efficacy statistic to a weighted log-rank or max-combo test. None of this needs a new function; it is the existing building blocks recombined. ## References Ministry of Health, Labour and Welfare. (2007). *Basic principles on global clinical trials* (Notification No. 0928010). Teng, Z., Lin, J., & Zhang, B. (2018). Practical recommendations for regional consistency evaluation in multi-regional clinical trials with different endpoints. *Statistics in Biopharmaceutical Research*, 10(1), 50-56. Homma, G. (2024). Cautionary note on regional consistency evaluation in multiregional clinical trials with binary outcomes. *Pharmaceutical Statistics*, 23(3), 385-398.