R/sensitivityplot.R

In aursiber/DRomics: Dose Response for Omics

# Sensitivity plot
# Plot of a summary of BMD values per group of items 
# (e.g. from biological annotation), 
# and optionally per molecular level 
# (or per another additional grouping level )
sensitivityplot <- function(extendedres, BMDtype = c("zSD", "xfold"),
                            group, ECDF_plot = TRUE, colorby,
                            BMDsummary = c("first.quartile", "median" , "median.and.IQR"),
                            BMD_log_transfo = TRUE,
                            line.size = 0.5, line.alpha = 0.5, point.alpha = 0.5)
{
  BMDtype <- match.arg(BMDtype, c("zSD", "xfold"))
  BMDsummary <- match.arg(BMDsummary, c("first.quartile", "median", "median.and.IQR"))
  
  if (missing(extendedres) | !is.data.frame(extendedres))
    stop("The first argument of sensitivityplot must be a dataframe 
    (see ?sensitivityplot for details).")
  
  cnames <- colnames(extendedres)
  
  if (BMDtype == "zSD")
  {  
    if (any(!is.element(c("BMD.zSD"), cnames)))
      stop("The first argument of sensitivityplot must be a dataframe
      containing a column named BMD.zSD and other columns coding for group of items.")
    variable <- extendedres[, "BMD.zSD"]
  } else 
  {
    if (any(!is.element(c("BMD.xfold"), cnames)))
      stop("The first argument of sensitivityplot must be a dataframe
      containing a column named BMD.xfold and other columns coding for groups of items.")
    variable <- extendedres[, "BMD.xfold"]
  }
  
  if (!is.character(group)) 
    stop("group should be a character string for the name of the column defining groups.")
  if (!is.element(group, cnames))
    stop("group should be a character string corresponding to the name of a column of
           extendedres, the dataframe given in input.")
  groupby <- as.factor(extendedres[, group])
  
  if (!missing(colorby))
  {
    if (!is.character(colorby)) 
      stop("colorby should be a character string for the name of the column coding for the point color.")
    if (!is.element(colorby, cnames))
      stop("colorby should be a character string corresponding to the name of a column of
           extendedres, the dataframe given in input.")
  }
  
  
   firstquartilefun <- function(x) stats::quantile(x, probs = 0.25, na.rm = TRUE)
   secondquartilefun <- function(x) stats::quantile(x, probs = 0.5, na.rm = TRUE)
   thirdquartilefun <- function(x) stats::quantile(x, probs = 0.75, na.rm = TRUE)
  
  if (missing(colorby))
  {
    dnb <- as.data.frame(table(groupby))
    colnames(dnb) <- c("groupby", "nb_of_items")
    dnb$firstquartile <- tapply(variable, groupby, firstquartilefun)
    dnb$secondquartile <- tapply(variable, groupby, secondquartilefun)
    dnb$thirdquartile <- tapply(variable, groupby, thirdquartilefun)

  } else
  {
    level <- extendedres[, colorby]
    dnb <- as.data.frame(table(groupby, level))
    colnames(dnb) <- c("groupby", "level","nb_of_items")
    dnb$firstquartile <- tapply(variable, level:groupby, firstquartilefun)
    dnb$secondquartile <- tapply(variable, level:groupby, secondquartilefun)
    dnb$thirdquartile <- tapply(variable, level:groupby, thirdquartilefun)
  }
  dnb <- dnb[dnb$nb_of_items != 0, ]
  
  if (!missing(colorby)) ECDF_plot <- FALSE
  
  if (ECDF_plot)
  {
    # order by chosen summary
    if (BMDsummary == "median" | BMDsummary == "median.and.IQR")
      dnb <- dnb[order(dnb$secondquartile), ] else 
    if (BMDsummary == "first.quartile")
      dnb <- dnb[order(dnb$firstquartile), ]
          
    # fix the order of the modalities of by as in the ordered data set
    dnb$groupby <- factor(dnb$groupby, levels = dnb$groupby)
  }
  
  if (BMDsummary == "first.quartile")
  {
    if (missing(colorby))
    {
      gg <- ggplot(dnb, aes(x = .data$groupby, y = .data$firstquartile, 
                             size = .data$nb_of_items))
    } else
    {  
      gg <- ggplot(dnb, aes(x = .data$groupby, y = .data$firstquartile, 
                             color = .data$level, 
                             size = .data$nb_of_items))
    }
    
    gg <- gg + geom_point(stat = 'identity', alpha = point.alpha)  +  coord_flip() +
      labs(x = "", y = "BMD 25th quantiles") 
    
  } else {
    if (missing(colorby))
    {
      gg <- ggplot(dnb, aes(x = .data$groupby, y = .data$secondquartile, 
                             size = .data$nb_of_items)) 
    } else {
      if (BMDsummary == "median") {
        gg <- ggplot(dnb, aes(x = .data$groupby, y = .data$secondquartile, 
                               color = .data$level, 
                               size = .data$nb_of_items))
      } else {
        gg <- ggplot(dnb, aes(x = .data$groupby, y = .data$secondquartile, 
                               color = .data$level, 
                               size = .data$nb_of_items))
      }
    }
    gg <- gg + geom_point(stat = 'identity', alpha = point.alpha) + coord_flip() 
    
    if (BMDsummary == "median")
    {
      gg <- gg + labs(x = "", y = "BMD medians") 
    } else {
      gg <- gg + geom_errorbar(aes(ymin = .data$firstquartile, 
                                    ymax = .data$thirdquartile),
                               linewidth = line.size, 
                               alpha = line.alpha,
                               width = 0) +
        # line to remove lines on the size legend
        guides(size = guide_legend(override.aes = list(linetype = 0))) +
        labs(x = "", y = "BMD medians and IQRs")
      }
  }    
  if (BMD_log_transfo)
    gg <- gg + scale_y_log10()
  
  if (!missing(colorby))
  {
    gg <- gg + labs(color = colorby)
  }
  
  round.quartiles.minmax <- unique(round(quantile(dnb$nb_of_items, probs = c(0, 0.25, 0.5, 0.75, 1))))
  gg <- gg + scale_size_continuous(breaks = as.numeric(round.quartiles.minmax)) + 
            labs(size = "nb. of items")
  return(gg)
}