R/plot.trtsel.R

#' plot risk curves, treatment effect curves or cdf of risk for a trtsel
#' object.
#' 
#' Plots a treatment selection marker.  Estimated "risk curves" and "treatment
#' effect curves" are plotted.  An object of class "trtsel" must first be
#' created using the function "trtsel" by supplying a data.frame containing
#' marker, treatment, and event status information.
#' 
#' 
#' @aliases plot.trtsel plot
#' @param x An object of class "trtsel", created by using the function
#' "trtsel."
#' @param bootstraps Number of bootstrap replicates for creating confidence
#' intervals and bands. Default value is 500. Set bootstraps=0 if no confidence
#' intervals or bands are desired.
#' @param plot.type Type of plot to produce. Options are "risk" (default) for
#' risk curves which show risk of adverse event as a function of marker and
#' treatment, "treatment effect" for the treatment effect curve which shows the
#' difference in marker-specific risk of adverse event with vs. without
#' treatment, "cdf" for the cumulative distribution function of the treatment
#' effect, or "selection impact" to show a plot of the estimated event rate if
#' different proportions of the population are treated.  Cdf plots cannot be
#' generated for a discrete marker.
#' @param ci Indication of whether horizontal or vertical confidence intervals
#' are desired. Character string of either "horizontal" (default) or
#' "vertical." For a discrete marker only vertical confidence bands can be
#' produced. See note below for more information.
#' @param alpha (1-alpha)*100\% confidence intervals are calculated. Default
#' value is alpha = 0.05 which yields 95\% CI's.
#' @param fixed.values A numeric vector indicating fixed values on the x- or
#' y-axes at which bootstrap confidence intervals are provided. If
#' "fixed.values" are provided, point-wise confidence intervals will be printed
#' (i.e. conf.bands will be taken as FALSE). Only applies to a continuous
#' marker.
#' @param offset If confidence intervals are to be plotted for specified
#' fixed.values, offset is the amount of distance to offset confidence
#' intervals so that they do not overlap on the plot. The default value is
#' 0.01.  Only applies to a continuous marker.
#' @param conf.bands Indication of whether pointwise confidence intervals are
#' shown for the curve(s).  TRUE (default) or FALSE. If "fixed.values" are
#' input, this option is ignored and no confidence bands are produced.
#' @param conf.bandsN If conf.bands = TRUE, the number of points along the x-
#' or y-axis at which to calculate the pointwise confidence intervals. The
#' default is 100.
#' @param trt.names A vector of length 2 indicating the names for the two
#' treatment options, T= 1 and T = 0, respectively, for the plot legend. This
#' option is only used when plot.type="risk". The default value is
#' c("Treatment", "No Treatment").
#' @param xlab A label for the x-axis. Default values depend on plot.type.
#' @param ylab A label for the y-axis. Default values depend on plot.type.
#' @param xlim The limits for the x-axisof the plot, in the form c(x1,x2)
#' @param ylim The limits for the y-axis of the plot, in the form c(y1,y2)
#' @param main The main title for the plot.
#' @param show.marker.axis For risk curve plots with a continuous marker,
#' should an axis showing marker values be displayed? default is TRUE.
#' @param \dots N/A
#' @return
#' 
#' Generates a plot. In addition, if the function call is assigned to a
#' variable (as in "myplot <- plot(mytrtsel)"), plot.trtsel returns a list
#' composed of: \item{plot}{ ggplot object containing the printed plot. }
#' \item{ci.bounds}{ A data.frame containing the bounds of the bootstrap-based
#' confidence intervals, along with the fixed.values they are centered around,
#' if applicable. }
#' @note For a discrete marker, vertical confidence intervals are generated by
#' fixing the marker value and generating a bootstrap distribution for risk
#' given marker. Bands are then generated based on the percentiles of the
#' empirical bootstrap distribution. For a continuous marker, vertical
#' confidence intervals are generated by fixing the empirical marker quantiles,
#' and generating the empirical distribution of risk given marker percentile.
#' @seealso \code{\link{trtsel}} for creating trtsel objects,
#' \code{\link{evaluate.trtsel}} for evaluating marker performance,
#' \code{\link{calibrate.trtsel}} for assessing model calibration, and
#' \code{\link{compare.trtsel}} to compare two trtsel object.

