R/plot_kb.R

In Keyboard: Bayesian Designs for Early Phase Clinical Trials

#' Plot Objects Generated by the Keyboard Design Package
#'
#' This function plots the objects returned by the other functions in the Keyboard package. \cr\cr
#'
#' (1) For get.oc.kb() and get.oc.comb.kb(), the operating characteristics of
#' the Keyboard design, which include: \cr
#' (i) the overall dose selection percentage of the simulations generated, \cr
#' (ii) the number of patients treated at each dose, and \cr
#' (iii) the toxicities presented at each dose. \cr\cr
#' (2) For get.oc.obd.kb() and get.oc.obd2.kb(), the operating characteristics of the Keyboard design for phase I/II, which includes: \cr
#' (i) the overall dose selection percentage of the simulations generated, \cr
#' (ii) the number of patients treated at each dose,  \cr
#' (iii) the toxicities presented at each dose, and \cr
#' (iv) the efficacies presented at each dose. \cr\cr
#' (3) For select.mtd.kb() and select.mtd.comb.kb(): an infographic showing the
#' estimates of the toxicity probability for each dose, with corresponding
#' 95\% credible intervals.
#'
#'
#'
#' @param x The object (returned by another function) to be plotted.
#' @param ... Ignored arguments.
#' @param name The name of the object to be plotted.\cr
#'             (Users should not input this parameter).
#'
#' @return This function returns a figure.
#' @import Rcpp methods graphics stats ggplot2
#' @author Xiaomeng Yuan, Chen Li, Hongying Sun, Li Tang and Haitao Pan
#' @examples
#' \donttest{
#' ### Single-agent trial ###
#'
#' ## Get the operating characteristics for a Keyboard single-agent trial
#' oc <- get.oc.kb(target=0.3, p.true=c(0.05, 0.15, 0.3, 0.45, 0.6),
#'                 ncohort=10, cohortsize=3, ntrial=1000)
#' summary_kb(oc)
#'
#' plot_kb(oc)
#' plot_kb(oc$selpercent)
#' plot_kb(oc$npatients)
#' plot_kb(oc$ntox)
#'
#' ## Select the MTD based on a trial's data
#' n <- c(3, 3, 15, 9, 0)
#' y <- c(0, 0, 4, 4, 0)
#' selmtd <- select.mtd.kb(target=0.3, npts=n, ntox=y)
#' summary_kb(selmtd)
#'
#' plot_kb(selmtd)
#'
#'
#' ### Drug-combination trial ###
#'
#' ## Get the operating characteristics for a Keyboard drug-combination trial
#' p.true <- matrix(c(0.01, 0.03, 0.10, 0.20, 0.30,
#'                    0.03, 0.05, 0.15, 0.30, 0.60,
#'                    0.08, 0.10, 0.30, 0.60, 0.75), byrow=TRUE, ncol=5)
#' oc.comb <- get.oc.comb.kb(target=0.3, p.true=p.true, ncohort=20,
#'  cohortsize=3, n.earlystop=12, startdose=c(1, 1), ntrial=100)
#' summary_kb(oc.comb)
#'
#' plot_kb(oc.comb) # use previous plot or next plot button to switch plots
#' plot_kb(oc.comb$selpercent)
#' plot_kb(oc.comb$nptsdose)
#' plot_kb(oc.comb$ntoxdose)
#'
#' ## Select the MTD based on a trial's data
#' n <- matrix(c(3, 5, 0, 0, 0,
#'               7, 6, 15, 0, 0,
#'               0, 0, 4, 0, 0), ncol=5, byrow=TRUE)
#' y <- matrix(c(0, 1, 0, 0, 0,
#'               1, 1, 4, 0, 0,
#'               0, 0, 2, 0, 0), ncol=5, byrow=TRUE)
#' sel.comb <- select.mtd.comb.kb(target=0.3, npts=n, ntox=y)
#' summary_kb(sel.comb)
#'
#' plot_kb(sel.comb)
#'
#' ### oc.obd.kb
#' toxicity.low <- 0.15
#' toxicity.moderate <- 0.33
#' toxicity.high <- 0.40
#' efficacy.low <- 0.20
#' efficacy.moderate <- 0.40
#' efficacy.high <- 0.60
#' target.toxicity <- 0.30
#' target.efficacy <- 0.40
#' p.true <-c(0.08,0.30,0.60,0.80)
#' q.true <- c(0.25,0.40,0.25,0.50)
#' oc.obd.kb <- get.oc.obd.kb(toxicity.low = toxicity.low,
#'              toxicity.moderate= toxicity.moderate,
#'              toxicity.high = toxicity.high,
#'              efficacy.low = efficacy.low,
#'              efficacy.moderate = efficacy.moderate,
#'              efficacy.high = efficacy.high,
#'              target.toxicity=target.toxicity,
#'              target.efficacy= target.efficacy,
#'              p.true= p.true, q.true= q.true)
#' summary_kb(oc.obd.kb)
#' plot_kb(oc.obd.kb)
#' plot_kb(oc.obd.kb$selpercent1)
#' plot_kb(oc.obd.kb$selpercent2)
#' plot_kb(oc.obd.kb$selpercent3)
#' plot_kb(oc.obd.kb$npatients)
#' plot_kb(oc.obd.kb$ntox)
#' plot_kb(oc.obd.kb$neff)
#' }
#' @seealso \code{\link{summary_kb}}
#' @export

