R/visit_analysis.R

In visit: Vaccine Phase I Design with Simultaneous Evaluation of Immunogenicity and Toxicity

#' Obtain decision map information
#'
#' Summarize the posterior distribution of \eqn{\theta^{(l)}_{00},
#'     \theta^{(l)}_{01}, \theta^{(l)}_{10}, \theta^{(l)}_{11}} and get
#'     information for making dose escalation decisions
#'
#' @rdname vtDecMap
#'
#' @inheritParams parameters
#'
#' @param thetas Posterior samples of \eqn{\theta}, a class \code{VTPOST} matrix
#'     generated by \code{\link{vtPost}}
#'
#' @details
#'
#' This function summarizes the posterior distribution of the
#'     \eqn{\theta^{(l)}_{00}, \theta^{(l)}_{01}, \theta^{(l)}_{10},
#'     \theta^{(l)}_{11}} and sequentially get the conditional probabilities of
#'     each decision map region. See \code{\link{visit}} for details of the
#'     decision map regions.
#'
#' @return
#'
#' A class \code{VTDEC} list. See the return value from \code{\link{vtInterim}}
#' for details.
#'
#' @examples
#' etas     <- c(0.1, 0.3)
#' dec.cut  <- c(0.6,0.6,0.6)
#' obs.y    <- rbind(c(5, 2, 0, 0))
#' rst.post <- vtPost(obs.y,  prob.mdl = "NONPARA", nsmp = 2000)
#' dec.map  <- vtDecMap(rst.post, etas = etas, dec.cut = dec.cut)
#'
#' @export
#'
vtDecMap <- function(thetas, etas, prev.res=0, dec.cut=0.6) {

    stopifnot(get.const()$CLSPOST %in% class(thetas));

    ## allow 3 cuts
    dec.cut <- rep(dec.cut, 3)[1:3];

    ## posterior region probabilities
    cur.smp <- thetas;
    cur.dlt <- cur.smp[,3] + cur.smp[,4];
    cur.res <- cur.smp[,2] + cur.smp[,4];

    rst    <- rep(0,4);
    ##d.reg.1: over toxic
    rst[1] <- mean(cur.dlt >= etas[2]);
    ##d.reg.2: less effective
    rst[2] <- mean(cur.dlt < etas[2] & cur.res < prev.res);
    ##d.reg.3: safe and effective
    rst[3] <- mean(cur.dlt < etas[1] & cur.res >= prev.res);
    ##d.reg.4: not over toxic, not safe, effective
    rst[4] <- 1 - sum(rst[1:3]);

    ##conditional approach
    cond.prob    <- rep(0,2);
    cond.prob[1] <- rst[2]/(1 - rst[1] + 1e-10); ## avoid Nan
    cond.prob[2] <- rst[3]/(1 - sum(rst[1:2]) + 1e-10);

    if (rst[1] > dec.cut[1]) {
        region         <- 1;
        cond.prob[1:2] <- 0;
    } else if (cond.prob[1] > dec.cut[2]) {
        region <- 2;
        cond.prob[2] <- 0;
    } else if (cond.prob[2] > dec.cut[3]) {
        region <- 3;
    } else  {
        region <- 4;
    }

    rst.dec <- list(prob      = rst,
                    region    = region,
                    ptox      = mean(cur.dlt),
                    pres      = mean(cur.res),
                    cond.prob = cond.prob,
                    prev.res  = mean(prev.res),
                    etas      = etas,
                    dec.cut   = dec.cut);

    class(rst.dec) <- get.const()$CLSDEC;

    rst.dec
}

