R/implementation_functions.R

In gsearly: Creates Group Sequential Trial Designs when Early Outcomes are Available

#-------------------------------------------------------------------------------
# Implementation functions for gsearly
#-------------------------------------------------------------------------------
# 1st July 2026
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 1.  expectSampsize
#-------------------------------------------------------------------------------
expectSampsize <- function(mod, signif = 3) {
  ## Using lower and upper stopping probabilities
  ninttotal <- c(mod$rdata$intnumbers[, 1], mod$rdata$n["total"])
  nintcontl <- mod$rdata$vphi * ninttotal
  ninttreat <- (1 - mod$rdata$vphi) * ninttotal
  ncontl <- mod$power$lowerror["prob", ] %*% nintcontl + mod$power$upperror["prob",
    ] %*% nintcontl
  ntreat <- mod$power$lowerror["prob", ] %*% ninttreat + mod$power$upperror["prob",
    ] %*% ninttreat
  ntotal <- ncontl + ntreat
  n <- c(control = round(as.numeric(ncontl), signif), treat = round(as.numeric(ntreat),
    signif), total = round(as.numeric(ntotal), signif))
  names(n) <- names(mod$rdata$n)
  return(n)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 6.  fixedSampsize
#-------------------------------------------------------------------------------
fixedSampsize <- function(mod, pow = NULL, roundup = TRUE, signif = 3) {

  ## Set power or take from original design
  if (is.null(pow) == TRUE) {
    if (is.na(mod$power$setpow) == FALSE) {
      pow <- mod$power$setpow
    } else {
      pow <- mod$power$totalerror["efficacy"]
    }
  } else {
    if (pow <= 0 | pow >= 1) {
      stop("power: pow must be in range (0,1)")
    }
  }
  nlooks <- mod$power$nlooks
  theta <- mod$power$theta
  contrat <- mod$rdata$contrat
  if (theta <= 0) {
    stop("theta: theta must be > 0")
  }

  ## Calculate sample size
  mid_term <- ((qnorm(1-mod$power$fp[nlooks]) + qnorm(pow))/theta)^2
  n1 <- (1/mod$rdata$vphi) * (mod$idata$sd^2) * mid_term
  n0 <- (mod$rdata$vphi/(1 - mod$rdata$vphi)) * n1
  n <- n0 + n1
  quot <- n %/% contrat[2]
  if(roundup ==  FALSE){
   n0 <- as.integer(quot * contrat[1])
   n <- as.integer(quot * contrat[2])
   n1 <- n - n0
  } else {
   n0 <- as.integer((quot + 1) * contrat[1])
   n <- as.integer((quot + 1) * contrat[2])
   n1 <- n - n0
  }
  ## Determine power for rounded design
  roundpow <- pnorm(theta * sqrt(mod$rdata$vphi * n1/(mod$idata$sd^2)) -
    qnorm(1-mod$power$fp[nlooks]), lower.tail = FALSE)
  roundpow <- 1 - roundpow
  n <- c(n0, n1, n = n)
  names(n) <- names(mod$rdata$n)

  ## Output
  return(list(n = n, power = round(roundpow, signif)))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 12.  gsearlyModel
#-------------------------------------------------------------------------------
gsearlyModel <- function(rmodel = "fix", trecruit, s, tfu, tinterims, pow = 0.9,
  n = NULL, theta, tref = c(1, 2), contrat = c(1,2), roundup = TRUE, m = 2, 
  cmodel = "uniform", sd = 1, rho = 0.5, fp, tn, 
  treatnames = c("control", "treat"), sopt = list(r = 18, 
  bisect = list(min = 20, max = 10000, niter = 1000, tol = 0.001))) {

  ## Validate rmodel and design inputs
  .valContrat(contrat)
  vphi <- contrat[1]/contrat[2]
  .valRmodel(rmodel)
  .valRecruitmod(rmodel = rmodel, trecruit = trecruit, s = s, tfu = tfu,
    theta = theta, tinterims = tinterims, tref = tref, vphi = vphi, 
    roundup = roundup, pow = pow, m = m, n = n, sopt = sopt)
  if (!is.null(n) && n%%2 != 0) {
    n <- as.integer(ceiling(n) + 1)
  }
  if (!is.null(n)) {
    n <- as.integer(n)
  }
  tfu <- tfuStandard(tfu = tfu, tref = tref)
  ## Force positive theta
  theta <- abs(theta)

  ## Validate information model inputs
  .valCorrmod(rmodel = rmodel, cmodel = cmodel, sd = sd, rho = rho, s = s)

  ## Information fraction
  tint <- c(tinterims, trecruit + tfu$tfu[s])
  if (cmodel == "uniform") {
    iint <- sapply(tint, .tau, cmodel = "uniform", tfu = tfu, trecruit = trecruit,
      sd = sd, rho = rho, vphi = vphi, rmodel = rmodel, m = m)
    relvar <- sapply(tint, .rvarUnif, tfu = tfu, alpha = rho, trecruit = trecruit,
      rmodel = rmodel, m = m)
    corrmat <- corrUnif(rho = rho, tfu = tfu)
  } else if (cmodel == "exponential") {
    iint <- sapply(tint, .tau, cmodel = "exponential", tfu = tfu, trecruit = trecruit,
      sd = sd, rho = rho, vphi = vphi, rmodel = rmodel, m = m)
    relvar <- sapply(tint, .rvarExp, tfu = tfu, alpha = rho, trecruit = trecruit,
      rmodel = rmodel, m = m)
    corrmat <- corrExp(rho = rho, tfu = tfu)
  }
  colnames(corrmat) <- rownames(corrmat) <- tfu$tfu

  ## Get bounds
  nlooks <- length(iint)
  .valFptn(fp = fp, tn = tn, nlooks = nlooks)
  dfp <- as.numeric(diff(c(0, fp)))
  dtn <- as.numeric(diff(c(0, tn)))
  ggbound <- gsBound(I = iint, trueneg = dtn, falsepos = dfp,
    r = as.numeric(sopt$r))

  ## Power function
  largs <- list(tfu = tfu, trecruit = trecruit, sd = sd, vphi = vphi,
    rmodel = rmodel, m = m, tint = tint, relvar = relvar, theta = theta,
    nlooks = nlooks, ggbound = ggbound, r = as.numeric(sopt$r))
  ffn <- function(x, pow, largs, powonly = TRUE) {
    vb0 <- sapply(largs$tint, .vbeta0, n = x, tfu = largs$tfu, trecruit = largs$trecruit,
      sd = largs$sd, vphi = largs$vphi, rmodel = largs$rmodel, m = largs$m)
    vvb <- vb0 * largs$relvar
    inform <- 1/vvb
    getpower <- gsProbability(k = largs$nlooks, theta = largs$theta,
      n.I = inform, a = largs$ggbound$a, b = largs$ggbound$b, r = largs$r)
    estpow <- sum(getpower$upper$prob)
    if (powonly == TRUE) {
      return(estpow - pow)
    } else {
      return(getpower)
    }
  }

  ## Get n
  if (is.null(n)) {
    getn <- .bisection(ffn, pow = pow, largs = largs, a = as.numeric(sopt$bisect$min),
      b = as.numeric(sopt$bisect$max), niter = as.numeric(sopt$bisect$niter),
      tol = as.numeric(sopt$bisect$tol))
    ## Target power
    tpow <- ffn(getn, pow = 0, largs = largs)
    quot <- getn %/% contrat[2]
    if(roundup ==  FALSE){
     n0 <- as.integer(quot * contrat[1])
     n <- as.integer(quot * contrat[2])
     n1 <- n - n0
    } else {
     n0 <- as.integer((quot + 1) * contrat[1])
     n <- as.integer((quot + 1) * contrat[2])
     n1 <- n - n0
    }

    ## Get rdata and idata and power
    rdata <- .recruitData(rmodel = rmodel, trecruit = trecruit, s = s,
      tfu = tfu, tinterims = tinterims, n = n, vphi = vphi, contrat = contrat, 
      m = m, intonly = TRUE)
    irecruit_mod <- as.integer(match(rmodel, c("dilin", "idlin", "diquad",
      "idquad", "filin", "fdlin", "ilinf", "dlinf"), -1))
    if (irecruit_mod < 1) {
      rdata$m <- NA
    }
    idata <- .informData(rdata = rdata, cmodel = cmodel, sd = sd, rho = rho,
      intonly = TRUE)

    ## Get power
    getpower <- ffn(n, pow = 0, largs = largs, powonly = FALSE)
    lowerror <- matrix(c(getpower[["lower"]]$bound, getpower[["lower"]]$prob),
      nrow = 2, ncol = nlooks, byrow = TRUE)
    upperror <- matrix(c(getpower[["upper"]]$bound, getpower[["upper"]]$prob),
      nrow = 2, ncol = nlooks, byrow = TRUE)
    colnames(lowerror) <- colnames(upperror) <- idata$tlooks
    rownames(lowerror) <- rownames(upperror) <- c("bound", "prob")
    totpow <- c(futility = sum(lowerror["prob", ]), efficacy = sum(upperror["prob",
      ]))
    lgetpower <- list(nlooks = nlooks, theta = theta, setpow = tpow,
      fp = fp, tn = tn, lowerror = lowerror, upperror = upperror,
      totalerror = totpow)

    ## Add group sizes to rdata
    rdata$n <- c(n0, n1, n)
    names(rdata$n) <- c(treatnames, "total")

  } else {

    ## Target power
    tpow <- NA
    ## Group sizes
    quot <- n %/% contrat[2]
    if(roundup ==  FALSE){
     n0 <- as.integer(quot * contrat[1])
     n <- as.integer(quot * contrat[2])
     n1 <- n - n0
    } else {
     n0 <- as.integer((quot + 1) * contrat[1])
     n <- as.integer((quot + 1) * contrat[2])
     n1 <- n - n0
    }

    ## Get rdata and idata and power
    rdata <- .recruitData(rmodel = rmodel, trecruit = trecruit, s = s,
      tfu = tfu, tinterims = tinterims, n = n, vphi = vphi, contrat = contrat,
      m = m, intonly = TRUE)
    idata <- .informData(rdata = rdata, cmodel = cmodel, sd = sd, rho = rho,
      intonly = TRUE)

    ## Get power
    getpower <- ffn(n, pow = 0, largs = largs, powonly = FALSE)
    lowerror <- matrix(c(getpower[["lower"]]$bound, getpower[["lower"]]$prob),
      nrow = 2, ncol = nlooks, byrow = TRUE)
    upperror <- matrix(c(getpower[["upper"]]$bound, getpower[["upper"]]$prob),
      nrow = 2, ncol = nlooks, byrow = TRUE)
    colnames(lowerror) <- colnames(upperror) <- idata$tlooks
    rownames(lowerror) <- rownames(upperror) <- c("bound", "prob")
    totpow <- c(futility = sum(lowerror["prob", ]), efficacy = sum(upperror["prob",
      ]))
    lgetpower <- list(nlooks = nlooks, theta = theta, setpow = tpow,
      fp = fp, tn = tn, lowerror = lowerror, upperror = upperror,
      totalerror = totpow)

    ## Add group sizes to rdata
    rdata$n <- c(n0, n1, n)
    names(rdata$n) <- c(treatnames, "total")

  }
  ## Data for cmodel
  idata$cmodel <- list(type = cmodel, rho = rho, corrmat = corrmat)
  ## Output
  names(rdata) <- c("rmodel", "m", "trecruit", "s", "tfu", "n", "vphi",     
                     "contrat", "tinterims", "intnumbers")
  names(idata) <- c("cmodel", "sd", "tlooks", "inform")
  packdetail <- "gsearly: 2026-07-01 version 1.1.0"
  gsearlydata <- list(title = packdetail, call = match.call(), rdata = rdata,
    idata = idata, power = lgetpower)
  class(gsearlydata) <- "gsearly"
  return(gsearlydata)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 14.  gsearlyUser
#-------------------------------------------------------------------------------
gsearlyUser <- function(trecruit, s, tfu, tinterims, ninterims, pow = 0.9,
  n = NULL, tref = c(1, 2), contrat = c(1,2), roundup = TRUE, 
  cmodel = "uniform", sd = 1, rho = 0.5, theta, fp, tn, 
  treatnames = c("control", "treat"), sopt = list(r = 18, 
  bisect = list(min = 20, max = 10000, niter = 1000, tol = 0.001))) {

  ## Validate rmodel and design inputs
  rmodel <- "none"
  m <- 2
  .valContrat(contrat)
  vphi <- contrat[1]/contrat[2]
  .valRecruitmod(rmodel = rmodel, trecruit = trecruit, s = s, tfu = tfu,
    theta = theta, tinterims = tinterims, ninterims = ninterims, tref = tref,
    vphi = vphi, roundup = roundup, pow = pow, m = m, n = n, sopt = sopt)
  if (!is.null(n) && n%%2 != 0) {
    n <- as.integer(ceiling(n) + 1)
  }
  if (!is.null(n)) {
    n <- as.integer(n)
  }
  tfu <- tfuStandard(tfu = tfu, tref = tref)
  ## Force positive theta
  theta <- abs(theta)

  ## Validate information model inputs
  .valCorrmod(rmodel = rmodel, cmodel = cmodel, sd = sd, rho = rho, s = s)

  ## Get covariance matrix
  if (is.matrix(cmodel) == TRUE) {
    corrmat <- cmodel
  } else {
    if (cmodel == "uniform") {
      corrmat <- corrUnif(rho = rho, tfu = tfu)
    } else if (cmodel == "exponential") {
      corrmat <- corrExp(rho = rho, tfu = tfu)
    }
  }
  colnames(corrmat) <- rownames(corrmat) <- tfu$tfu
  sdmat <- diag(rep(sd, s))
  vcovmat <- sdmat %*% corrmat %*% sdmat

  ## Validate fp and tn
  alltfu <- c(0, as.numeric(tfu$tfu))
  nlooks <- length(tinterims) + 1
  .valFptn(fp = fp, tn = tn, nlooks = nlooks)
  dfp <- as.numeric(diff(c(0, fp)))
  dtn <- as.numeric(diff(c(0, tn)))

  ## Power function
  tint <- c(tinterims, trecruit + tfu$tfu[s])
  largs <- list(tfu = tfu, trecruit = trecruit, sd = sd, vphi = vphi,
    tint = tint, ninterims = ninterims, s = s, theta = theta, nlooks = nlooks,
    vcovmat = vcovmat, dfp = dfp, dtn = dtn, r = as.numeric(sopt$r))
  ffn <- function(x, pow, largs, powonly = TRUE) {
    getinf <- .userinformData(x = x, largs = largs)
    ggbound <- gsBound(I = getinf["tau", ], trueneg = largs$dtn,
      falsepos = largs$dfp, r = largs$r)
    getpower <- gsProbability(k = largs$nlooks, theta = largs$theta,
      n.I = 1/getinf["vbeta", ], a = ggbound$a, b = ggbound$b, r = largs$r)
    estpow <- sum(getpower$upper$prob)
    if (powonly == TRUE) {
      return(estpow - pow)
    } else {
      return(getpower)
    }
  }

  ## Number n must be larger than max of ninterims
  if (is.matrix(ninterims) == TRUE) {
    minn <- max(ninterims)
  } else if (is.function(ninterims) == TRUE) {
    minn <- as.numeric(sopt$bisect$min)
  }

  ## Get n
  if (is.null(n)) {

    ## Optimisation can only work if pow is within attainable range
    minpow <- ffn(minn, pow = 0, largs = largs, powonly = TRUE)
    if (pow < minpow) {
      warning("Target power pow < min power; change pow or design",
        immediate. = FALSE)
    }
    if (is.function(ninterims) == FALSE) {
      getmaxpow <- ffn(minn, pow = 0, largs = largs, powonly = FALSE)
      maxpow <- 1 - sum(getmaxpow[["lower"]]$prob[1:(nlooks - 1)])
      if (pow > maxpow) {
        warning("Target power pow > max power; change pow or design",
          immediate. = FALSE)
      }
    } else {
      maxpow <- 1
    }
    ## Set pow to mid range if pow outside range
    rangepow <- c(minpow, maxpow)
    if (pow > maxpow || pow < minpow) {
      pow <- (maxpow + minpow)/2
    }

    ## Find n for target power
    getn <- .bisection(ffn, pow = pow, largs = largs, a = minn, b = as.numeric(sopt$bisect$max),
      niter = as.numeric(sopt$bisect$niter), tol = as.numeric(sopt$bisect$tol))

    ## Target power and n
    tpow <- ffn(getn, pow = 0, largs = largs)
    ## Actual power
    quot <- getn %/% contrat[2]
    if(roundup ==  FALSE){
     n0 <- as.integer(quot * contrat[1])
     n <- as.integer(quot * contrat[2])
     n1 <- n - n0
    } else {
     n0 <- as.integer((quot + 1) * contrat[1])
     n <- as.integer((quot + 1) * contrat[2])
     n1 <- n - n0
    }
    if (n == as.numeric(sopt$bisect$max)) {
      warning("Optimisation has hit max n; increase sopt$bisect$max",
        immediate. = FALSE)
    }

    ## Get rdata and idata and power
    if (is.matrix(ninterims) == TRUE) {
      rdata <- list(rmodel = rmodel, m = NA, trecruit = trecruit, s = s,
        tfu = tfu, n = n, vphi = vphi, contrat = contrat, tinterims = tinterims,
        interims = ninterims)
    } else if (is.function(ninterims) == TRUE) {
      rdata <- list(rmodel = rmodel, m = NA, trecruit = trecruit, s = s,
        tfu = tfu, n = n, vphi = vphi, contrat = contrat, tinterims = tinterims,
        interims = ninterims(n))
    }
    idata <- list(cmodel = cmodel, sd = sd, tinterims = tint, interims = .userinformData(n,
      largs = largs))
    colnames(idata$interims) <- as.character(tint)

    ## Get power
    getpower <- ffn(n, pow = 0, largs = largs, powonly = FALSE)
    lowerror <- matrix(c(getpower[["lower"]]$bound, getpower[["lower"]]$prob),
      nrow = 2, ncol = nlooks, byrow = TRUE)
    upperror <- matrix(c(getpower[["upper"]]$bound, getpower[["upper"]]$prob),
      nrow = 2, ncol = nlooks, byrow = TRUE)
    colnames(lowerror) <- colnames(upperror) <- colnames(idata$interims)
    rownames(lowerror) <- rownames(upperror) <- c("bound", "prob")
    totpow <- c(futility = sum(lowerror["prob", ]), efficacy = sum(upperror["prob",
      ]))
    lgetpower <- list(nlooks = nlooks, theta = theta, setpow = tpow,
      rangepow = rangepow, fp = fp, tn = tn, lowerror = lowerror,
      upperror = upperror, totalerror = totpow)

  } else {

    ## Target power
    tpow <- NA
    ## Group sizes
    if (n < minn) {
      stop("Need n>=max(niterims) to evaluate design")
    }
    quot <- n %/% contrat[2]
    if(roundup ==  FALSE){
     n0 <- as.integer(quot * contrat[1])
     n <- as.integer(quot * contrat[2])
     n1 <- n - n0
    } else {
     n0 <- as.integer((quot + 1) * contrat[1])
     n <- as.integer((quot + 1) * contrat[2])
     n1 <- n - n0
    }

    ## Get rdata and idata and power
    if (is.matrix(ninterims) == TRUE) {
      rdata <- list(rmodel = rmodel, m = NA, trecruit = trecruit, s = s,
        tfu = tfu, n = n, vphi = vphi, contrat = contrat, tinterims = tinterims,
        interims = ninterims)
    } else if (is.function(ninterims) == TRUE) {
      rdata <- list(rmodel = rmodel, m = NA, trecruit = trecruit, s = s,
        tfu = tfu, n = n, vphi = vphi, contrat = contrat, tinterims = tinterims,
        interims = ninterims(n))
    }
    idata <- list(cmodel = cmodel, sd = sd, tinterims = tint, interims = .userinformData(n,
      largs = largs))
    colnames(idata$interims) <- as.character(tint)

    ## Get power
    getpower <- ffn(n, pow = 0, largs = largs, powonly = FALSE)
    lowerror <- matrix(c(getpower[["lower"]]$bound, getpower[["lower"]]$prob),
      nrow = 2, ncol = nlooks, byrow = TRUE)
    upperror <- matrix(c(getpower[["upper"]]$bound, getpower[["upper"]]$prob),
      nrow = 2, ncol = nlooks, byrow = TRUE)
    colnames(lowerror) <- colnames(upperror) <- colnames(idata$interims)
    rownames(lowerror) <- rownames(upperror) <- c("bound", "prob")
    totpow <- c(futility = sum(lowerror["prob", ]), efficacy = sum(upperror["prob",
      ]))
    lgetpower <- list(nlooks = nlooks, theta = theta, setpow = tpow,
      fp = fp, tn = tn, lowerror = lowerror, upperror = upperror,
      totalerror = totpow)
  }

  ## Data for rmodel
  rdata$n <- c(n0, n1, n)
  names(rdata$n) <- c(treatnames, "total")
  colnames(rdata$interims) <- as.character(alltfu)
  rownames(rdata$interims) <- as.character(tinterims)

  ## Data for cmodel
  if (is.matrix(cmodel) == TRUE) {
    idata$cmodel <- list(type = "none", rho = NA, corrmat = corrmat)
  } else {
    idata$cmodel <- list(type = cmodel, rho = rho, corrmat = corrmat)
  }

  ## Output
  names(rdata) <- c("rmodel", "m", "trecruit", "s", "tfu", "n", "vphi",     
                     "contrat", "tinterims", "intnumbers")
  names(idata) <- c("cmodel", "sd", "tlooks", "inform")
  packdetail <- "gsearly: 2026-07-01 version 1.1.0"
  gsearlydata <- list(title = packdetail, call = match.call(), rdata = rdata,
    idata = idata, power = lgetpower)
  class(gsearlydata) <- "gsearly"
  return(gsearlydata)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 26.  roundInterims
#-------------------------------------------------------------------------------
roundInterims <- function(mod, roundup = TRUE, full = FALSE) {

  ## Set-up variables
  x <- mod$rdata$intnumbers
  x <- rbind(x, rep(as.integer(mod$rdata$n["total"]), 
      times = mod$rdata$s + 1))
  rownames(x) <- mod$idata$tlooks
  nlooks <- mod$power$nlooks
  s <- mod$rdata$s
  vphi <- mod$rdata$vphi
  contrat <- mod$rdata$contrat
  quotx <- x %/% contrat[2]
  if(roundup ==  FALSE){
   n0 <- quotx * contrat[1]; mode(n0) <- 'integer'
   n <- quotx * contrat[2]; mode(n1) <- 'integer'
   n1 <- n - n0
  } else {
   n0 <- (quotx + 1) * contrat[1]; mode(n0) <- 'integer'
   n <- (quotx + 1) * contrat[2]; mode(n) <- 'integer'
   n1 <- n - n0; mode(n1) <- 'integer'
  }

  ## Output
  if(full == TRUE){
   out <- list(n0, n1, n)
   names(out)<- names(mod$rdata$n)
  } else {
   out <- list(n = n)
   names(out)<- tail(names(mod$rdata$n),n=1)
  }
  return(out)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# end
#-------------------------------------------------------------------------------