plot_kb <- function(x, ..., name = deparse(substitute(x))) {

    suppressWarnings(par(new=TRUE))

    # UseMethod("plot_kb")
    new.obj = unlist(strsplit(name, split = "\\$"))
    strpattern = "none"
    assign("objectPlot", get(new.obj[1])) #as.name() will convert a string into a reference
    if (length(new.obj) >= 2) {
        strpattern = new.obj[2]
    }
    ## If invalid option, e.g. plot_kb(oc.comb$flowchart)
    if (!is.element(strpattern, c("none", names(objectPlot)))) {
        stop("Please double check and specify the variable to be plotted...\n")
        ## edit: changed from warning()
    }

######################################################################################
    ## plotting for get.oc.kb()
    if (!is.null(objectPlot$percentstop)) {
        # dev.new(); op <- par(no.readonly = TRUE); dev.off(); par(op)
        opar <- par(no.readonly = TRUE)
        on.exit(par(opar))
        par(mfrow = c(1,1), mar = c(5,5,4,2))

        if (length(new.obj) == 2) { ## WHEN TRUE??? Shouldn't it be length(x)?
            if (strpattern == "selpercent") {
                bplot = barplot(objectPlot$selpercent,
                                ylab = "selection percentage (%)", ylim = c(0,100),
                                cex.names = 1, xaxt = "n", cex.lab = 1.3)
                axis(1, at = bplot, labels = seq(1, length(objectPlot$selpercent)))
                mtext("dose level", 1, line = 2, cex = 1.3)
            }
            else if (strpattern == "npatients") {
                bplot = barplot(objectPlot$npatients,
                                ylab = "number of patients", ylim = c(0, max(objectPlot$npatients)),
                                cex.names = 1, beside = FALSE, xaxt = "n", cex.lab = 1.3)
                axis(1, at = bplot, labels = seq(1, length(objectPlot$npatients)))
                mtext("dose level", 1, line = 2, cex = 1.3)
            }
            else if (strpattern == "ntox") {
                bplot = barplot(objectPlot$ntox,
                                ylab = "number of toxicities", ylim = c(0, max(objectPlot$ntox)),
                                cex.names = 1, beside = FALSE, xaxt = "n", cex.lab = 1.3)
                axis(1, at = bplot, labels = seq(1, length(objectPlot$ntox)))
                mtext("dose level", 1, line = 2, cex = 1.3)
            }
        }
        else {
            opar <- par(no.readonly = TRUE)
            on.exit(par(opar))
            par(mfrow = c(2,1), mar = c(5,5,4,2))

            bplot = barplot(objectPlot$selpercent,
                            ylab = "selection percentage (%)", ylim = c(0,100),
                            cex.names = 1, xaxt = "n", cex.lab = 1.3)
            mtext("selection percentage", 3, line = 2, cex = 1.3)
            axis(1, at = bplot, labels = seq(1, length(objectPlot$selpercent)))
            mtext("dose level", 1, line = 2, cex = 1.3)

            bplot = barplot(objectPlot$npatients,
                            ylab = "number of patients", ylim = c(0, max(objectPlot$npatients)),
                            cex.names = 1, beside = FALSE, xaxt = "n", cex.lab = 1.3)
            axis(1, at = bplot, labels = seq(1, length(objectPlot$npatients)))
            mtext("patient allocation", 3, line = 2, cex = 1.3)
            mtext("dose level", 1, line = 2, cex = 1.3)

            ## why not "number of toxicities" too?
        }
    }

    ######################################################################################
    ## Plotting for get.oc.comb.kb()
    else if(!is.null(objectPlot$pcs)) { # | !is.null(objectPlot$pcs.contour)) {

        # if (!requireNamespace("epade", quietly = TRUE)) {
        #             stop("Package \"epade\" needed for this function to work.",
        #      "Please install it.", call. = FALSE)
        #
        # }

        # dev.new(); op <- par(no.readonly = TRUE); dev.off(); par(op)
        #par(mar=c(0,0,0,0),oma=c(0,0,0,0))
         opar <- par(no.readonly = TRUE)
         on.exit(par(opar))
        # par(mfrow = c(1,1), mar = c(5,5,2,2))

        if (strpattern == "selpercent") {
            # requireNamespace("epade", quietly = TRUE)
            # epade::bar3d.ade(objectPlot$selpercent,
            #                  wall = 6, zw = 1, xw = 0.2,
            #                  yticks = seq(0,100,20),
            #                  zticks = seq(1:nrow(objectPlot$selpercent)),
            #                  xticks = seq(1:ncol(objectPlot$selpercent)),
            #                  xlab = "Drug B",
            #                  ylab = "selection percentage (%)",
            #                  zlab = "Drug A",
            #                  col = 'lavender', alpha = 0.5)

          yv <- as.vector(objectPlot[["selpercent"]])
          # drug A
          Drug_A <- rep(c(1: nrow(objectPlot[["selpercent"]])),ncol(objectPlot$selpercent))
          # drug B
          Drug_B <- rep(c(1: ncol(objectPlot[["selpercent"]])),each=nrow(objectPlot$selpercent))

          df <- as.data.frame(cbind(yv, Drug_A, Drug_B))

          print(ggplot(df, aes(Drug_A, yv)) +
                  ylab('selection percentage (%)')  +
                  xlab("Drug A")+
                  labs(title="selection percentage")+
                  geom_bar(stat="identity") +
                  facet_wrap(~ Drug_B,ncol=ncol(objectPlot$selpercent), labeller = label_both)+
                  theme(strip.text.x = element_text(size = 10, angle = 0)))


        }
        else if (strpattern == "nptsdose") { ## changed from BOIN's "npatients" (to match var in get.oc.comb.kb)
            # requireNamespace("epade", quietly = TRUE)
            # epade::bar3d.ade(objectPlot$nptsdose,
            #                  wall = 6, zw = 1, xw = 0.2,
            #                  yticks = round(seq(1,max(objectPlot$nptsdose),
            #                                 max(objectPlot$nptsdose)/7),0),
            #                  zticks = seq(1:nrow(objectPlot$nptsdose)),
            #                  xticks = seq(1:ncol(objectPlot$nptsdose)),
            #                  xlab = "Drug B",
            #                  ylab = "number of patients",
            #                  zlab = "Drug A",
            #                  col  = 'lavender', alpha = 0.5)

          yv <- as.vector(objectPlot[["nptsdose"]])
          # drug A
          Drug_A <- rep(c(1: nrow(objectPlot[["nptsdose"]])),ncol(objectPlot$nptsdose))
          # drug B
          Drug_B <- rep(c(1: ncol(objectPlot[["nptsdose"]])),each=nrow(objectPlot$nptsdose))

          df <- as.data.frame(cbind(yv, Drug_A, Drug_B))

          print(ggplot(df, aes(Drug_A, yv)) +
                  ylab('number of patients')  +
                  xlab("Drug A")+
                  labs(title="patient allocation")+
                  geom_bar(stat="identity") +
                  facet_wrap(~ Drug_B,ncol=ncol(objectPlot$nptsdose), labeller = label_both)+
                  theme(strip.text.x = element_text(size = 10, angle = 0)))

        }
        else if (strpattern == "ntoxdose"){ ## changed from BOIN's "ntox" (to match var in get.oc.comb.kb)
            # requireNamespace("epade", quietly = TRUE)
            # epade::bar3d.ade(objectPlot$ntoxdose,
            #                  wall = 6, zw = 1, xw = 0.2,
            #                  yticks = round(seq(1, max(objectPlot$ntoxdose),
            #                                 max(objectPlot$ntoxdose)/7),0),
            #                  zticks = seq(1:nrow(objectPlot$ntoxdose)),
            #                  xticks = seq(1:ncol(objectPlot$ntoxdose)),
            #                  xlab = "Drug B",
            #                  ylab = "number of toxicities",
            #                  zlab = "Drug A",
            #                  col = 'lavender', alpha = 0.5)

          yv <- as.vector(objectPlot[["ntoxdose"]])
          # drug A
          Drug_A <- rep(c(1: nrow(objectPlot[["ntoxdose"]])),ncol(objectPlot$ntoxdose))
          # drug B
          Drug_B <- rep(c(1: ncol(objectPlot[["ntoxdose"]])),each=nrow(objectPlot$ntoxdose))

          df <- as.data.frame(cbind(yv, Drug_A, Drug_B))

          print(ggplot(df, aes(Drug_A, yv)) +
                  ylab('number of toxicities')  +
                  xlab("Drug A")+
                  labs(title="number of toxicities")+
                  geom_bar(stat="identity") +
                  facet_wrap(~ Drug_B,ncol=ncol(objectPlot$ntoxdose), labeller = label_both)+
                  theme(strip.text.x = element_text(size = 10, angle = 0)))

        }

        else {
            if (strpattern == "none") {
                # dev.new(); op <- par(no.readonly = TRUE); dev.off(); par(op)
                # opar <- par(no.readonly = TRUE)
                # on.exit(par(opar))
                #
                # par(mfrow = c(3,1), mar = c(4,3,2,2))

                # epade::bar3d.ade(objectPlot$selpercent,
                #                  wall = 6, zw = 1, xw = 0.2,
                #                  main = "selection percentage",
                #                  yticks = seq(0,100,20),
                #                  zticks = seq(1:nrow(objectPlot$selpercent)),
                #                  xticks = seq(1:ncol(objectPlot$selpercent)),
                #                  xlab = "Drug B",
                #                  ylab = "selection percentage (%)",
                #                  zlab = "Drug A",
                #                  col = 'lavender', alpha = 0.5)

                yv <- as.vector(objectPlot[["selpercent"]])
                # drug A
                Drug_A <- rep(c(1: nrow(objectPlot[["selpercent"]])),ncol(objectPlot$selpercent))
                # drug B
                Drug_B <- rep(c(1: ncol(objectPlot[["selpercent"]])),each=nrow(objectPlot$selpercent))

                df <- as.data.frame(cbind(yv, Drug_A, Drug_B))

                print(ggplot(df, aes(Drug_A, yv)) +
                        ylab('selection percentage (%)')  +
                        xlab("Drug A")+
                        labs(title="selection percentage")+
                        geom_bar(stat="identity") +
                        facet_wrap(~ Drug_B,ncol=ncol(objectPlot$selpercent), labeller = label_both)+
                        theme(strip.text.x = element_text(size = 10, angle = 0)))


                # epade::bar3d.ade(objectPlot$nptsdose,
                #                  wall = 6, zw = 1, xw = 0.2,
                #                  main = "patient allocation",
                #                  yticks = round(seq(1,max(objectPlot$nptsdose),
                #                                 max(objectPlot$nptsdose)/7),0),
                #                  zticks = seq(1:nrow(objectPlot$nptsdose)),
                #                  xticks = seq(1:ncol(objectPlot$nptsdose)),
                #                  xlab = "Drug B",
                #                  ylab = "number of patients",
                #                  zlab = "Drug A",
                #                  col = 'lavender', alpha = 0.5)

                yv <- as.vector(objectPlot[["nptsdose"]])
                # drug A
                Drug_A <- rep(c(1: nrow(objectPlot[["nptsdose"]])),ncol(objectPlot$nptsdose))
                # drug B
                Drug_B <- rep(c(1: ncol(objectPlot[["nptsdose"]])),each=nrow(objectPlot$nptsdose))

                df <- as.data.frame(cbind(yv, Drug_A, Drug_B))

                print(ggplot(df, aes(Drug_A, yv)) +
                        ylab('number of patients')  +
                        xlab("Drug A")+
                        labs(title="patient allocation")+
                        geom_bar(stat="identity") +
                        facet_wrap(~ Drug_B,ncol=ncol(objectPlot$nptsdose), labeller = label_both)+
                        theme(strip.text.x = element_text(size = 10, angle = 0)))

                ## why not "number of toxicities" too?

                yv <- as.vector(objectPlot[["ntoxdose"]])
                # drug A
                Drug_A <- rep(c(1: nrow(objectPlot[["ntoxdose"]])),ncol(objectPlot$ntoxdose))
                # drug B
                Drug_B <- rep(c(1: ncol(objectPlot[["ntoxdose"]])),each=nrow(objectPlot$ntoxdose))

                df <- as.data.frame(cbind(yv, Drug_A, Drug_B))

                print(ggplot(df, aes(Drug_A, yv)) +
                        ylab('number of toxicities')  +
                        xlab("Drug A")+
                        labs(title="number of toxicities")+
                        geom_bar(stat="identity") +
                        facet_wrap(~ Drug_B,ncol=ncol(objectPlot$ntoxdose), labeller = label_both)+
                        theme(strip.text.x = element_text(size = 10, angle = 0)))

            }
            else {
                warning(paste("The variable (", strpattern,") to be plotted cannot be found in data object: ",
                        new.obj[1],".\n",sep=""))
            }
        }
    }

    ######################################################################################
    ## Plotting for select.mtd.kb() or select.mtd.comb.kb()
    ## [the latter function only finds a single MTD, not a MTD contour]
    else if (!is.null(objectPlot$MTD)) {
        if (objectPlot$MTD[1] == 99) {
            stop("All tested doses are overly toxic. No MTD is selected!\n")
            ## edit: changed from warning()
        }
        else {
            if (!is.null(objectPlot$p_est)) {
                # dev.new(); op <- par(no.readonly = TRUE); dev.off(); par(op)
                opar <- par(no.readonly = TRUE)
                on.exit(par(opar))
                par(mfrow = c(1,1), mar = c(5,5,4,2))

                ## Note that this if-else is switched from that of the BOIN code.
                if (length(objectPlot$MTD) < 2) { ## select.mtd.kb()
                    p_est = objectPlot$p_est; p_hat = p_est[,2]
                    ci = p_est[,3]; ci = gsub("[\\(\\)]", "", ci)
                    conf.intv = matrix(unlist(strsplit(ci, ",")), byrow = TRUE, ncol = 2)
                    if (p_est[1,2] == "----") {
                        stop("The trial is stopped since the lowest dose is too toxic.\n")
                        ## edit: changed from warning()
                    }
                    else {
                        numbs = ifelse(sum(p_hat == "----") == 0,
                                       length(p_hat), min(which(p_hat == "----")) - 1)
                        numbs2 = length(p_hat)
                        phatx = as.numeric(as.character(p_hat[1:numbs]))
                        lwr = as.numeric(as.character(conf.intv[1:numbs,1]))
                        upr = as.numeric(as.character(conf.intv[1:numbs,2]))
                        opar <- par(no.readonly = TRUE)
                        on.exit(par(opar))
                        par(mar = c(5,5,4,2))
                        plot(1:numbs2, ylim = c(0,1), xlab = 'dose level',
                             ylab = 'DLT rate', pch = '', xaxt = "n", cex.lab = 1.3)
                        axis(1, at = 1:numbs2, labels = 1:numbs2)
                        abline(h = objectPlot$target, lty = 2, col = 2)
                        points(1:numbs, phatx, pch = 19)
                        arrows(x0 = 1:numbs, x1 = 1:numbs, y0 = lwr, y1 = upr,
                               code = 3, angle = 90, length = 0.1)
                        if (numbs < numbs2) {
                            points((numbs + 1):numbs2,
                                    seq(min(1, max(phatx, na.rm = T) + 0.05),
                                        min(max(phatx, na.rm = T) + 0.20, 1),
                                        length = numbs2 - numbs),
                                    pch = "*", cex = 1.5)
                            legend('topleft', '*   no patient treated')
                        }
                    }
                }
                else{ ## select.mtd.comb.kb()
             #        if (!requireNamespace("epade", quietly = TRUE)) {
             #                    stop("Package \"epade\" needed for this function to work.",
             # "Please install it.", call. = FALSE)
             #
             #        }
             #
             #        requireNamespace("epade", quietly = TRUE)
                    # epade::bar3d.ade(objectPlot$p_est,
                    #                  wall = 6, zw = 1, xw = 0.2,
                    #                  yticks = round(seq(0,1,0.1),1),
                    #                  zticks = seq(1:nrow(objectPlot$p_est)),
                    #                  xticks = seq(1:ncol(objectPlot$p_est)),
                    #                  xlab = "Drug B",
                    #                  ylab = "DLT rate",
                    #                  zlab = "Drug A",
                    #                  col = 'lavender', alpha = 0.5)


                  yv <- as.vector(objectPlot[["p_est"]])
                  # drug A
                  Drug_A <- rep(c(1: nrow(objectPlot[["p_est"]])),ncol(objectPlot$p_est))
                  # drug B
                  Drug_B <- rep(c(1: ncol(objectPlot[["p_est"]])),each=nrow(objectPlot$p_est))

                  df <- as.data.frame(cbind(yv, Drug_A, Drug_B))

                  print(ggplot(df, aes(Drug_A, yv)) +
                          ylab('DLT rate')  +
                          xlab("Drug A")+
                          geom_bar(stat="identity") +
                          facet_wrap(~ Drug_B,ncol=ncol(objectPlot$p_est), labeller = label_both)+
                          theme(strip.text.x = element_text(size = 10, angle = 0)))


                }
            }
        }
    }

    ######################################################################################
        ## plotting for get.oc.obd.kb()
    else if (!is.null(objectPlot$name) && objectPlot$name == "get.oc.obd.kb") {
        # dev.new(); op <- par(no.readonly = TRUE); dev.off(); par(op)
        opar <- par(no.readonly = TRUE)
        on.exit(par(opar))
        par(mfrow = c(1,1), mar = c(5,5,4,2))

        if (length(new.obj) == 2) { ## WHEN TRUE??? Shouldn't it be length(x)?
            if (strpattern == "selpercent1") {
                bplot = barplot(objectPlot$selpercent1,
                                ylab = "selection percentage (%) using utility function 1", ylim = c(0,100),
                                cex.names = 1, xaxt = "n", cex.lab = 1.3)
                # mtext("selection percentage using utility function 1", 3, line = 2, cex = 1.3)
                axis(1, at = bplot, labels = seq(1, length(objectPlot$selpercent1)))
                mtext("dose level", 1, line = 2, cex = 1.3)
            }

            else if (strpattern == "selpercent2") {
                bplot = barplot(objectPlot$selpercent2,
                                ylab = "selection percentage (%) using utility function 2", ylim = c(0,100),
                                cex.names = 1, xaxt = "n", cex.lab = 1.3)
                axis(1, at = bplot, labels = seq(1, length(objectPlot$selpercent2)))
                mtext("dose level", 1, line = 2, cex = 1.3)
            }
            else if (strpattern == "selpercent3") {
                bplot = barplot(objectPlot$selpercent3,
                                ylab = "selection percentage (%) using utility function 3", ylim = c(0,100),
                                cex.names = 1, xaxt = "n", cex.lab = 1.3)
                axis(1, at = bplot, labels = seq(1, length(objectPlot$selpercent3)))
                mtext("dose level", 1, line = 2, cex = 1.3)
            }
            else if (strpattern == "npatients") {
                bplot = barplot(objectPlot$npatients,
                                ylab = "number of patients", ylim = c(0, max(objectPlot$npatients)),
                                cex.names = 1, beside = FALSE, xaxt = "n", cex.lab = 1.3)
                axis(1, at = bplot, labels = seq(1, length(objectPlot$npatients)))
                mtext("dose level", 1, line = 2, cex = 1.3)
            }
            else if (strpattern == "ntox") {
                bplot = barplot(objectPlot$ntox,
                                ylab = "number of toxicities", ylim = c(0, max(objectPlot$ntox)),
                                cex.names = 1, beside = FALSE, xaxt = "n", cex.lab = 1.3)
                axis(1, at = bplot, labels = seq(1, length(objectPlot$ntox)))
                mtext("dose level", 1, line = 2, cex = 1.3)
            }
            else if (strpattern == "neff") {
                bplot = barplot(objectPlot$neff,
                                ylab = "number of efficacies", ylim = c(0, max(objectPlot$ntox)),
                                cex.names = 1, beside = FALSE, xaxt = "n", cex.lab = 1.3)
                axis(1, at = bplot, labels = seq(1, length(objectPlot$neff)))
                mtext("dose level", 1, line = 2, cex = 1.3)
            }
        }
        else {
            par(mfrow = c(3,1), mar = c(5,5,4,2))
            bplot = barplot(objectPlot$selpercent1,
                            ylab = "selection percentage (%) ", ylim = c(0,100),
                            cex.names = 1, xaxt = "n", cex.lab = 1.3)
            mtext("selection percentage using utility function 1", 3, line = 2, cex = 1.3)
            axis(1, at = bplot, labels = seq(1, length(objectPlot$selpercent1)))
            mtext("dose level", 1, line = 2, cex = 1.3)

            bplot = barplot(objectPlot$selpercent2,
                            ylab = "selection percentage (%) ", ylim = c(0,100),
                            cex.names = 1, xaxt = "n", cex.lab = 1.3)
            mtext("selection percentage using utility function 2", 3, line = 2, cex = 1.3)
            axis(1, at = bplot, labels = seq(1, length(objectPlot$selpercent2)))
            mtext("dose level", 1, line = 2, cex = 1.3)

            bplot = barplot(objectPlot$selpercent3,
                            ylab = "selection percentage (%) ", ylim = c(0,100),
                            cex.names = 1, xaxt = "n", cex.lab = 1.3)
            mtext("selection percentage using utility function 3", 3, line = 2, cex = 1.3)
            axis(1, at = bplot, labels = seq(1, length(objectPlot$selpercent3)))
            mtext("dose level", 1, line = 2, cex = 1.3)

            bplot = barplot(objectPlot$npatients,
                            ylab = "number of patients", ylim = c(0, max(objectPlot$npatients)),
                            cex.names = 1, beside = FALSE, xaxt = "n", cex.lab = 1.3)
            axis(1, at = bplot, labels = seq(1, length(objectPlot$npatients)))
            mtext("patient allocation", 3, line = 2, cex = 1.3)
            mtext("dose level", 1, line = 2, cex = 1.3)

            ## number of toxicities
            bplot = barplot(objectPlot$ntox,
                            ylab = "number of toxicities",
                            ylim = c(0, max(objectPlot$ntox)),
                            cex.names = 1, beside = FALSE, xaxt = "n", cex.lab = 1.3)
            axis(1, at = bplot, labels = seq(1, length(objectPlot$ntox)))
            mtext("patient allocation for toxicities", 3, line = 2, cex = 1.3)
            mtext("dose level", 1, line = 2, cex = 1.3)

            bplot = barplot(objectPlot$neff,
                            ylab = "number of efficacies",
                            ylim = c(0, max(objectPlot$neff)),
                            cex.names = 1, beside = FALSE, xaxt = "n", cex.lab = 1.3)
            axis(1, at = bplot, labels = seq(1, length(objectPlot$ntox)))
            mtext("patient allocation for efficacies", 3, line = 2, cex = 1.3)
            mtext("dose level", 1, line = 2, cex = 1.3)
        }
    }


    else {
        stop("Please double check and specify the variable to be plotted...\n")

    }

    par(opar)
}