## Supplementary R code for the "Letter to the Editor Regarding the Paper `Step-Down Parametric Procedures for Testing ## Correlated Endpoints in a Group-Sequential Trial' by Yiyong Fu" ## ## Last change: 18June2018 library(mvtnorm) library(gsDesign) #--- Calculation of group-sequential boundaries for the intersection test of the example in Fu (2018, Section 6.1) # Parameters used in the example w1 <- 0.4 w2 <- 1-w1 size_sg <- 0.7 alpha0 <− 0.025 interim_timing <- c(0.5,1) # functions provided by Fu (2018, appendix) (modified to use algorithm Miwa with maximal precision) dunn <- function(c,w1,w2, alpha, rho){ mean <- rep(0, 2) corr <- diag(2) corr[lower.tri(corr)] <- rho corr[upper.tri(corr)] <- rho z <- pmvnorm(lower = -Inf , upper = c(qnorm(1-alpha*w1*c), qnorm(1-alpha*w2*c)), mean = mean, corr = corr,algorithm =Miwa(steps=4096)) nonlin <- 1 - z - alpha return(nonlin) } y2 <- function(w1,w2,alpha,rho){ z <- uniroot(dunn, lower = 1, upper = 9,w1 = w1,w2 = w2,alpha=alpha,rho = rho)$root return(z) } # calculate upscaling factor xsi and corresponding boundaries as per Fu (2018) ( xsi <- y2(w1 = w1,w2=w2, alpha = 0.025, rho = sqrt(size_sg)) ) # O'Brien-Fleming boundary for the intersection test with this upscaling as per Fu (2018) b_1 <- gsDesign(k = length(interim_timing), test.type = 1,timing=interim_timing,sfu="OF",alpha = w1*xsi*alpha0) b_2 <- gsDesign(k = length(interim_timing), test.type = 1,sfu="OF",timing=interim_timing,alpha = w2*xsi*alpha0) # Extra note: Boundary for hypothesis test of individual hypotheses should be at level alpha=alpha0 and NOT xsi*alpha0 as in Fu 2018 (appendix) # # Pocock boundary for the intersection test with this upscaling as per Fu (2018) # b_1 <- gsDesign(k = length(interim_timing), test.type = 1,timing=interim_timing,sfu="Pocock",alpha = w1*xsi*alpha0) # b_2 <- gsDesign(k = length(interim_timing), test.type = 1,sfu="Pocock",timing=interim_timing,alpha = w2*xsi*alpha0) # corresponding critical nominal levels on Z and P scale z_crit_1 <- b_1$upper$bound z_crit_2 <- b_2$upper$bound ( p_crit_1 <- pnorm(-z_crit_1) ) ( p_crit_2 <- pnorm(-z_crit_2) ) ( p_ratio <- p_crit_1/p_crit_2 ) ##--- Calculate exact type 1 error under the globa null for this boundary # calculate full correlation structure between endpoints and interim analysis cor.u <- matrix(c(1,sqrt(size_sg),sqrt(size_sg),1),nrow=2) # correlation between subgroup and full population cor.t <- outer(interim_timing,interim_timing,function(x,z) sqrt(pmin(x,z)/pmax(x,z))) # correlation between interim data at 50% and final data tot.cor <- cor.t %x% cor.u # calculate type 1 error under the global null using algorithm Miwa with high accuracy 1-pmvnorm(lower=rep(-Inf,nrow(tot.cor)), upper=c(rbind(z_crit_1,z_crit_2)), corr=tot.cor,algorithm =Miwa(steps=4096)) # calculate type 1 error under the global null again using algorithm GenzBretz with high accuracy # (takes a while...) set.seed(1) # for reproducibility 1-pmvnorm(lower=rep(-Inf,nrow(tot.cor)), upper=c(rbind(z_crit_1,z_crit_2)), corr=tot.cor,algorithm =GenzBretz(maxpts = 1E8, abseps = 1E-12, releps = 0))