R/ecdfplotwithCI.R

In DRomics: Dose Response for Omics

ecdfplotwithCI <- function(variable, CI.lower, CI.upper, by, CI.col = "blue", CI.alpha = 1,
                           add.point = TRUE, point.size = 1, point.type = 16) {
  if (!is.numeric(variable))
    stop("Arguments variable, CI.lower and CI.upper must be numeric vectors of the same size.")
  leng <- length(variable)
  if (!is.numeric(CI.lower) || (length(CI.lower) != leng) ||
      !is.numeric(CI.upper) || (length(CI.upper) != leng))
    stop("Arguments variable, CI.lower and CI.upper must be numeric vectors of the same size.")
  d <- data.frame(variable = variable, lower = CI.lower, upper = CI.upper)
  
  if (!missing(by)) {
    if (length(by) != leng)
      stop("Argument by must be a factor of the same length as argument variable.")
    d$by <- as.factor(by)
  }
  if (is.factor(CI.col))
    d$CI.col <- CI.col
  if (is.factor(point.type))
    d$point.type <- point.type
  
  if (!missing(by)) {
    ntot <- nrow(d)
    uniqueby <- unique(d$by)
    n.uniqueby <- length(uniqueby)
    d$ECDF <- rep(0, ntot) # initialization
    for (i in 1:n.uniqueby) {
      indi <- which(d$by == uniqueby[i])
      ntoti <- length(indi)
      d$ECDF[indi] <- (rank(d$variable[indi], ties.method = "first") - 0.5) / ntoti
      # not strictly equivalent with ecdf (i / n)
      # d$ECDF[indi] <- ecdf(d$variable[indi])(d$variable[indi])
    }
    g <- ggplot(data = d, mapping = aes(x = variable, y = .data$ECDF)) + facet_wrap(~ by)
    
    if (is.factor(CI.col)) {
      g <- g +
        geom_errorbarh(aes(xmin = .data$lower, xmax = .data$upper, color = .data$CI.col),
                       alpha = CI.alpha, height = 0)
    } else {
      g <- g +
        geom_errorbarh(aes(xmin = .data$lower, xmax = .data$upper), color = CI.col,
                       alpha = CI.alpha, height = 0)
    }
  } else {
    ntot <- nrow(d)
    d$ECDF <- (rank(d$variable, ties.method = "first") - 0.5) / ntot
    # not strictly equivalent with ecdf (i / n)
    # d$variable <- ecdf(d$variable)(d$variable)
    g <- ggplot(data = d, mapping = aes(x = variable, y = .data$ECDF))
    if (is.factor(CI.col)) {
      g <- g +
        geom_errorbarh(aes(xmin = .data$lower, xmax = .data$upper, color = .data$CI.col),
                       alpha = CI.alpha, height = 0)
    } else {
      g <- g +
        geom_errorbarh(aes(xmin = .data$lower, xmax = .data$upper), color = CI.col,
                       alpha = CI.alpha, height = 0)
    }
  }
  if (add.point) {
    if (is.factor(point.type)) {
      g <- g + geom_point(aes(shape = .data$point.type),
                          size = point.size)
    } else {
      g <- g + geom_point(shape = point.type,
                          size = point.size)
    }
  }
  return(g)
}