#' Conduct interim analysis
#'
#' Conduct an interim analysis for determining dose escalation actions
#'
#' @rdname vtInterim
#'
#' @inheritParams parameters
#'
#' @param cur.obs.y Observed data from the current level, which is a vector of
#'     length 4. The numbers correspond to \code{obs.y} in \code{\link{vtPost}}.
#' @param prev.obs.y Observed data from previous levels, which has the same
#'     structure as \code{obs.y} in \code{\link{vtPost}}.
#' @param ... Additional arguments for \code{\link{vtPost}}
#'
#' @details
#'
#' Using data from previous levels and the current level to conduct Bayesian
#' analysis, get the decision map information and make decision about dose
#' escalation actions. The actions include stop the trial, escalate to the next
#' higher dose level, or enroll more patients in the current level. See \code{\link{visit}}
#' for details.
#'
#' @return
#' A class \code{VTDEC} list containing
#' \itemize{
#' \item{prob: }{Probabilities of each decision map region}
#' \item{region: }{The region selected based on the sequential procedure described in \code{\link{visit}} }
#' \item{ptox: }{Mean risk of DLT, \eqn{E(p^{(l)})}}
#' \item{pres: }{Mean immune response rate, \eqn{E(q^{(l)})}}
#' \item{con.prob: }{Conditional probabilities of each decision map region}
#' \item{prev.res: }{Function parameter}
#' \item{etas: }{Function parameter}
#' \item{dec.cut: }{Function parameter}
#' }
#'
#' @examples
#'
#' etas       <- c(0.1, 0.3)
#' dec.cut    <- c(0.6,0.6,0.6)
#' cur.obs.y  <- c(3, 2, 1, 1)
#' prev.obs.y <- c(5, 2, 0, 0)
#' rst.inter  <- vtInterim(cur.obs.y,  prev.obs.y = prev.obs.y,
#'                         prob.mdl = "NONPARA", etas = etas,
#'                         dec.cut = dec.cut,
#'                         nsmp = 2000);
#'
#' @export
#'
vtInterim <- function(cur.obs.y, prev.obs.y = NULL, prev.res = NULL,
                      etas        = c(0.1,0.3),
                      dec.cut     = 0.65,
                      priors      = NULL,
                      prob.mdl    = c("NONPARA", "NONPARA+", "PARA", "PARA+"),
                      seed        = NULL,
                      ...) {

    if (!is.null(seed)) {
        old_seed <- set.seed(seed);
    }

    prob.mdl <- match.arg(prob.mdl);

    ##bayesian
    post.smp  <- vtPost(rbind(cur.obs.y), prob.mdl, priors, ...);

    if (!is.null(prev.obs.y)) {
        post.prev <- vtPost(rbind(prev.obs.y), prob.mdl, priors, ...);
        prev.res  <- apply(post.prev[,c(2,4)],1,sum);
    } else if (is.null(prev.res)){
        stop("Please provide either prev.obs.y or prev.res");
    }

    ## decision
    rst.dec  <- vtDecMap(post.smp, etas, prev.res=prev.res, dec.cut=dec.cut);

    ## reset random see
    if (!is.null(seed)) {
        set.seed(old_seed);
    }

    rst.dec
}