#' @examples
#' 
#' 
#' data(tsdata)
#' 
#' ###########################
#' ## Create trtsel objects 
#' ###########################
#' 
#' trtsel.Y1 <- trtsel(event ~ Y1*trt, 
#'                    treatment.name = "trt", 
#'                    data = tsdata, 
#'                    study.design = "RCT",
#'                    default.trt = "trt all")
#'
#' 
#' 
#' 
#' ##########################
#' ## Use the plot function 
#' ##########################
#' 
#' # Plot risk curves
#' plot(trtsel.Y1, main = "Marker Y1", 
#'      plot.type = "risk", bootstraps = 10, #set higher in practice
#'      trt.names=c("trt","no trt"))
#'      
#' # Now with confidence intervals around fixed.values 
#'  plot(trtsel.Y1, main = "Marker Y1",
#'                         plot.type = "risk", ci = "horizontal", 
#'                         fixed.values = c(.2, .5), 
#'                         offset = .01, bootstraps = 10,
#'                         trt.names=c("trt","no trt"))
#'                         
#' # Plot treatment effect curves
#' #with confidence intervals around fixed.values
#' plot(trtsel.Y1, main = "Marker Y1", 
#'                        plot.type = "treatment effect",
#'                        ci = "vertical", conf.bands = FALSE,
#'                        fixed.values = c(10, 20), bootstraps = 10)
#' 
#' @import ggplot2
#' @import grid 
#' @import graphics
#' @method plot trtsel
#' @export


plot.trtsel <-
function(x, bootstraps = 500,  
            plot.type = c( "treatment effect","risk", "cdf", "selection impact"), 
            ci = c("default", "horizontal", "vertical", "none"), 
            alpha = .05, 
            fixed.values = NULL,  
            offset = 0.01, 
            conf.bands = TRUE, 
            conf.bandsN = 100, 
            trt.names = c("Treatment", "    No \nTreatment"),
            xlab = NULL, ylab = NULL, xlim = NULL, ylim = NULL, main = NULL, show.marker.axis = TRUE, 
            ...)

