R/utility.R

In tsdf: Two-/Three-Stage Designs for Phase 1&2 Clinical Trials

#' Plot a decision table from a "dec.table" object
#' @description \code{plot} method for class \code{"dec.table"}.
#' @param x An object of class \code{"dec.table"}, typically returned by
#'   \code{\link{dec.table}}.
#' @param ... Unused arguments.
#' @details \code{plot.dec.table} plots the decision boundaries.
#' @import graphics
#' @export
#' @examples
#' truep <- c(0.3, 0.45, 0.5, 0.6)
#' out <- dec.table(0.6,0.4,0.2,0.3,c(3,3,3))
#' plot(out)
plot.dec.table <- function(x, ...) {
  n <- x$n
  nc <- cumsum(n)
  r <- x$E
  s <- x$D
  su <- x$DU
  col <- c("cadetblue1", "darkorchid2", "indianred2", "seagreen2")
  plot(nc, r, xaxt = "n", ylab = "Boundary", xlab = 'Sample Size', ylim = c(0, max(su)+2.3), type = "o", col=col[1], main = "Decision Plot", panel.first = grid())
  points(nc, s + 1, type = "o", col = col[2])
  points(nc, su + 1, type = "o", col = col[3])
  axis(1, at = nc, labels = nc)
  polygon(c(nc, rev(nc)), c(r + 0.05, rev(s + 1 - 0.05)), col = col[4], border = NA)
  polygon(c(nc, rev(nc)), c(rep(0,length(nc)), rev(r - 0.05)), col = col[1], border = NA)
  polygon(c(nc, rev(nc)), c(s + 1 + 0.05, rev(su + 1 - 0.05)), col = col[2], border = NA)
  polygon(c(nc, rev(nc)), c(su + 1 + 0.05, rev(rep(max(su)+2, length(nc)))), col = col[3], border = NA)
  legend("top", c("E (<=)", "D (>=)", "DU (>=)", "S"), horiz = TRUE, fill = col)
}


#' Print a decision table from a "dec.table" object
#' @description \code{print} method for class \code{"dec.table"}.
#' @param x An object of class \code{"dec.table"}, typically returned by
#'   \code{\link{dec.table}}.
#' @param ... Unused arguments.
#' @details \code{print.dec.table} prints the decision table together with a
#'   legend.
#' @export
#' @examples
#' print(dec.table(0.6,0.4,0.2,0.3,c(3,3,3)))
print.dec.table <- function(x, ...) {
  print(x$table)
  cat("\n")
  cat("	  Row : Number of DLTs", "\n")
  cat("Column : Number of subjects", "\n")
  cat("     E : Escalate to the next higher dose", "\n")
  cat("     S : Stay at the same dose", "\n")
  cat("     D : De-escalate to the previous lower dose", "\n")
  cat("    DU : De-escalate and never use this dose again", "\n")
  cat("\n")
  cat("Type 1 error : ", "\n")
  print(x$alpha.two)
  cat("Type 2 error : ", x$beta,"\n")
  cat("Right-side type 1 error : ", "\n")
  print(x$alpha.one)
}

