R/plot_functions.R

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

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

#-------------------------------------------------------------------------------
# 13.  .gsearlyPlots
#-------------------------------------------------------------------------------
.gsearlyPlots <- function(plottype) {
  plots <- list(plotRecruit = c("1", "recruit"), plotInform = c("2",
    "inform"), plotBoundary = c("3", "boundary"), plotPower = c("4",
    "power"))
  plottype <- match.arg(tolower(as.character(plottype)), as.vector(unlist(plots)))
  pickplot <- which(unlist(lapply(plots, function(mod, type) {
    is.element(type, mod)
  }, type = plottype)))
  return(names(plots)[pickplot])
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 18.  plot.gsearly
#-------------------------------------------------------------------------------
plot.gsearly <- function(x, plottype = 1, ...) {
  do.call(.gsearlyPlots(plottype), list(x, ...))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 19.  plotBoundary
#-------------------------------------------------------------------------------
plotBoundary <- function (mod, xlim = c(0,1), ylim = NA, xlab = NA, ylab = NA,
          title = NULL, col = 1, lty = c(3, 3), pch = 3, las = 1, concol = grey(0.9),
          reglab = TRUE, signif = 3, bounds = TRUE, pcol = 1, labpos = c(2, 2), ...)
{
  finform <- mod$idata$inform["tau", ]
  lowb <- lbound <- mod$power$lowerror["bound", ]
  uppb <- ubound <- mod$power$upperror["bound", ]
  tn <- diff(c(0,mod$power$tn))
  fp <- diff(c(0,mod$power$fp))
  lowb[tn == 0] <- NA
  uppb[fp == 0] <- NA
  lbound[tn == 0] <- -.Machine$integer.max
  ubound[fp == 0] <- .Machine$integer.max
  nint <- length(finform)
  if (length(labpos) < 2) {
    labpos <- rep(2, 2)
  }
  if (is.na(xlab) == TRUE) {
    xlab <- "Information"
  }
  if (is.na(ylab) == TRUE) {
    ylab <- "Normal Critical Values"
  }
  if (any(is.na(xlim)) == TRUE) {
    xlim <- c(0, 1)
  }
  if (xlim[2] > 1) {
    xlim[2] <- 1
  }
  if (xlim[1] < 0) {
    xlim[1] <- 0
  }
  if (any(is.na(ylim))) {
    ylim <- c(floor(min(c(lowb, uppb), na.rm = TRUE) - 0.25),
              ceiling(max(c(lowb, uppb), na.rm = TRUE) + 0.25))
  }
  else {
    if (ylim[1] < (-100)) {
      ylim[1] <- -100
    }
    if (ylim[2] > (100)) {
      ylim[2] <- 100
    }
  }
  plot(x = NULL, y = NULL, type = "n", ylim = ylim,
                 xlim = xlim, las = las, xlab = xlab, ylab = ylab, main = title,
                 las = las, yaxs = "i")
  polygon(x = c(finform, rev(finform)),
          y = c(lbound, rev(ubound)), border = NA, col = concol)
  box(lty = 1, col = 1)
  abline(h = 0)
  if (length(lty) == 1) {
    lty <- rep(lty, 2)
  }
  lines(x = finform, y = ubound, lty = lty[1], col = col,
                  ...)
  lines(x = finform, y = lbound, lty = lty[2], col = col,
                  ...)
  if (reglab == TRUE) {
    text(x = c(xlim[1], xlim[1]), y = c(ylim[2], ylim[1]),
        expression(paste("Reject ", H[0], sep = ""),
                   paste("Accept ", H[0], sep = "")), pos = c(2, 4),
                   col = pcol, srt=90,offset=c(0,1))
  }
  points(x = finform, y = ubound, pch = pch, col = pcol,
                   ...)
  points(x = finform, y = lbound, pch = pch, col = pcol,
                   ...)
  if (bounds == TRUE) {
    text(x = finform, y = ubound,
                   as.character(round(ubound, signif)), pos = labpos[1], col = col)
    text(x = finform, y = lbound,
                   as.character(round(lbound, signif)), pos = labpos[2], col = col)
  }
  invisible(list(interims = finform, bounds = list(low = lbound,
                                                   upp = ubound)))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 20.  plotInform
#-------------------------------------------------------------------------------
plotInform <- function(mod, xlim = NA, ylim = NA, xlab = NA, ylab = NA,
  title = NULL, col = 1, lty = 1, pch = 3, las = 1, tlag = c(0, 0),
  wopcol = grey(0.9), intlab = NA, labpos = 1, legpos = "bottomleft",
  ptype = "var", pcol = 1, ...){

  ## Find recruitmant data and information
  if (length(tlag) == 1) {
    tlag <- rep(tlag, 2)
  }
  if (tlag[1] < 0 | tlag[2] < 0) {
    tlag <- c(0, 0)
  }
  if (is.null(labpos) == TRUE) {
    labpos <- NA
  }
  if (mod$rdata$rmodel == "none") {
    rdata <- mod$rdata
    n <- as.numeric(rdata$n["total"])
    idata <- mod$idata
  } else {
    rdata <- .recruitData(tfu = mod$rdata$tfu, trecruit = mod$rdata$trecruit,
      s = mod$rdata$s, vphi = mod$rdata$vphi, contrat = mod$rdata$contrat, 
      n = mod$rdata$n["total"], tinterims = mod$rdata$tinterims, 
      rmodel = mod$rdata$rmodel, m = mod$rdata$m, tlag = tlag)
    n <- as.numeric(rdata$n)
    idata <- .informData(rdata = rdata, cmodel = mod$idata$cmodel$type,
      sd = mod$idata$sd, rho = mod$idata$cmodel$rho)
  }
  s <- rdata$s
  ss <- s + 1
  ninterims <- length(rdata$tinterims)
  trecruit <- rdata$trecruit
  tfu <- rdata$tfu$tfu
  alltfu <- c(0, tfu)
  if (length(tlag) == 1) {
    tlag <- rep(tlag, 2)
  }
  if (tlag[1] < 0 | tlag[2] < 0) {
    tlag <- c(0, 0)
  }
  if (any(is.na(xlim)) == TRUE) {
    xlim <- c(-tlag[1], trecruit + tfu[s] + tlag[2])
  }
  if (is.na(xlab) == TRUE) {
    xlab <- "Time"
  }

  ## Get plot data
  if (ptype == "inform") {

    ## Information data
    if (mod$rdata$rmodel == "none") {
      ymain <- c(0, 0, idata$inform["tau", ], tail(idata$inform["tau",
        ], n = 1))
      yint <- idata$inform["tau", ]
      xmain <- c(-tlag[1], tfu[s], idata$tlooks, tail(idata$tlooks,
        n = 1) + tlag[2])
    } else {
      ymain <- idata$data$inform
      yint <- mod$idata$inform["tau", ]
      xmain <- rdata$data$t
    }
    ypoly <- c(0, 1, 1, 0, 0)
    if (any(is.na(ylim)) == TRUE) {
      ylim <- c(0, 1)
    }
    yleg <- "Information"
    if (is.na(ylab) == TRUE) {
      ylab <- "Information fraction"
    }
  } else if (ptype == "var") {

    ## Variance data
    if (mod$rdata$rmodel == "none") {
      ymain <- c(idata$inform["vbeta", ], tail(idata$inform["vbeta",
        ], n = 1))
      yint <- idata$inform["vbeta", ]
      xmain <- c(idata$tlooks, tail(idata$tlooks, n = 1) +
        tlag[2])
    } else {
      ymain <- idata$data$var
      yint <- mod$idata$inform["vbeta", ]
      xmain <- rdata$data$t
      ymain[xmain < idata$tinterims[1]] <- NA
    }
    lowy <- (mod$idata$sd^2)/(mod$rdata$vphi * (1 - mod$rdata$vphi) *
      n)
    ypoly <- c(sqrt(lowy), sqrt(max(ymain, na.rm = TRUE)), sqrt(max(ymain,
      na.rm = TRUE)), sqrt(lowy), sqrt(lowy))
    if (any(is.na(ylim)) == TRUE) {
      ylim <- c(sqrt(lowy), sqrt(max(ymain, na.rm = TRUE)))
    }
    yleg <- "SE"
    if (is.na(ylab) == TRUE) {
      ylab <- "SE"
    }
  }

  # Plot
  if (ptype == "var") {

    ## Variance plot
    plot(x = xmain, y = sqrt(ymain), type = "n", col = 1,
      ylim = ylim, xlim = xlim, las = las, xlab = xlab, ylab = ylab,
      main = title)
    polygon(x = c(rep(tfu[s], 2), rep(trecruit, 2), tfu[s]),
      y = ypoly, border = NA, col = wopcol)
    if (mod$rdata$rmodel == "none") {
      lines(x = xmain, y = sqrt(ymain), col = col, type = "l",
        lty = lty, ...)
    } else {
      lines(x = xmain, y = sqrt(ymain), col = col, type = "l",
        lty = lty, ...)
    }
    for (i in 1:ninterims) {
      points(x = rdata$tinterims[i], y = sqrt(yint[i]),
        pch = pch, col = pcol, ...)
      if (is.na(labpos) != TRUE) {
        text(x = rdata$tinterims[i], y = sqrt(yint[i]),
          if (any(is.na(intlab) == TRUE)) {
          substitute(t[j], list(j = i))
          } else {
          intlab[i]
          }, pos = labpos)
      }
    }
    if (mod$rdata$rmodel == "none") {
      legend(legpos, legend = c(yleg, "Interims"), lty = c(lty,
        NA), pch = c(NA, pch), col = c(col, pcol), cex = 0.9, bty = "n")
    } else {
      legend(legpos, legend = c(yleg, "Interims"), lty = c(lty,
        NA), pch = c(NA, pch), col = c(col, pcol), cex = 0.9, bty = "n")
    }

  } else if (ptype == "inform") {

    ## Information plot
    plot(x = xmain, y = ymain, type = "n", col = 1, ylim = ylim,
      xlim = xlim, las = las, xlab = xlab, ylab = ylab, main = title)
    polygon(x = c(rep(tfu[s], 2), rep(trecruit, 2), tfu[s]),
      y = ypoly, border = NA, col = wopcol)
    if (mod$rdata$rmodel == "none") {
      lines(x = xmain, y = ymain, col = col, type = "l",
        lty = lty, ...)
    } else {
      lines(x = xmain, y = ymain, col = col, type = "l",
        lty = lty, ...)
    }
    for (i in 1:ninterims) {
      points(x = rdata$tinterims[i], y = yint[i], pch = pch,
        col = pcol, ...)
      if (is.na(labpos) != TRUE) {
        text(x = rdata$tinterims[i], y = yint[i], if (any(is.na(intlab) ==
          TRUE)) {
          substitute(t[j], list(j = i))
        } else {
          intlab[i]
        }, pos = labpos)
      }
    }
    if (mod$rdata$rmodel == "none") {
      legend(legpos, legend = c(yleg, "Interims"), lty = c(lty,
        NA), pch = c(NA, pch), col = c(col, pcol), cex = 0.9, bty = "n")
    } else {
      legend(legpos, legend = c(yleg, "Interims"), lty = c(lty,
        NA), pch = c(NA, pch), col = c(col, pcol), cex = 0.9, bty = "n")
    }
  }
  ## Plot data
  summ_plotdata <- data.frame(xmain, ymain)
  names(summ_plotdata) <- c("t", ptype)
  invisible(list(data = summ_plotdata, interims = yint))

}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 21.  plotPower
#-------------------------------------------------------------------------------
plotPower <- function(mod, xlim = NA, ylim = NA, xlab = NA, ylab = NA,
  title = NULL, col = c(1, 2), lty = NA, pch = 3, las = 1,
  legpos = c("bottomright", "topright"), xtype = "theta",
  delta = seq(0, 1, 0.05), legstudy = NA, ...){

  ## Probability, theta and delta data
  theta <- delta * mod$idata$sd
  ntheta <- length(theta)
  nlooks <- mod$power$nlooks
  low.p <- upp.p <- matrix(NA, ncol = nlooks, nrow = ntheta)
  for (i in 1:ntheta) {
    nmod <- update(mod, theta = theta[i], n = mod$rdata$n["total"])
    low.p[i, ] <- cumsum(nmod$power$lowerror["prob", ])
    upp.p[i, ] <- cumsum(nmod$power$upperror["prob", ])
  }
  colnames(upp.p) <- colnames(low.p) <- colnames(nmod$power$lowerror)
  rownames(upp.p) <- rownames(low.p) <- 1:ntheta
  probs <- list(upp.p, low.p)
  pmodel <- list(cumsum(mod$power$upperror["prob", ]), cumsum(mod$power$lowerror["prob",
    ]))

  ## Select plot x-axis
  if (is.na(ylab) == TRUE) {
    ylab <- "Cumulative boundary crossing probability"
  }
  if (xtype == "delta") {
    xrange <- delta
    xpoint <- rep(mod$power$theta/mod$idata$sd, nlooks)
    if (any(is.na(xlim)) == TRUE) {
      xlim <- c(min(delta), max(delta))
    }
    if (is.na(xlab) == TRUE) {
      xlab <- "Standardized effect size"
    }
  } else {
    xrange <- theta
    xpoint <- rep(mod$power$theta, nlooks)
    if (any(is.na(xlim)) == TRUE) {
      xlim <- c(min(theta), max(theta))
    }
    if (is.na(xlab) == TRUE) {
      xlab <- "Effect size"
    }
  }
  if (length(col) != 2) {
    col <- 1:2
  }
  if (any(is.na(ylim)) == TRUE) {
    ylim <- c(0, 1)
  }

  ## Plot
  if (is.na(legstudy) == TRUE) {
    legstudy <- "study"
  }
  if (length(lty) != nlooks) {
    lty <- 1:nlooks
  }
  plot(x = NULL, y = NULL, type = "n", ylim = ylim, xlim = xlim,
    las = las, xlab = xlab, ylab = ylab, main = title, las = las)
  for (j in 1:nlooks) {
    lines(x = xrange, y = probs[[1]][, j], lty = rev(lty)[j],
      col = col[1], ...)
    lines(x = xrange, y = probs[[2]][, j], lty = rev(lty)[j],
      col = col[2], ...)
  }
  points(x = c(xpoint, xpoint), y = c(pmodel[[1]], pmodel[[2]]),
    pch = pch, col = col[1], ...)
  if (!is.na(legpos[1])) {
    legend(legpos[1], lty = c(NA, lty), legend = c(legstudy,
      rev(colnames(mod$power$lowerror))), title = "Interims", pch = c(pch,
      rep(NA, nlooks)), bty = "n", inset = 0.05, col = 1, text.col = 1,
      cex = 0.9)
  }
  if (!is.na(legpos[2])) {
    legend(legpos[2], lty = 1, legend = c("Upper bound",
      "Lower bound"), title = "Probability", pch = NA, bty = "n",
      inset = 0.05, col = col, text.col = 1, cex = 0.9)
  }

  ## Plot data
  names(probs) <- names(pmodel) <- c("upper", "lower")
  invisible(list(data = list(x = xrange, y = probs), interims = list(x = xpoint,
    y = pmodel)))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 22.  plotRecruit
#-------------------------------------------------------------------------------
plotRecruit <- function(mod, xlim = NA, ylim = NA, xlab = NA, ylab = NA,
  title = NULL, col = NA, lty = c(1, 2), pch = 3, las = 1, wopcol = grey(0.9),
  intlab = NA, labpos = 3, legpos = "topleft", tlag = c(0, 0), pcol = 1, ...){

  ## Find recruitment data
  if (length(tlag) == 1) {
    tlag <- rep(tlag, 2)
  }
  if (tlag[1] < 0 | tlag[2] < 0) {
    tlag <- c(0, 0)
  }
  if (is.null(labpos) == TRUE) {
    labpos <- NA
  }
  if (mod$rdata$rmodel == "none") {
    rdata <- mod$rdata
    n <- as.numeric(rdata$n["total"])
  } else {
    rdata <- .recruitData(tfu = mod$rdata$tfu, trecruit = mod$rdata$trecruit,
      s = mod$rdata$s, vphi = mod$rdata$vphi, contrat = mod$rdata$contrat, 
      n = mod$rdata$n["total"], tinterims = mod$rdata$tinterims, 
      rmodel = mod$rdata$rmodel, m = mod$rdata$m, tlag = tlag)
    n <- as.numeric(rdata$n)
  }
  s <- rdata$s
  ss <- s + 1
  ninterims <- length(rdata$tinterims)
  trecruit <- rdata$trecruit
  tfu <- rdata$tfu$tfu
  alltfu <- c(0, tfu)
  if (any(is.na(xlim)) == TRUE) {
    xlim <- c(-tlag[1], trecruit + tfu[s] + tlag[2])
  }
  if (any(is.na(ylim)) == TRUE) {
    ylim <- c(0, n)
  }
  if (any(is.na(col)) == TRUE) {
    col <- 1:ss
  } else {
    if (length(col) != ss) {
      col <- rep(col, ss)
      col <- col[1:ss]
    }
  }
  pcol <- pcol[1]
  if (is.na(ylab) == TRUE) {
    ylab <- "Number Recruited"
  }
  if (is.na(xlab) == TRUE) {
    xlab <- "Time"
  }

  ## Plot
  summ_plotdata <- data.frame()
  plot(x = rdata$data$t, y = rdata$data$n[[1]], type = "n",
    col = 1, ylim = ylim, xlim = xlim, las = las, xlab = xlab, ylab = ylab,
    main = title)
  windowop <- c(low = tfu[s], upp = trecruit)
  polygon(x = c(rep(tfu[s], 2), rep(trecruit, 2), tfu[s]),
    y = c(0, n, n, 0, 0), border = NA, col = wopcol)
  if (length(lty) == 1) {
    lty <- rep(lty, 2)
  }
  if (mod$rdata$rmodel == "none") {
    userPlotdata <- function(x, rdata, s, alltfu, n, tlag) {
      x.dat <- c(-tlag[1], alltfu[x], rdata$tinterims, trecruit +
        alltfu[x], trecruit + alltfu[s + 1] + tlag[2])
      y.dat <- c(0, 0, as.numeric(rdata$intnumbers[, x]), n, n)
      return(data.frame(x = x.dat, y = y.dat))
    }
    for (i in 1:ss) {
      get_plotdata <- userPlotdata(x = i, rdata = rdata, s = s, alltfu = alltfu,
        n = n, tlag = tlag)
      lines(get_plotdata, type = "l", lty = lty[1], col = col[i],
        ...)
      if (i == 1) {
        summ_plotdata <- get_plotdata
      } else {
        summ_plotdata <- cbind(summ_plotdata, get_plotdata[, 2])
      }
    }
  } else {
    for (i in 1:ss) {
      get_plotdata <- data.frame(x = rdata$data$t, y = rdata$data$n[[i]])
      lines(get_plotdata, type = "l", lty = lty[1], col = col[i],
        ...)
      if (i == 1) {
        summ_plotdata <- get_plotdata
      } else {
        summ_plotdata <- cbind(summ_plotdata, get_plotdata[, 2])
      }
    }
  }
  for (i in 1:ninterims) {
    points(x = rep(rdata$tinterims[i], ss), y = rdata$intnumbers[i,
      ], pch = pch, col = pcol, ...)
    lines(x = rep(rdata$tinterims[i], 2), y = c(rdata$intnumbers[i,
      ss], rdata$intnumbers[i, 1]), lty = lty[2], col = pcol, ...)
    if (is.na(labpos) != TRUE) {
      text(x = rdata$tinterims[i], y = rdata$intnumbers[i,
        1], if (any(is.na(intlab) == TRUE)) {
        substitute(t[j], list(j = i))
      } else {
        intlab[i]
      }, pos = labpos)
    }
  }
  if (mod$rdata$rmodel == "none") {
    legend(legpos, legend = c("Recruitment", as.character(tfu),
      "Interims"), lty = c(rep(lty[1], ss), lty[2]), pch = c(rep(NA,
      ss), pch), col = c(col, pcol), cex = 0.9, bty = "n")
  } else {
    legend(legpos, legend = c("Recruitment", as.character(tfu),
      "Interims"), lty = c(rep(lty[1], ss), lty[2]), pch = c(rep(NA,
      ss), pch), col = c(col, pcol), cex = 0.9, bty = "n")
  }
  ## Plot data
  names(summ_plotdata) <- c("t", alltfu)
  invisible(list(interims = rdata$interims, windowop = windowop, data = summ_plotdata))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

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