#' Plot decision map
#'
#' Plot a decision map based on a class \code{VTDEC} object that contains the
#' current posterior analysis results
#'
#' @rdname plot.VTDEC
#'
#' @param x A class \code{VTDEC} list generated by \code{\link{vtDecMap}}
#' @param margin Margin between regions in the decision map
#' @param nms Labels of the regions on a decision map. Defaults are:
#' \itemize{
#' \item{\code{TT}:}{Too Toxic}
#' \item{\code{NME}:}{No More Effective}
#' \item{\code{SE}:}{Safe and Effective}
#' \item{\code{UN}:}{Uncertain}
#' }
#' @param col.reg Background color of the selected region
#' @param col.prob Text color of the selected region.
#' @param cex.prob Text size of the probabilities
#' @param cex.nms Text size of the region labels
#' @param ... Optional arguments for \code{plot}.
#'
#' @examples
#'
#' etas       <- c(0.1, 0.3)
#' dec.cut    <- c(0.6,0.6,0.6)
#' cur.obs.y  <- c(3, 2, 1, 1)
#' prev.obs.y <- c(5, 2, 0, 0)
#' rst.inter  <- vtInterim(cur.obs.y,  prev.obs.y = prev.obs.y,
#'                         prob.mdl = "NONPARA", etas = etas, dec.cut = dec.cut,
#'                         nsmp = 2000);
#' plot(rst.inter)
#'
#' @method plot VTDEC
#'
#' @export
#'
plot.VTDEC <- function(x,
                       margin = 0.003, nms = c("TT", "NME", "SE", "UN"),
                       col.reg = "pink", col.prob = "blue", cex.prob = 0.9, cex.nms = 1,
                       ...) {

    f.reg <- function(x1, x2, y1, y2, labels, cols, tt = margin) {

        x1 <- x1 + tt;
        x2 <- x2 - tt;
        y1 <- y1 + tt;
        y2 <- y2 - tt;

        rect(x1, y1, x2, y2, col = cols[1], lwd = 1);
        text(x1+(x2-x1)/2, y1+(y2-y1)/2, labels = labels[1], col=cols[2], cex=cex.prob);
        text(x1+(x2-x1)/2, y1+(y2-y1)/2, labels = labels[2], col=cols[2], cex=cex.prob, pos = 1);
        text(x1+(x2-x1)/2, y2, labels = labels[3], cex=cex.nms, pos = 1);
    }

    ##
    etas      <- x$etas;
    prev.res  <- round(x$prev.res,2);
    probs     <- x$prob * 100;
    cond.prob <- x$cond.prob * 100;

    col.regs  <- rep("white", 4);
    col.ps    <- rep("black", 4);
    col.regs[x$region] <- col.reg;
    col.ps[x$region] <- col.prob;

    ##plot
    par(xaxs="i", yaxs="i");
    plot(NULL, xlim=c(0,1), ylim=c(0,1),
         xlab="DLT Risk", ylab="Immune Response Rate",
         axes=FALSE, ...);
    ##box(lwd = 3);
    ##axis(1, at = round(etas,2));

    ##region I
    lbls <- c(sprintf("%.2f", probs[1]), "", nms[1]);
    f.reg(etas[2], 1, 0, 1, labels = lbls, cols=c(col.regs[1], col.ps[1]));

    ##region II
    if (prev.res > 0) {
        axis(2, at = prev.res);
        l2 <- "";
        if (0 < cond.prob[1]) {
            l2 <- sprintf("(%.2f)", cond.prob[1]);
        }
        lbls <- c(sprintf("%.2f", probs[2]), l2, nms[2]);
        f.reg(0, etas[2], 0, prev.res,
              labels = lbls, cols=c(col.regs[2], col.ps[2]));
    }

    ##region III
    l2 <- "";
    if (0 < cond.prob[2]) {
        l2 <- sprintf("(%.2f)", cond.prob[2]);
    }
    lbls <- c(sprintf("%.2f", probs[3]), l2, nms[3]);
    f.reg(0, etas[1], prev.res, 1,
          labels = lbls, cols=c(col.regs[3], col.ps[3]));

    ##region IV
    lbls <- c(sprintf("%.2f", probs[4]), "", nms[4]);
    f.reg(etas[1], etas[2], prev.res, 1,
          labels = lbls, cols=c(col.regs[4], col.ps[4]));
}


