--- title: "Group-sequential design with correlated PFS and OS under the Fleischer model" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Group-sequential design with correlated PFS and OS under the Fleischer model} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6.5, fig.height = 4.2, fig.align = "center" ) library(FastSurvival) has_gsd <- requireNamespace("gsDesign", quietly = TRUE) ``` ## Overview Many confirmatory oncology trials test progression-free survival (PFS) as the primary endpoint and overall survival (OS) as a key secondary endpoint, with the two endpoints tested in a fixed sequence to control the family-wise error rate. Because every PFS event is a death or a progression and every death is also a PFS event, the two endpoints are positively correlated, and that correlation is needed to characterize the operating characteristics of the sequential procedure. This vignette simulates such a design with the simulation trio `simdata_fast`, `analysis_fast`, and `simsummary_fast`. It generates correlated PFS and OS times from the Fleischer maximal-independence model, runs an event-driven group-sequential design with one futility analysis followed by two efficacy analyses, evaluates the power of each endpoint under the alternative, and estimates the correlation of the standardized log-rank statistics `Corr(Z_PFS, Z_OS)`. The empirical correlation reported here is the quantity for which a closed-form expression is derived in a companion methodological note; the simulation provides an independent check of that derivation and of the design's operating characteristics. ### The Fleischer model For a subject in group `g` the Fleischer maximal-independence model (Fleischer, Gaschler-Markefski, and Bluhmki, 2009) takes an independent time-to-progression `TTP ~ Exp(lam1)` and an overall survival time `OS ~ Exp(lam2)`, and sets ``` PFS = min(TTP, OS), OS = OS. ``` Since `TTP` and `OS` are independent exponentials, `PFS` is exponential with hazard `Lambda = lam1 + lam2`. Given the medians of PFS and OS, the OS hazard is `lam2 = log(2) / median_OS` and the PFS hazard is `Lambda = log(2) / median_PFS`, so the implied TTP hazard is `lam1 = Lambda - lam2`, which must be positive (median OS at least median PFS in each group). The dependence between the two endpoints is induced entirely through the shared `OS` component. ### The design The trial is event-driven with three synchronized analyses indexed by cumulative PFS event counts. The first analysis is a non-binding futility look on PFS; the second and third are the efficacy interim and the efficacy final. PFS is the primary endpoint and OS is the key secondary endpoint, and the two are tested in the order PFS then OS: OS is declared significant only if PFS has already been declared significant at the same or an earlier look. The analyses are timed by PFS events, and the OS analysis at each look reuses the same calendar cutoff as the PFS analysis. The PFS efficacy boundary uses an O'Brien-Fleming Lan-DeMets alpha-spending function and the OS efficacy boundary uses a Pocock-type Lan-DeMets spending function, each at a one-sided level of 0.025 and each spent over the two efficacy looks. All tests are the ordinary unweighted log-rank test, one-sided in the direction of treatment benefit. ## Design parameters ```{r design} nsim <- 5000 seed <- 20260611 # Per-group sample sizes (group 1 = control, group 2 = treatment) n_control <- 300 n_treat <- 300 r_alloc <- n_treat / n_control # allocation ratio treatment:control # Piecewise-uniform accrual: 20 per month, then 40 per month (600 enrolled by month 18) a_time <- c(0, 6) a_rate <- c(20, 40) # Median PFS and OS by group, in months mst_PFS_C <- 6 mst_PFS_T <- 9 mst_OS_C <- 14 mst_OS_T <- 18 # Exponential hazards implied by the medians lam_PFS_C <- log(2) / mst_PFS_C lam_PFS_T <- log(2) / mst_PFS_T lam_OS_C <- log(2) / mst_OS_C lam_OS_T <- log(2) / mst_OS_T # Fleischer maximal-independence: PFS hazard = TTP hazard + OS hazard, # so the latent TTP hazard is the PFS hazard minus the OS hazard. lam_TTP_C <- lam_PFS_C - lam_OS_C lam_TTP_T <- lam_PFS_T - lam_OS_T stopifnot(lam_TTP_C > 0, lam_TTP_T > 0) # Non-informative exponential dropout, 10 percent per year on the month scale eta <- -log(1 - 0.10) / 12 # Event-driven looks (cumulative PFS event counts): # look 1 = futility, look 2 = efficacy interim, look 3 = efficacy final d_PFS <- c(200, 300, 400) L <- length(d_PFS) # One-sided overall significance level alpha_total <- 0.025 c(HR_PFS = lam_PFS_T / lam_PFS_C, HR_OS = lam_OS_T / lam_OS_C) ``` The PFS hazard ratio is the stronger effect and the OS hazard ratio is more modest, which is the usual pattern when a treatment delays progression more than it extends survival. ## Data generation A single `simdata_fast` call generates the two correlated endpoints directly from the illness-death structure. The intermediate event (progression) has hazard `h01 = lam_TTP` and the terminal event (death) has hazard `h02 = lam_OS`; the post-progression death hazard is left at its default, equal to `h02`, which reproduces the maximal-independence model described above. The first endpoint `e1` is PFS and the terminal endpoint `e2` is OS. The dropout time is shared between the two endpoints, so a subject who drops out is censored for both PFS and OS at the same time. ```{r generate} # A single call returns both correlated endpoints. e1 is PFS (the state-0 exit: # progression or death), e2 is OS (the terminal event), and dropout_time is the # shared non-informative dropout. The observed columns e1_tte / e1_event and # e2_tte / e2_event are each already censored at the shared dropout time. df <- simdata_fast( nsim = nsim, n = c(n_control, n_treat), a.time = a_time, a.rate = a_rate, h01.hazard = list(lam_TTP_C, lam_TTP_T), h02.hazard = list(lam_OS_C, lam_OS_T), d.hazard = eta, seed = seed ) # Observed PFS endpoint (e1). pfs_dat <- data.frame( sim = df$sim, group = df$group, accrual_time = df$accrual_time, tte = df$e1_tte, event = df$e1_event ) # Observed OS endpoint (e2). os_dat <- data.frame( sim = df$sim, group = df$group, accrual_time = df$accrual_time, os_tte = df$e2_tte, os_event = df$e2_event ) ``` ## Primary endpoint analysis (PFS) PFS is analyzed at the three event-driven looks. With `event.looks` the calendar cutoff at look `l` is the time of the `d_PFS[l]`-th PFS event over the whole trial, and `analysis_fast` returns that cutoff per simulation. These cutoffs are exactly the analysis times of the design, and the OS analysis below reuses them. ```{r pfs-analysis} pfs_res <- analysis_fast( pfs_dat, control = 1, event.looks = d_PFS, stat = "logrank", side = 1 ) # Per-simulation calendar cutoffs A_l (rows = simulation, columns = look) A_mat <- matrix(pfs_res$cutoff, nrow = nsim, ncol = L, byrow = TRUE) # With this sample size and these event targets the looks are reached in # essentially every simulation; warn if any target is not reached. if (!all(pfs_res$reached)) { warning("Some PFS event targets were not reached; increase n or lower d_PFS.") } # Standardized PFS log-rank statistics (natural sign: benefit is negative) Z_PFS <- matrix(pfs_res$logrank.z, nrow = nsim, ncol = L, byrow = TRUE) ``` ## Key secondary endpoint analysis (OS) OS is analyzed at the same per-simulation cutoffs `A_l`. For each look the OS times are administratively censored at `A_l`, and the one-sided log-rank statistic is computed with `analysis_fast` at a single, deliberately large calendar time so that no further censoring is applied to the already-cut data. This reproduces the design rule that the OS analysis at a look uses the PFS-driven analysis time. ```{r os-analysis} big <- max(A_mat, na.rm = TRUE) + 1 os_blocks <- vector("list", L) for (l in seq_len(L)) { A_l <- A_mat[os_dat$sim, l] # per-row cutoff at look l follow <- pmax(0, A_l - os_dat$accrual_time) # administrative follow-up obs_t <- pmin(os_dat$os_tte, follow) # observed OS time at the look obs_e <- os_dat$os_event * as.integer(os_dat$accrual_time + os_dat$os_tte <= A_l) cut_l <- data.frame( sim = os_dat$sim, group = os_dat$group, accrual_time = 0, tte = obs_t, event = obs_e ) res_l <- analysis_fast( cut_l, control = 1, time.looks = big, stat = "logrank", side = 1 ) # Relabel as look l and report the PFS-driven cutoff for the timing summary; # keep only the columns that are meaningful for the censored OS data. res_l$look <- l res_l$look.value <- d_PFS[l] res_l$cutoff <- A_mat[, l] os_blocks[[l]] <- res_l[, c("sim", "look", "look.value", "cutoff", "n.event", "logrank.z", "logrank.chisq", "logrank.p")] } os_res <- do.call(rbind, os_blocks) os_res <- os_res[order(os_res$sim, os_res$look), ] # Standardized OS log-rank statistics Z_OS <- matrix(NA_real_, nsim, L) for (l in seq_len(L)) Z_OS[, l] <- os_blocks[[l]]$logrank.z ``` The mean PFS and OS event counts at each look summarize the information available to each endpoint. PFS reaches its event targets by construction, while OS, being the slower endpoint, accrues far fewer events at the same analysis times. ```{r event-counts} d_PFS_emp <- colMeans(matrix(pfs_res$n.event, nrow = nsim, ncol = L, byrow = TRUE)) d_OS_emp <- colMeans(matrix(os_res$n.event, nrow = nsim, ncol = L, byrow = TRUE)) A_mean <- colMeans(A_mat) data.frame( look = seq_len(L), d_PFS_target = d_PFS, d_PFS_mean = round(d_PFS_emp), d_OS_mean = round(d_OS_emp), cutoff_months = round(A_mean, 1) ) ``` ## Efficacy and futility boundaries Efficacy is assessed at the two later looks only. Two-look Lan-DeMets alpha-spending boundaries are computed with `gsDesign` at the information fractions implied by the event counts: the PFS information fractions are the exact event-target ratios, and the OS information fractions use the mean OS event counts. `gsDesign` reports its boundaries on the convention that a positive Z favors the treatment, so they are negated to match the natural sign of `logrank.z`, where treatment benefit is a negative value. The first look carries a non-binding futility rule on PFS, expressed as a threshold on the standardized statistic. ```{r boundaries, eval = has_gsd} timing_PFS_eff <- c(d_PFS[2] / d_PFS[3], 1) timing_OS_eff <- c(d_OS_emp[2] / d_OS_emp[3], 1) gs_PFS <- gsDesign::gsDesign( k = 2, test.type = 1, alpha = alpha_total, timing = timing_PFS_eff, sfu = gsDesign::sfLDOF ) gs_OS <- gsDesign::gsDesign( k = 2, test.type = 1, alpha = alpha_total, timing = timing_OS_eff, sfu = gsDesign::sfLDPocock ) b_PFS <- gs_PFS$upper$bound # positive z-boundaries (interim, final) b_OS <- gs_OS$upper$bound # Efficacy boundaries on the natural-sign scale; NA at the futility-only look 1. eff_PFS <- c(NA, -b_PFS[1], -b_PFS[2]) eff_OS <- c(NA, -b_OS[1], -b_OS[2]) # Non-binding PFS futility at look 1: stop if the observed PFS hazard ratio is at # or above futility_HR. On the natural-sign scale this is logrank.z at or above # log(futility_HR) * sqrt(r * d) / (1 + r). futility_HR <- 1.0 fut1 <- log(futility_HR) * sqrt(r_alloc * d_PFS[1]) / (1 + r_alloc) fut_PFS <- c(fut1, NA, NA) data.frame( look = seq_len(L), PFS_efficacy = round(eff_PFS, 3), PFS_futility = round(fut_PFS, 3), OS_efficacy = round(eff_OS, 3) ) ``` ## Operating characteristics The marginal operating characteristics of each endpoint are obtained with `simsummary_fast`. For PFS the efficacy boundary and the look-1 futility boundary are applied together on `logrank.z` with `direction = "lower"`, so a trial stops for efficacy when the statistic is at or below the efficacy boundary and for futility when it is at or above the futility boundary. For OS the efficacy boundary is applied alone. The `cum.reject` value on the overall row is the power of that endpoint considered on its own. ```{r oc-pfs, eval = has_gsd} oc_PFS <- simsummary_fast( pfs_res, eff.col = "logrank.z", efficacy = eff_PFS, fut.col = "logrank.z", futility = fut_PFS, direction = "lower" ) oc_PFS ``` ```{r oc-os, eval = has_gsd} oc_OS <- simsummary_fast( os_res, eff.col = "logrank.z", efficacy = eff_OS, direction = "lower" ) oc_OS ``` Under the alternative the PFS futility rule rarely stops the trial, which is the intended behavior: a futility look is meant to protect against continuing an ineffective treatment, not to interrupt a genuinely effective one. ### Sequential testing of PFS then OS The procedure tests OS only after PFS has been declared significant, and OS may be claimed only at the same look as, or a later look than, the PFS claim. The per-simulation rejection looks are recovered from the boundary-crossing logic, and the joint power is the probability of declaring both endpoints under this gatekeeping rule. ```{r joint, eval = has_gsd} # First efficacy-crossing look for a matrix of statistics, honouring an optional # futility rule that stops the trial without a rejection. reject_look <- function(Z, efficacy, futility = NULL) { nl <- ncol(Z) ns <- nrow(Z) stop_flag <- logical(ns) out <- rep(NA_integer_, ns) for (l in seq_len(nl)) { active <- !stop_flag if (!is.na(efficacy[l])) { hit <- active & !is.na(Z[, l]) & Z[, l] <= efficacy[l] out[hit] <- l stop_flag[hit] <- TRUE } if (!is.null(futility) && !is.na(futility[l])) { fut_hit <- !stop_flag & active & !is.na(Z[, l]) & Z[, l] >= futility[l] stop_flag[fut_hit] <- TRUE } } out } pfs_look <- reject_look(Z_PFS, eff_PFS, fut_PFS) os_look <- reject_look(Z_OS, eff_OS) pfs_reject <- !is.na(pfs_look) os_reject <- !is.na(os_look) # Hierarchical gatekeeping: OS is claimed only if PFS is claimed at the same or # an earlier look. joint_reject <- pfs_reject & os_reject & (os_look >= pfs_look) power_table <- data.frame( quantity = c("PFS power (marginal)", "OS power (marginal)", "OS power (after PFS, hierarchical)", "Joint power (both endpoints)"), value = c(mean(pfs_reject), mean(os_reject), mean(joint_reject), mean(joint_reject)) ) power_table$value <- round(power_table$value, 3) power_table ``` The marginal OS power is the probability of crossing the OS efficacy boundary ignoring the gatekeeping, and the hierarchical OS power is the probability of claiming OS within the procedure. The latter cannot exceed the PFS power, since OS is reachable only through a PFS rejection. ## Correlation of the log-rank statistics The standardized log-rank statistics for the two endpoints are positively correlated because they share the OS component and are evaluated on overlapping risk sets at the same calendar cutoffs. The full empirical correlation matrix of `Z = (Z_PFS_1, Z_PFS_2, Z_PFS_3, Z_OS_1, Z_OS_2, Z_OS_3)` summarizes both the within-endpoint correlation across looks and the cross-endpoint correlation. ```{r correlation} Zc <- cbind(Z_PFS, Z_OS) colnames(Zc) <- c(paste0("PFS_", seq_len(L)), paste0("OS_", seq_len(L))) cor_mat <- cor(Zc, use = "pairwise.complete.obs") round(cor_mat, 3) ``` The cross-endpoint correlation at matching looks is the diagonal of the PFS-by-OS block. ```{r cross-endpoint} cross_block <- cor_mat[paste0("PFS_", seq_len(L)), paste0("OS_", seq_len(L)), drop = FALSE] data.frame( look = seq_len(L), cor_PFS_OS = round(diag(cross_block), 3) ) ``` The within-endpoint correlation across looks follows the canonical group-sequential form, in which the correlation between two looks of the same endpoint is the square root of the ratio of their event counts, `sqrt(d_min / d_max)`. The empirical values reproduce this closed form, with the PFS event counts equal to the exact targets and the OS event counts taken from the simulation. ```{r within-endpoint} pairs_idx <- rbind(c(1, 2), c(1, 3), c(2, 3)) within_tab <- function(Zmat, d_counts, tag) { emp <- apply(pairs_idx, 1, function(p) { stats::cor(Zmat[, p[1]], Zmat[, p[2]], use = "complete.obs") }) theo <- apply(pairs_idx, 1, function(p) { sqrt(min(d_counts[p]) / max(d_counts[p])) }) data.frame( endpoint = tag, look_pair = paste0(pairs_idx[, 1], "-", pairs_idx[, 2]), empirical = round(emp, 3), closed_form = round(theo, 3) ) } rbind( within_tab(Z_PFS, d_PFS, "PFS"), within_tab(Z_OS, d_OS_emp, "OS") ) ``` A scatter of the two final-look statistics shows the positive association directly. Each point is one simulated trial; the cloud is tilted, reflecting the shared OS component. ```{r scatter, fig.cap = "Standardized log-rank statistics at the final look."} plot( Z_PFS[, L], Z_OS[, L], pch = 16, col = grDevices::adjustcolor("steelblue", alpha.f = 0.25), xlab = expression(Z[PFS]), ylab = expression(Z[OS]), main = "Final-look log-rank statistics" ) abline(h = 0, v = 0, col = "grey70", lty = 3) ``` ## Remarks The cross-endpoint correlation at matching looks is moderate, around 0.6, and is largest at the first look, easing slightly as the trial matures. This cross-endpoint correlation is the input a sequential PFS-then-OS procedure needs to characterize its joint operating characteristics, and the value estimated here by simulation matches the closed-form expression derived for the Fleischer model in the companion methodological note. The within-endpoint correlations match the canonical group-sequential form to simulation error, which serves as a check on the event-driven timing and the shared-cutoff OS analysis. The simulation uses the validated Rcpp primitives behind `analysis_fast` for every log-rank computation, so the entire study of `r format(nsim, big.mark = ",")` trials with three looks per endpoint runs in a few seconds. Raising `nsim` tightens the correlation and power estimates at a proportional cost. ### References Fleischer, F., Gaschler-Markefski, B., and Bluhmki, E. (2009). A statistical model for the dependence between progression-free survival and overall survival. *Statistics in Medicine*, 28(21), 2669-2686. Lan, K. K. G. and DeMets, D. L. (1983). Discrete sequential boundaries for clinical trials. *Biometrika*, 70(3), 659-663. O'Brien, P. C. and Fleming, T. R. (1979). A multiple testing procedure for clinical trials. *Biometrics*, 35(3), 549-556. Pocock, S. J. (1977). Group sequential methods in the design and analysis of clinical trials. *Biometrika*, 64(2), 191-199.