#' Plot simulation results from a "dec.sim" object
#' @description Available plot types are: true toxicity at each dose level
#'   (\code{type = "s"}); a bar plot of the probability of selecting each dose
#'   as the MTD (\code{type = "prob"}); a bar plot of the average number of
#'   patients treated at each dose level (\code{type = "np"}); and a bar plot
#'   of the average number of DLTs at each dose level (\code{type = "dlt"}).
#'   Setting \code{type = "all"} produces all four plots.
#' @param x An object of class \code{"dec.sim"} or \code{"sl.sim"}, returned
#'   by \code{\link{dec.sim}} or \code{\link{sl.sim}}.
#' @param pt A vector of target toxicity values, one for each scenario.
#' @param s The scenario to plot. Defaults to 1.
#' @param type Plot type. See the description above.
#' @param label Logical; if \code{TRUE}, values are displayed on the plot.
#' @param col Graphical parameter \code{col}; see \code{\link{par}}.
#' @param text.col Color used for text labels.
#' @param cex Graphical parameter \code{cex}; see \code{\link{par}}.
#' @param ... Arguments passed to plotting functions.
#' @import graphics
#' @export
#' @examples 
#' # generate decision table
#' dt <- dec.table(0.6,0.4,0.2,0.3,c(3,3,3))
#' # Simulate trials from test data
#' test.file <- system.file("extdata", "testS.csv", package = "tsdf")
#' out <- sl.sim(dt$table, test.file)
#' plot(out, pt=rep(0.3,2), s=1, type="all")
#' plot(out, pt=rep(0.3,2), s=2, type="prob")
#' plot(out, pt=rep(0.3,2), s=1, type="np")
#' plot(out, pt=rep(0.3,2), s=2, type="dlt")
plot.dec.sim <- function(x, pt, s = 1, type = c("all", "s", "prob", "np", "dlt"), label = TRUE, col = "cornflowerblue", text.col = "darkblue", cex = 1, ...) {
  type <- match.arg(type)
  if (class(x)[1] == "sl.sim") {
    obj <- x[[s]]
    ns <- length(x)
  } else {
    obj <- x
    ns <- 1
  }
  truep <- obj$truep
  n_dose <- length(truep)
  names.arg <- 1:n_dose
  scen <- paste("( Scenario", s, ")")
  if (missing(pt)) {
    print("No target toxicity/mtd shown")
  }
  else if (length(pt) != ns) {
    stop("true toxicity is a vector with length = # of cenarios")
  } else {
    if(max(truep) < pt[s]) {
      print("MTD is higher than highest dose")
    } else if(min(truep) > pt[s]) {
      print("MTD is lower than lowest dose")
    } 
    else {
      mtd <- max(which(truep <= pt[s]))
      names.arg[mtd] <- paste(mtd, "(MTD)")
    }
  }
  if (type == "all") {
    par(mfrow = c(2, 2), cex = cex)
    plot(x, pt, s, "s", label, col, text.col, cex, ...)
    plot(x, pt, s, "prob", label, col, text.col, cex, ...)
    plot(x, pt, s, "dlt", label, col, text.col, cex, ...)
    plot(x, pt, s, "np", label, col, text.col, cex, ...)		
  }
  if (type == "prob") {
    ans <- obj$mtd.prob[1:n_dose]
    names(ans) <- NULL
    out <- barplot(rep(NA, n_dose), xlab = "Dose level", names.arg = names.arg, panel.first = box(), ylim = c(0, max(ans) * 1.1))
    grid()
    barplot(ans, add = TRUE, col = col)
    title(paste("Probability of selection", scen))
    if (label) {
      text(out, ans, labels = ans, pos = 3, cex = cex, col = text.col)
    }
  }
  if (type == "dlt") {
    ans <- obj$dlt
    out <- barplot(rep(NA, n_dose), xlab = "Dose level", names.arg = names.arg, panel.first = box(), ylim = c(0, max(ans) * 1.1))
    grid()
    barplot(ans, add = TRUE, col = col)
    title(paste("Average number of DLTs", scen))
    if (label) {
      text(out, ans, labels = ans, pos = 3, cex = cex, col = text.col)
    }
  }
  if (type == "np") {
    ans <- obj$n.patients
    out <- barplot(rep(NA, n_dose), xlab = "Dose level", names.arg = names.arg, panel.first = box(), ylim = c(0, max(ans) * 1.1))
    grid()
    barplot(ans, add = TRUE, col = col)
    title(paste("Average number of patients", scen))
    if (label) {
      text(out, ans, labels = ans, pos = 3, cex = cex, col = text.col)
    }
  }
  if (type == "s") {
    out <- plot(1:n_dose, truep, xlab = "Dose level", ylab = "Toxicity", col = col, type = "b", ylim = c(min(truep) * 0.9, max(truep) * 1.1), xaxt = "n", pch = 19, panel.first = c(box(), grid()))
    axis(1, at = 1:n_dose, labels = names.arg)
    title(paste("True toxicity", scen))
    if (!missing(pt)) {
      abline(pt[s], 0, lty = 5)
      text(n_dose/2, pt[s], labels = paste("Target =", pt[s]), pos = 3, cex = cex, col = text.col)
    }
    if (label) {
      text(1:n_dose, truep, labels = truep, pos = 3, cex = cex, col = text.col)
    }
  }
}


