R/op_power_surv.R

In gsMAMS: Group Sequential Designs of Multi-Arm Multi-Stage Trials

#' @title Provides operating characteristics of group sequential MAMS trial for survival outcome
#' @description Computes power and other characteristics for group-sequential MAMS trial for survival outcome.
#' @param m0 numeric Median survival time of control group.
#' @param alpha numeric Type I error.
#' @param beta numeric Type II error.
#' @param p numeric Number of treatment arms.
#' @param frac numeric Vector of fractions for information time at each look.
#' @param hr0 numeric Hazard ratio of ineffective treatment group vs control.
#' @param hr1 numeric Hazard ratio of effective treatment group vs control.
#' @param nsim numeric Number of simulations.
#' @param ta numeric Accrual time.
#' @param tf numeric Follow-up time.
#' @param kappa numeric Shape parameter (kappa=1 for exponential distribution).
#' @param eta numeric  Rate of loss to follow-up.
#' @param seed numeric Random seed number.
#' @return A list of power, stage-wise probability of success, stopping probability, probability of futility, average number of events happened per arm, average duration of trial.
#' @examples
#' op_power_surv(m0 = 20,
#'               alpha = 0.05,
#'               beta = 0.1,
#'               p = 4,
#'               frac = c(1 / 2, 1),
#'               hr0 = 1,
#'               hr1 = 0.74,
#'               ta = 12,
#'               tf = 40,
#'               nsim = 20,
#'               kappa = 1,
#'               eta = 0,
#'               seed = 12)
#' @export