#' Plot the track plot of dose escalation
#'
#' Generate a plot representing the observed data and dose escalation decisions.
#'
#' @rdname vtTrack
#'
#' @param obs.all All observations collected in a matrix with 5 columns. Column
#'     1 is the index of interim analysis starting from 1. Columns 2-5
#'     correspond to columns 1-4 in \code{obs.y} for \code{\link{vtPost}}.
#' @param cex.txt Text size of numbers in the plot
#' @param decision Dose escalation decision. The options are
#' \itemize{
#' \item{\code{1}: }{Escalate}
#' \item{\code{2}: }{Continue at the same level}
#' \item{\code{3}: }{Stop the trial}
#' }
#'
#' @param max.level Maximum number of dose levels shown in the plot
#' @param letters Labels for dose escalation actions 1-3. Default values are
#'     "E", "C", "S"
#' @param colors Possible colors in the last action box
#' @param height Height of each individual box
#' @param end.width Width of the last action box
#' @param end.height Height of the last action box
#' @param cex.roman Text size of the roman numerals
#' @param cex.end Text size of the letter in the last action box
#' @param ... Optional arguments for \code{plot}.
#'
#'
#' @examples
#'
#' obs.all  <- rbind(c(1, 5, 2, 0, 0),
#'                   c(2, 3, 4, 0, 0),
#'                   c(3, 1, 6, 0, 0));
#' vtTrack(obs.all, end.width = 0.8, max.level = 3, decision = 3);
#'
#' @export
#'
vtTrack <- function(obs.all,
                    cex.txt = 0.9, decision = 1, max.level = NULL,
                    letters = c("E", "C", "S"),
                    colors = c("green", "yellow", "red"), height = 0.5,
                    end.width = 2, end.height = height, cex.roman = 0.9, cex.end = 0.9,
                    ...) {

    f.rec <- function(x, y, width, height, ys, cex.txt, margin = 0.1) {
        width <- width/4;
        cys   <- rbind(c(0,0), c(0,1), c(1,0), c(1,1));
        for (i in 1:4) {
            cur.x   <- x + width *(i-1);
            xleft   <- cur.x;
            xright  <- xleft + 0.5 * width;
            ytop    <- y;
            ybottom <- ytop - 0.5 * height;

            if (1 == cys[i,1]) {
                r.txt <- "T";
                r.col <- "gray80";
            } else {
                r.txt <- "NoT";
                r.col <- "white";
            }
            rect(xleft, ybottom, xright, ytop, col=r.col, lwd=0.1);
            text(xleft + 0.25 * width, ytop - 0.25 * height, labels = r.txt, cex = 0.5 * cex.txt);

            xleft   <- cur.x + 0.5 * width;
            xright  <- xleft + 0.5 * width;
            ytop    <- y;
            ybottom <- ytop - 0.5 * height;

            if (1 == cys[i,2]) {
                r.txt <- "R";
                r.col <- "pink";
            } else {
                r.txt <- "NoR";
                r.col <- "white";
            }
            rect(xleft, ybottom, xright, ytop, col=r.col, lwd=0.1);
            text(xleft + 0.25 * width, ytop - 0.25 * height, labels = r.txt, cex = 0.5 * cex.txt);

            xleft   <- cur.x;
            xright  <- xleft + width;
            ytop    <- y - 0.5 * height;
            ybottom <- ytop - 0.5 * height;

            rect(xleft, ybottom, xright, ytop, lwd=0.1);
            text(xleft + 0.5 * width, ytop - 0.25 * height, labels = ys[i], cex = cex.txt);

            rect(cur.x, y - height, cur.x + width, y, lwd=1);
        }
    }

    f.arrow <- function(x, y, width, height, direction, labels = NULL) {

        end.xy <- switch(direction,
                         "flat" = c(width, 0, 3),
                         "up"   = c(0, height, 2),
                         "down" = c(0, -height, 4));

        arrows(x, y, x + end.xy[1], y + end.xy[2], length = 0.1);
        text(x + end.xy[1]/2, y + end.xy[2]/2,
             labels = as.roman(labels), pos = end.xy[3], col = "brown", cex = cex.roman);
    }

    f.end <- function(x, y, width, height, letter, color) {
        rect(x - 0.5*width, y - 0.5*height, x + 0.5*width, y + 0.5*height, lwd = 1);
        rect(x - 0.45*width, y - 0.45*height, x + 0.45*width, y + 0.45*height,
             lwd = 1, col = color);
        text(x, y, labels = letter, cex = cex.end);
    }

    if (is.null(max.level))
        max.level <- max(obs.all[,1]);

    tbl.level <- table(obs.all[,1]);
    x.tot     <- 0;
    for (i in 1:length(tbl.level)) {
        x.tot <- x.tot + 4*tbl.level[i] + tbl.level[i] - 1;
    }
    x.tot <- x.tot - length(tbl.level) + 1;

    x.max <- (x.tot + 3) * 1;
    y.max <- max.level * 1;

    plot(NULL, NULL, xlim=c(0, x.max), ylim=c(height/2, y.max),
         axes = FALSE, xlab = "", ylab = "", ...);
    box();

    last.level <-  0;
    cur.x      <-  1;
    cur.y      <-  1;
    inx        <-  0;
    for (i in 1:nrow(obs.all)) {
        cur.level <- obs.all[i,1];
        if (cur.level != last.level) {
            last.level <- cur.level;
            if (cur.level > 1) {
                inx <- inx + 1;
                f.arrow(cur.x + 3.5, cur.y, 0, 1-height, "up", labels = inx);
                cur.x <- cur.x + 3;
                cur.y <- cur.y + 1;
            }
        } else {
            inx  <- inx + 1;
            f.arrow(cur.x + 4, cur.y - 0.5*height, 1, 0, "flat", labels = inx);
            cur.x <- cur.x + 5;
        }
        f.rec(cur.x, cur.y, 4, height = height, ys = obs.all[i, 2:5], cex.txt = cex.txt);
    }

    inx  <- inx + 1;
    f.arrow(cur.x + 4, cur.y - 0.5*height, 1, 0, "flat", labels = inx);
    f.end(cur.x + 5 + end.width/2, cur.y - 0.5*height, width = end.width, height = end.height,
          letter = letters[decision], color = colors[decision]);
}