#' Summarize simulation results from a "dec.sim" object
#' @description \code{summary} method for class \code{"dec.sim"}.
#' @details \code{summary} formats key statistics from dose-finding
#'   simulations. Given the target toxicity, it reports the probability of
#'   selecting each dose level as the MTD, the probability of selecting doses
#'   above the true MTD, the probability of selecting the true MTD, and the
#'   probability that subjects are treated at or below the true MTD. The true
#'   MTD is defined as the highest dose level with toxicity probability less
#'   than or equal to the target toxicity. If the target is below the smallest
#'   toxicity probability, the lowest dose level is treated as the MTD. For
#'   example, if the target is 0.3 and the true toxicity probabilities for five
#'   doses are 0.1, 0.25, 0.35, 0.40, and 0.50, then the MTD is dose 2.
#' @param object An object of class \code{"dec.sim"}, returned by
#'   \code{\link{dec.sim}} or \code{\link{sl.sim}}.
#' @param pt Target toxicity for each scenario.
#' @param ... Unused arguments.
#' @examples 
#' test.file <- system.file("extdata", "testS.csv", package = "tsdf")
#' dt <- dec.table(0.6,0.4,0.2,0.3,c(3,3,3))
#' out <- sl.sim(dt$table, test.file)
#' pt <- c(0.3, 0.4)
#' summary(out, pt)
#' @import stats
#' @export
summary.dec.sim <- function(object, pt, ...) {
  if(missing(pt)) {
    warning("Missing true toxicity")
  }
  if(class(object)[1] == "sl.sim" & length(object) != length(pt)) {
    warning("pt length not equal to number of scenarios; only returns the first scenario stats")
  }
  ns <- length(pt)
  avg.np <- rep(0, ns)
  avg.dose <- rep(0, ns)
  avg.prob <- rep(0, ns)
  out <- vector("list", ns)
  for(i in 1:ns) {
    if(class(object)[1] == "sl.sim"){
      ans <- object[[i]]
    } else {
      ans <- object
    }
    truep <- ans$truep
    n_dose <- length(truep)
    res <- matrix(0, n_dose, 6)
    colnames(res) <- c("Level", "Truth", "MTD", "Over","DLT", "NP")
    res[, 1] <- 1:n_dose 
    res[, 2] <- truep
    res[, 3] <- ans$mtd.prob[1:n_dose]
    res[, 4] <- ans$over.prob[1:n_dose]
    res[, 5] <- ans$dlt
    res[, 6] <- ans$n.patients
    # calculate stats
    avg.np[i] <- sum(res[, 6])
    if(max(truep) < pt[i]) {
      mtd.dose <- "higher than highest dose"
      avg.dose[i] <- ans$mtd.prob[n_dose+2]
      avg.prob[i] <- 1
    } else if(min(truep) > pt[i]) {
      mtd.dose <- "lower than lowest dose"
      avg.dose[i] <- ans$mtd.prob[n_dose+1]
      avg.prob[i] <- 0
    } else {
      mtd.dose <- max(which(truep <= pt[i]))
      avg.dose[i] <- res[, 3][mtd.dose]
      avg.prob[i] <-sum(res[, 6][truep <= pt[i]])/sum(res[, 6])
    }
    out[[i]] <- list(dose.stats = res, prob.select = avg.dose[i], at.below.mtd = avg.prob[i], mtd = mtd.dose, nsim = ans$nsim, pt = pt[i], avg.np = avg.np[i], start.level = ans$start.level)
  }
  class(out) <- "summary.dec.sim"
  out
}

#' @export
print.summary.dec.sim <- function(x, ...) {
  cat("What does each column represent ?", "\n\n")
  cat("Level : Dose level", "\n")
  cat("Truth : True toxicity probability", "\n")
  cat("MTD   : The probability of selecting current dose level as the MTD", "\n")
  cat("Over  : The probability of selecting current dose level as over the MTD", "\n")
  cat("DLT   : The average number of subjects experienced DLT at current dose level", "\n")
  cat("NP    : The average number of subjects treated at current dose level", "\n\n")
  
  for(i in 1:length(x)){
    # table legends
    cat("Scenario ", i, "\n\n")
    cat("Simulation results are based on", x[[i]]$nsim, "simulated trials","\n")
    cat("Starting dose level is", x[[i]]$start.level, "; MTD is dose", x[[i]]$mtd, "\n")
    cat("Target toxicity probability =", x[[i]]$pt, "\n")
    cat("Average number of subjects =", x[[i]]$avg.np, "\n")
    cat("Probability of selecting the true MTD =", x[[i]]$prob.select, "\n")
    cat("Probability of subjects treated at or below the true MTD =", x[[i]]$at.below.mtd, "\n\n")
    print(as.data.frame(x[[i]]$dose.stats))
    cat("\n")
  }
}


#' Print Zhong's design from an "opt.design" object
#' @description \code{print} method for class \code{"opt.design"}.
#' @param x An object of class \code{"opt.design"}, typically returned by
#'   \code{\link{opt.design}}.
#' @param ... Unused arguments.
#' @export
#' @examples
#' alpha1 <- 0.20
#' alpha2 <- 0.1
#' beta <- 0.20
#' pc <- 0.5
#' pt <- pc + 0.2
#' out <- opt.design(alpha1, alpha2, beta, pc, pt, stage = 2, sf.param = 1)
#' print(out)
print.opt.design <- function(x, ...) {
  cat("\n Zhong's", x$stage,"stage Phase II design \n\n")
  cat("Minimal response rate: ", unique(x$pe), "\n")
  cat("Postulate response rate: ", x$pt, "\n")
  cat("Left-side type 1 error: ",x$alpha1, "\n")
  cat("Right-side type 1 error: ",x$alpha2, "\n")
  cat("Type 2 error: ",x$beta, "\n\n")
  cat("Notation:", "\n")
  cat("Left-side rejection region at stage i is response <= ri", "\n")
  cat("Right-side rejection region at stage", x$stage, "is response > s", "\n")
  cat("Sample size used at stage i is ni", "\n\n")
  cat("Optimal design : ", "\n")
  print(c(x$bdry, x$n))
  cat("True errors : ", "\n")
  print(x$error)
  cat("\n")
}