op_power_surv <- function(m0, alpha, beta, p, frac, hr0, hr1, nsim, ta, tf, kappa, eta, seed) {
  K<-p
  if (K <= 1) {
    stop("K should be greater than 1.")
  }
  if (length(frac) == 1) {
    stop("The length of frac should be greater than 1.")
  }
  j <- length(frac)

  bound <- scprt(alpha = alpha, k = K, frac = frac)

  lambda0 <- log(2) / m0^kappa

  lambda <- numeric()

  lambda[1] <- lambda0 * hr1

  for (i in 2:K) {
    lambda[i] <- lambda0 * hr0
  }

  scale <- vector()
  scale[1] <- 1 / lambda0^(1 / kappa)


  for (i in 2:(K + 1)) {
    scale[i] <- 1 / lambda[i - 1]^(1 / kappa)
  }


  s2 <- numeric(length = j)




  n <- size_surv(m0 = m0, hr0 = hr0, hr1 = hr1, ta = ta, tf = tf, k = K, beta = beta, alpha = alpha, kappa = kappa, eta = eta, frac = frac)[2]
  d <- size_surv(m0 = m0, hr0 = hr0, hr1 = hr1, ta = ta, tf = tf, k = K, beta = beta, alpha = alpha, kappa = kappa, eta = eta, frac = frac)[1]

  frac <- frac
  d1 <- numeric(length = j)
  for (i in 1:j) {
    d1[i] <- ceiling(d * frac[i])
  }



  a <- data.frame()
  for (i in 1:K) {
    a <- rbind.data.frame(a, bound$lshape)
  }

  b <- data.frame()
  for (i in 1:K) {
    b <- rbind.data.frame(b, bound$ushape)
  }





  tstar <- 1 / j



  tau <- ta + tf
  sp <- numeric()
  pf <- numeric()
  smn <- numeric()
  smd <- numeric()
  dur <- numeric()

   
  set.seed(seed)
  for (e in 1:nsim) {
    Q <- rep(0, j)
     
    datagen <- NULL
    for (i in 1:(K + 1)) {
      w <- stats::rweibull(n, kappa, scale[i])
      u <- stats::runif(n, 0, ta) ## generate accrual time
      if (eta != 0) {
        g <- stats::rexp(n, rate = eta)
      }
      if (eta == 0) {
        g <- tau - u
      }
      x <- pmax(0, pmin(w, tau - u, g)) ## observed survival time
      cens <- as.numeric(w < pmin((tau - u), g)) ## censoring indicator
      group <- rep((i - 1), n)
      datagen[[i]] <- cbind(x, cens, group, u)
    }


    data <- NULL
    for (i in 1:K) {
      dat <- rbind.data.frame(datagen[[1]], datagen[[i + 1]])
      colnames(dat) <- c("time", "event", "group", "accrualtime")
       
      dat$calendar_T <- dat$accrualtime + dat$time
      dat <- dat[order(dat$calendar_T), ]
      dat$cum_event <- cumsum(dat$event)
      data[[i]] <- dat
       
    }


    loc <- matrix(NA, K, j)
    z <- matrix(NA, K, j)
    var_result <- matrix(NA, K, j)
    N_result <- matrix(NA, K + 1, j)
    d_result <- matrix(NA, K + 1, j)
    calendarT_look <- matrix(NA, K, j)
    stagensizeT <- matrix(NA, K, j)
    stagensizeC <- matrix(NA, K, j)
    stagedsizeT <- matrix(NA, K, j)
    stagedsizeC <- matrix(NA, K, j)


    for (h in (1:j)) {
      for (k in (1:K)) {
        if (h < j) {
          dlook <- ceiling(tstar * h * 2 * d)
          loc <- min(which(data[[k]]$cum_event == dlook))
          calendarT_look[k, h] <- data[[k]]$calendar_T[loc]
          data_compare_part1 <- data[[k]][1:loc, ]
          data_compare_part2 <- data[[k]][(loc + 1):(2 * n), ]
          data_compare_part2 <- data_compare_part2[
            data_compare_part2$accrualtime < data_compare_part1$calendar_T[loc],
          ]
          data_compare_part2$event <- 0
          data_compare_part2$time <- data_compare_part1$calendar_T[loc] - data_compare_part2$accrualtime
           
          data_compare <- rbind(data_compare_part1, data_compare_part2)
          stagensizeT[k, h] <- nrow(data_compare[data_compare$group != 0, ])
          stagensizeC[k, h] <- nrow(data_compare[data_compare$group == 0, ])

          stagedsizeT[k, h] <- nrow(data_compare[data_compare$group != 0 & data_compare$event == 1, ])
          stagedsizeC[k, h] <- nrow(data_compare[data_compare$group == 0 & data_compare$event == 1, ])
        }

        
        if (h == j) {
          calendarT_look[k, h] <- tau
          data_compare <- data[[k]]
          stagensizeT[k, h] <- nrow(data_compare[data_compare$group != 0, ])
          stagensizeC[k, h] <- nrow(data_compare[data_compare$group == 0, ])
          stagedsizeT[k, h] <- nrow(data_compare[data_compare$group != 0 & data_compare$event == 1, ])
          stagedsizeC[k, h] <- nrow(data_compare[data_compare$group == 0 & data_compare$event == 1, ])
        }

        temp <- survival::survdiff(survival::Surv(time, event) ~ group, data = data_compare)
        z[k, h] <- sign(temp$obs[1] - temp$exp[1]) * sqrt(temp$chisq)
         
      }
    }


    m1 <- as.numeric()
    for (i in 1:K) {
      if (i == 1) {
        m1[i] <- z[1, 1] > b[i, 1]
      } else {
        m1[i] <- z[1, 1] > z[i, 1]
      }
    }
    pk <- prod(m1)
    if (pk >= 1) {
      s2[1] <- s2[1] + 1
    }


    w <- K
    g <- data.frame(matrix(ncol = w - 1, nrow = 0))

    for (q in 2:(length(g) + 1)) {
      # j<-3
      p <- numeric(length = (j - 1) * 3)
      k <- seq(1, 25, by = 3)[1:(j - 1)]
      sk <- list(a = 1, b = c(2, 3))
       

      for (i in 2:(j)) {
        if (i == 2) {
           
          p[k[i - 1]] <- z[q, (i - 1)] < a[q, (i - 1)]
          p[k[i - 1] + 1] <- z[q, (i - 1)] > a[q, (i - 1)] & z[q, (i - 1)] < b[q, (i - 1)]
          p[k[i - 1] + 2] <- z[1, (i)] > z[q, (i)]

          g[i, q - 1] <- 0
          for (o in seq_len(length(sk))) {
            g[i, q - 1] <- g[i, q - 1] + prod(p[sk[[o]]])
          }
          next
        }

        tk <- sk[[length(sk)]]

        sk[[length(sk)]][length(sk[[length(sk)]])] <- sk[[length(sk)]][length(sk[[length(sk)]])] + 1

        sk[[length(sk) + 1]] <- tk
        sk[[length(sk)]][length(sk[[length(sk)]])] <- sk[[length(sk)]][length(sk[[length(sk)]])] + 2
        sk[[length(sk)]][length(sk[[length(sk)]]) + 1] <- sk[[length(sk)]][length(sk[[length(sk)]])] + 1

        p[k[i - 1]] <- z[q, (i - 1)] < a[q, (i - 1)]
        p[k[i - 1] + 1] <- z[q, (i - 1)] > a[q, (i - 1)] & z[q, (i - 1)] < b[q, (i - 1)]
        p[k[i - 1] + 2] <- z[1, i] > z[q, i]

        g[i, q - 1] <- 0
        for (o in seq_len(length(sk))) {
          g[i, q - 1] <- g[i, q - 1] + prod(p[sk[[o]]])
        }
      }
    }


    sj <- data.frame(matrix(ncol = 1, nrow = 0))

    for (q in 1:1) {
       
      p <- numeric(length = (j - 1) * 2)
      k <- seq(1, 20, by = 2)[1:(j - 1)]
      sk <- list(a = c(1, 2))
       

      for (i in 2:(j)) {
        if (i == 2) {
          p[k[i - 1]] <- z[q, (i - 1)] > a[q, (i - 1)] & z[q, (i - 1)] < b[q, (i - 1)]
          p[k[i - 1] + 1] <- z[1, (i)] > b[q, (i)]

          sj[i, 1] <- prod(p[sk[[length(sk)]]])
          next
        }

         

        sk[[length(sk) + 1]] <- sk[[length(sk)]]
        sk[[length(sk)]][length(sk[[length(sk)]])] <- sk[[length(sk)]][length(sk[[length(sk)]])] + 1
        sk[[length(sk)]][length(sk[[length(sk)]]) + 1] <- sk[[length(sk)]][length(sk[[length(sk)]])] + 1

        p[k[i - 1]] <- z[q, (i - 1)] > a[q, (i - 1)] & z[q, (i - 1)] < b[q, (i - 1)]
        p[k[i - 1] + 1] <- z[1, (i)] > b[q, (i)]


        sj[i, 1] <- prod(p[sk[[length(sk)]]])
      }
    }






    ff <- cbind.data.frame(g, sj)

    pl <- numeric(length = j)

    for (y in 2:j) {
      if (prod(ff[y, ]) >= 1) {
        s2[y] <- s2[y] + 1
      }
    }
 
    stopprob <- rep(0, j)
    probfut <- rep(0, j)
    stop <- rep(0, j)
    stopF <- matrix(NA, K, j)
    stopE <- matrix(NA, K, j)
    Fflag <- rep(0, j)
    Eflag <- rep(0, j)

    for (k in 1:K) {
      for (h in 1:j) {
        stopF[k, h] <- z[k, h] <= a[k, h]
        stopE[k, h] <- z[k, h] >= b[k, h]
      }
    }


    for (k in 1:K) {
      if (any(stopF[k, ] == TRUE)) {
        s <- min(which(stopF[k, ] == TRUE))
        stopF[k, s:j] <- TRUE
      }
    }


    for (h in 1:j) {
      if (sum(stopF[, h] == TRUE) == K) {
        Fflag[h] <- 1
      }
      if (sum(stopE[, h] == TRUE) >= 1) {
        Eflag[h] <- 1
      }
      stop[h] <- Eflag[h] + Fflag[h] # stop is rep (0,J) originally, if we stop at stage 3, we will have (0,0,1,0), thus stop records at which stage we stop for current simulation
    }



    stopstage <- min(which(stop >= 1)) # the stage we should stop, no matter it is due to efficacy or futility
    stopFstage <- ifelse(all(Fflag == 0), 0, min(which(Fflag == 1))) # the stage we  stop due to futility

    if (stopFstage == stopstage) {
      probfut[stopFstage] <- 1
    }

    stopprob[stopstage] <- 1

    samplesize <- rep(0, K + 1)
    samplesized <- rep(0, K + 1)

    for (k in 1:K) {
      stopFstage1 <- ifelse(any(stopF[k, ] == TRUE), min(which(stopF[k, ] == TRUE)), Inf)

      stopstagefork <- min(stopstage, stopFstage1)
      samplesize[k] <- stagensizeT[k, stopstagefork]
      samplesized[k] <- stagedsizeT[k, stopstagefork]
    }
    samplesize[K + 1] <- max(stagensizeC[, stopstage])
    samplesized[K + 1] <- max(stagedsizeC[, stopstage])



    duration <- max(calendarT_look[, stopstage])




    sp <- rbind.data.frame(stopprob, sp)
    pf <- rbind.data.frame(probfut, pf)
    smn[e] <- sum(samplesize) / (K + 1)
    smd[e] <- sum(samplesized) / (K + 1)
    dur[e] <- duration
  }

  power <- round(sum(s2) / nsim, 3)


  s2 <- rbind.data.frame(s2 / nsim)
  names(s2) <- paste0("look", 1:j)
  names(sp) <- paste0("look", 1:j)
  names(pf) <- paste0("look", 1:j)


  p <- list("Power" = power, "Stagewise Power" = colMeans(s2), "Stopping probability under alternative" = colMeans(sp), "Probability of futility under alternative" = colMeans(pf), "Average number of events happened per arm under alternative" = mean(smd), "Average duration of trial(months)" = mean(dur))
  return(p)
}