{
 
  if(!is.trtsel(x)) stop("x must be an object of class 'trtsel' created by using the function 'trtsel' see ?trtsel for more help")
  
  plot.type <- match.arg(plot.type) 
  
  if(!is.element(plot.type, c("risk", "treatment effect", "cdf", "selection impact"))){ 
    stop("plot.type must be one of  \"treatment effect\", \"risk\", \"cdf\" or \"selection impact\"")
  }
  stopifnot(length(plot.type) ==1)
  ci <- match.arg(ci)

  if(!is.element(ci, c("default", "horizontal", "vertical", "none"))){ 
    stop("ci must be one of \"default\",  \"horizontal\", or  \"vertical\" or \"none\" ")
  }
  

  if(alpha<0 | alpha > 1) stop("Error: alpha should be between 0 and 1")
  

  if(bootstraps < 2) warning("Number of bootstraps must be greater than 1, bootstrap confidence intervals will not be computed") 
  if(ci == "none") bootstraps = 0; 
  # theta curves can only be plotted for cohort data and with vertical ci's
  if(substr(plot.type, 1, 3)=="sel"){
   # if(substr(x$model.fit$study.design, 1, 3) != "ran") stop("theta curves cannot be created for subcohort designs")
    if(substr(ci, 1, 1) =="h") {
      
      warning("horizontal ci's are not available for selection impact curves, vertical ci's will be calculated")
      ci = "vertical"
    }
  }
  #set default ci's 

  if(ci =="default"){
    #continuous marker
    if(is.null(x$model.fit$disc.rec.no.trt)){
      if(substr(plot.type, 1, 3) =="ris") ci = "horizontal"
      if(substr(plot.type, 1, 3) =="tre") ci = "horizontal"
      if(substr(plot.type, 1, 3) =="cdf") ci = "vertical"
      if(substr(plot.type, 1, 3) =="sel") ci = "vertical"
    }else{
      
      if(substr(plot.type, 1, 3) =="ris") ci = "vertical"
      if(substr(plot.type, 1, 3) =="tre") ci = "vertical"
      if(substr(plot.type, 1, 3) =="cdf") ci = "horizontal"
    }
    print(paste(ci, "confidence intervals will be plotted."))

  }
  #no cdf plots for binary marker
  
  if(!is.null(x$model.fit$disc.rec.no.trt)){
  if(is.element(substr(plot.type, 1, 3), c("cdf", "sel"))) stop("cdf or selection impact plots cannot be created for a binary marker. Please choose plot.type to be \"risk\" or \"treatment effect\" ")
  }
  #save the current plot parameters
  #old.par <- par(no.readonly = TRUE)
 

#extract the needed data from x, which is our TrtSel object


  rho  <- x$model.fit$cohort.attributes
  if(length(rho) == 4) rho = c(rho, -9999, -9999, -9999) #accomodate the nested case-control design
  
  study.design <- x$model.fit$study.design
  link <- x$model.fit$family$family
  boot.sample <- x$functions$boot.sample
  get.F <- x$functions$get.F
  delta <- x$derived.data$trt.effect  

  if(link == "risks_provided") {
    marker = NULL
    show.marker.axis = FALSE
  }
  

 if(is.null(x$model.fit$disc.rec.no.trt)){  #continuous marker 
    plot.functions <- list(  predcurvePLOT_gg, trteffectPLOT_gg, CDFdeltaPLOT_gg, SelectionImpactPLOT_gg)
    
  if(length(fixed.values)!=0) conf.bands = FALSE
  if(conf.bands & length(fixed.values)==0 ){
   #create fixed values

   if(substr(ci, 1,1 )=="v"){
       
      if(is.element(substr(plot.type, 1,3), c("tre", "ris", "sel"))) fixed.values = seq(from = 1, to = 100, length.out = conf.bandsN)
      else if(substr(plot.type, 1,3)=="cdf") fixed.values = seq(from = min(delta), to = max(delta), length.out = conf.bandsN)
      
   }else{
      
      if(is.element(substr(plot.type, 1,3), c("cdf"))) fixed.values = seq(from = 1, to = 100, length.out = conf.bandsN)
      else if(substr(plot.type, 1,3)=="tre") fixed.values = seq(from = min(delta), to = max(delta), length.out = conf.bandsN)
      else if(substr(plot.type, 1,3)=="ris"){
       allrisks <- c(x$derived.data$fittedrisk.t0, 
                     x$derived.data$fittedrisk.t1)
       fixed.values = seq(from = min(allrisks), to = max(allrisks), length.out = conf.bandsN)

      }
      ## add case for theta curve

   }
   offset = 0

   }
    ##Bootstrap for confidence intervals...
    #bootstrapping done by fixing response and strapping marker 
    
    if((conf.bands & bootstraps>1) | (length(fixed.values)>0 & bootstraps > 1)){ 

      ci.bounds <- get.plot.ci(x, 
                               plot.type, 
                               ci, bootstraps, 
                               fixed.values =fixed.values,  
                               alpha = alpha)
    }else{ 
      ci.bounds <- NULL
      conf.bands = FALSE
    }
    ## end Bootstrap
    
    
    
    
  }else{
    #discrete marker...ignore fixed values and ci bands, we only need ci's around the observed points
    
    plot.functions <- list(  predcurvePLOT_gg_disc, trteffectPLOT_gg_disc, CDFdeltaPLOT_gg_disc)

    ##Bootstrap for confidence intervals...
    #much simpler for a discrete marker

    if(( conf.bands & bootstraps>1) ){ 
      
      ci.bounds <- get.plot.ci_disc(x, plot.type, 
                                    ci, bootstraps, 
                                    alpha)
      ci = ci.bounds$newci
      ci.bounds = ci.bounds$myconf.ints
      }else{ 
      ci.bounds <- NULL
      conf.bands = FALSE
    }
    ## end Bootstrap
    
    
  }

  tmp.plotfun <- plot.functions[[match(plot.type, c("risk", "treatment effect", "cdf", "selection impact"))]]
   

if(is.null(x$model.fit$disc.rec.no.trt)){  # continous marker 
  if(substring(plot.type, 1, 4) == "risk"){

  curves <- tmp.plotfun(x, ci, ci.bounds, get.F, fixed.values, conf.bands,  rho, trt.names, xlab, ylab, xlim, ylim, main, show.marker.axis, offset = offset)

    if(!is.null(ci.bounds)){
      ci.bounds <- data.frame(t(ci.bounds))
      ci.bounds <- cbind(fixed.values, ci.bounds)
      names(ci.bounds) <- c("fixed.values", "trt0.lower", "trt0.upper", "trt1.lower", "trt1.upper")
      }
    }else{

      curves <- tmp.plotfun(x, ci, ci.bounds, get.F, fixed.values, conf.bands, rho, xlab, ylab, xlim, ylim, main)
      if(!is.null(ci.bounds)){
        ci.bounds <- data.frame(t(ci.bounds))
        ci.bounds <- cbind(fixed.values, ci.bounds)
        names(ci.bounds) <- c("fixed.values", "lower", "upper")
      }
    }
  
}else{  # discrete marker
  

    
    curves <- tmp.plotfun(x=x, ci = ci, ci.bounds = ci.bounds, get.F = get.F, 
                          xlab = xlab, ylab=ylab, 
                          xlim = xlim, ylim=ylim, 
                          main=main, trt.names = trt.names)

    if(!is.null(ci.bounds)){
      ci.bounds <- data.frame(curves[[2]])

    }
   
}    

  #par(old.par)

invisible(list("plot" = curves[[1]], "ci.bounds" = ci.bounds))
}
mdbrown/TreatmentSelection documentation built on May 22, 2019, 3:23 p.m.