R/gsSimulate.R

In keaven/gsDesign: Group Sequential Design

# gsAdaptSim roxy [sinew] ----
#' @importFrom stats qnorm rbinom
# gsAdaptSim function [sinew] ----
gsAdaptSim <- function(SimStage, IniSim, TrialPar, SimPar, cp = 0, thetacp = -100, maxn = 100000, pdeltamin = 0, ...) {
  # simulate 1st k-1 stages of trial
  x <- gsSimulate(nstage = TrialPar$gsx$k - 1, SimStage, IniSim, TrialPar, SimPar)

  # for default,  set conditional power to be used
  if (cp == 0) {
    cp <- 1 - x$gsx$beta
  }

  # compute planned information between final interim and final analysisp
  Inew <- x$gsx$n.I[x$gsx$k] - x$gsx$n.I[x$gsx$k - 1]

  # compute final n per group assuming no adaptation
  nnewe <- ceiling(x$gsx$n.I[x$gsx$k] * x$ratio / (1 + x$ratio) - x$ne)

  if (length(nnewe) == 1) {
    nnewe <- rep(nnewe, TrialPar$nsim)
  }

  nnewc <- ceiling(x$gsx$n.I[x$gsx$k] / (1 + x$ratio) - x$nc)

  if (length(nnewc) == 1) {
    nnewc <- rep(nnewc, TrialPar$nsim)
  }

  # compute bounds for each simulation for z based on data between final interim and final analysis
  bnew <- (x$gsx$upper$bound[x$gsx$k] - x$z * sqrt(x$gsx$n.I[x$gsx$k - 1] / x$gsx$n.I[x$gsx$k])) / sqrt(Inew / x$gsx$n.I[x$gsx$k])
  lnew <- (x$gsx$lower$bound[x$gsx$k] - x$z * sqrt(x$gsx$n.I[x$gsx$k - 1] / x$gsx$n.I[x$gsx$k])) / sqrt(Inew / x$gsx$n.I[x$gsx$k])

  # compute estimated standardized treatment effect
  thetahat <- x$z / sqrt(x$gsx$n.I[x$gsx$k - 1])

  # for default,  set theta parameter for which conditional power is to be adjusted
  # to estimated value at interim
  if (thetacp == -100) {
    thetacp <- thetahat
  }

  # compute sample size to achieve desired conditional power
  ncp <- rep(maxn, TrialPar$nsim) - x$nc - x$ne
  ncp[thetacp > 0] <- as.integer(ceiling(((bnew[thetacp > 0] + stats::qnorm(cp)) / thetacp[thetacp > 0])^2))

  if (maxn > 0) {
    maxn <- maxn - x$nc - x$ne
    if (length(maxn) == 1) {
      maxn <- rep(maxn, x$nsim)
    }
    ncp[ncp > maxn] <- maxn[ncp > maxn]
  }

  # set flag for adapting sample size at final analysis based on conditional power,
  # min desired sample size,  and min observed delta
  flag <- 1 * (x$outcome == 0)
  flag <- flag * (x$xc / x$nc - x$xe / x$ne >= pdeltamin)
  flag <- flag * (ncp > ceiling(Inew))

  # compute updated sample sizes for final analysis
  ncpe <- as.integer(ceiling(ncp * x$ratio / (1 + x$ratio)))
  ncpc <- ncp - ncpe
  nnewc[flag == 1] <- ncpc[flag == 1]
  nnewe[flag == 1] <- ncpe[flag == 1]

  # set up analysis object for final stage
  zero <- rep(0, x$nsim)
  y <- list(xc = zero, xe = zero, nc = nnewc, ne = nnewe, xadd = x$xadd, nadd = x$nadd)

  # simulate data for final stage
  for (j in 1:x$nsim)
  {
    if (x$outcome[j] == 0) {
      y$xc[j] <- stats::rbinom(1, nnewc[j], x$pcsim)
      y$xe[j] <- stats::rbinom(1, nnewe[j], x$pesim)
      x$xc[j] <- x$xc[j] + y$xc[j]
      x$xe[j] <- x$xe[j] + y$xe[j]
    }
    # this is just to avoid 0 divide in zStat later
    else {
      y$xc[j] <- y$nc[j] / 2
      y$xe[j] <- y$ne[j] / 2
    }
  }

  if (length(x$nc) == 1) {
    x$nc <- rep(x$nc, TrialPar$nsim)
  }
  x$nc[x$outcome == 0] <- x$nc[x$outcome == 0] + nnewc[x$outcome == 0]

  if (length(x$ne) == 1) {
    x$ne <- rep(x$ne, TrialPar$nsim)
  }
  x$ne[x$outcome == 0] <- x$ne[x$outcome == 0] + nnewe[x$outcome == 0]

  # compute final z statistics
  y <- x$zStat(y)
  x$y <- y

  # update final outcome
  x$outcome <- x$outcome + (x$outcome == 0) * x$gsx$k * ((y$z >= bnew) - (y$z < lnew))

  x
}

# gsSimulate function [sinew] ----
gsSimulate <- function(nstage = 0, SimStage, IniSim, TrialPar, SimPar) {
  if (nstage == 0) {
    nstage <- TrialPar$gsx$k
  }
  x <- IniSim(TrialPar, SimPar)
  nim1 <- 0
  for (i in 1:nstage)
  {
    x <- SimStage(x$gsx$n.I[i], x)
    x$outcome <- x$outcome + (x$outcome == 0) * i *
      ((x$z >= x$gsx$upper$bound[i]) - (x$z < x$gsx$lower$bound[i]))
  }

  x
}