R/vioplotGP.R

Defines functions vioplotGP

Documented in vioplotGP

#' Violin plots for Groups
#'
#' Produce a Violin plot according to groups taking into account weightings and compute comparison tests
#'
#'
#' @param x a numeric vector of data
#' @param gp factors with the same length as x describing groups
#' @param weights a numeric vector with the same length as x
#' @param labels group names in the same order than \code{gp} levels
#' @param xlab label for the x-axis
#' @param xlim, ylim range of values used for x and y-axes
#' @param at position of each group in the y-axis
#' @param wex the maximum value for the density curve
#' @param h the kernel density smoothing parameter
#' @param nblim minimum length of x for computing the kernel density
#' @param Test a character value for group comparison tests "tukey" or "wilcox"
#' @param adj.Test,side.Test,line.Test position paramters for the test value and P.value
#' @param Fname ?
#' @param type.letters if Test is done, "latin" or "greek" letters are used to show significant differences between groups
#' @param pos.letters,col.letters position (1, 2, 3 or 4) and colors for letters
#' @param add if TRUE adds the plot to the current plot window
#' @param col ?
#' @param pch ?
#' @param lwd ?
#' @param lty ?
#' @param cex ?
#' @param cex.points ?
#' @param col.points ?
#' @param bg ?
#' @param las ?
#' @param boxplot logical value for adding a boxplot or not
#' @param side a character value describing if the distribution should be drawn
#'             in the "above", "below", "right", "left" or "both" side(s)
#' @param fill the color of the distribution area
#' @param pts logical value for adding x data on the distribution curve
#' @param boxplot.fill the color of the boxplot
#' @param horizontal logical value
#'
#' @examples
#'
#' data("ChickWeight")
#' vioplotGP(x = chickwts$weight, gp = chickwts$feed, las = 2,
#'           horizontal = FALSE, side = "both", wex = .4, fill = "grey70",
#'           boxplot.fill = "grey30", xlab = "Chick weights (mg)",
#'           Test = "tukey", pos.letters = 3, col.letters = "white",
#'           col.points = "grey50")
#'
#' @seealso
#'
#' For group comparison see vioplotGp
#'
#' @import graphics stats
#' @importFrom multcomp glht mcp
#'
#' @export

vioplotGP <- function(x, gp, weights = NULL, labels = NA, xlab = NA, xlim = NA,
                    ylim = NA, at = NA, wex =.75, h = NA, nblim = 10, Test = NA,
                    adj.Test = 0.01, line.Test = NA, side.Test=3, Fname=NA, type.letters="latin",
                    pos.letters = 1, col.letters = 1, add = F, col = 1, pch = 19,
                    lwd = 1, lty = 1, cex = 1, cex.points = 1, col.points = 1,
                    bg = 1, las = 0, boxplot = TRUE, side = "above", fill = NA,
                    pts = TRUE, boxplot.fill = "white", horizontal = TRUE){

  ##### check ####

  if (!is.na(Test) & !(Test%in%c("wilcox","tukey"))) {
    stop("Test argument should be 'wilcox' or 'tukey'")
  }
  if (!(type.letters%in%c("latin","greek"))) {
    stop("Test argument should be 'latin' or 'greek'")
  }
  if (!side %in% c("above", "below", "right", "left", "both")) {
    stop('side argument should be "above","below", "right", "left" or "both"')
  }
  if (length(x) != length(gp)) {
    stop("x and gp are not the same length")
  }

  #### remove NA values ####
  gp <- as.factor(gp[!is.na(x)])
  if(!is.null(weights)) weights <- weights[!is.na(x)]
  x <- x[!is.na(x)]

  #### data descriptors ####
  ngp <- nlevels(gp)
  xgp <- sapply(X = levels(gp),FUN = function(X){x[gp == X]},simplify = F)
  if(class(xgp)!="list") xgp <- list(xgp)
  quart <- lapply(xgp,quantile)

  #### weights ####
  w <- weights
  n <- length(x)
  if (is.null(w)) {
     w <- rep(1,n)
  } else if (length(w)!=n) {
    stop("Weight length not correct")
  }

  w <- w/sum(w)
  wgp <- sapply(X =levels(gp), FUN = function(X){w[gp==X]/sum(w[gp==X])}, simplify = F )
  if(class(wgp)!="list") wgp <- list(wgp)

  #### at ####
  if (all(is.na(at))) {
    at <- 1:ngp
  } else if (length(at)!=ngp) {
    stop("at length not correct")
  }

  ##### labels ####
  if (all(is.na(labels))) {
    labels <- levels(as.factor(gp))
  }

  #### if test ####
  testdone <- FALSE

  if (!is.na(Test) & Test == "wilcox"){
    test <- pairwise.wilcox.test(x, gp, p.adjust.methods = "bon")$p.value

    p.values <- cbind(rbind(rep(NA, ngp - 1), test), rep(NA, ngp))
    colnames(p.values) <- levels(gp)
    rownames(p.values) <- levels(gp)

    for (i in 1:(ngp - 1)){
      p.values[i, (i + 1):ngp] <- p.values[(i + 1):ngp, i]
    }

    diag(p.values) <- 1
    stat <- ""
    testdone <- TRUE
  } else if (!is.na(Test) & Test == "tukey") {

    test <- summary(glht(lm(x ~ gp, weights = w), mcp(gp = "Tukey")))$test

    comp <- as.numeric(test$pvalues)
    compnames <- do.call(cbind, strsplit(names(test$coefficients)," - "))
    p.values <- matrix(1, ngp, ngp)
    colnames(p.values) <- rownames(p.values) <- levels(gp)
    for (i in 1:length(comp)){
      p.values[compnames[1, i], compnames[2, i]] <- p.values[compnames[2, i],compnames[1, i]] <- comp[i]
    }
    Fval <- round(summary(aov(x ~ gp, weights = w))[[1]][1,"F value"],2)
    pval <- summary(aov(x ~ gp, weights = w))[[1]][1, "Pr(>F)"]
    pval <- ifelse(pval < 0.001, "***",
                   ifelse(pval < 0.01, "**",
                          ifelse(pval < 0.05, "*",
                                 ifelse(pval < 0.1, ".", "ns"))))

    Fname <- ifelse(is.na(Fname),"F",Fname)
    stat <- paste0(Fname, "==", Fval, "*plain('",pval,"')")
    testdone <- TRUE
  }
  if(testdone) {
    p.values[p.values > 0.05] <- 1
    p.values[p.values != 1] <- 0
    ordre <- order(tapply(x, gp, mean, na.rm = T))
    p.values <- p.values[ordre, ordre]
    p.values <- unique(p.values)
    p.values <- p.values[apply(p.values, MARGIN = 1, function(x) any(x == 0)),]

    if(nrow(p.values) == 0){
      lettres <- rep("a", ngp)
      names(lettres) <- levels(gp)
    }else{
      lettres <- rep("", ngp)
      for (k in 1:nrow(p.values)){
        lettres[p.values[k,] == 1] <- paste0(lettres[p.values[k,] == 1], letters[k])
      }
      names(lettres) <- colnames(p.values)
      lettres <- lettres[match(levels(gp), names(lettres))]
    }

    if(type.letters=="greek") {
      greekletters <- c("alpha", "beta", "gamma", "delta", "epsilon","zeta","eta","theta", "iota", "kappa")
      lettres <- lapply(strsplit(lettres, ""), function(X){
        mapply(function(x) {which(letters == x)},X)
      })
      lettres <- lapply(lettres, function(X) { greekletters[X] })
      lettres <- unlist(lapply(lettres, function(X){
        do.call(paste, c(as.list(X), sep="*"))
      }))
    }
  }

  #plot window and par
  par0 <- par()["cex"]
  par(cex = cex)

  if (!add) {
    if (horizontal) {
      plot.new()
      if (any(is.na(xlim))) {
        xlimi <- c(min(x),max(x))
      } else {
        xlimi = xlim
      }

      if (any(is.na(ylim))) {
        ylimi <- c(min(at)-1, max(at)+1)
      } else {
        ylimi = ylim
      }

      plot.window(xlim = xlimi, ylim = ylimi)
      box(which ="plot", col = 1, lwd = 1)
      axis(1, las = las)
      if (any(labels != "none")) {
        axis(2, at = at, labels = labels, las = las)
      }
      if (all(is.na(xlab))) {
        title(xlab = deparse(substitute(x)))
      } else if (any(xlab!="none")) {
        title(xlab = xlab)
      }
    } else {
      plot.new()

      if (any(is.na(xlim))) {
        ylimi <- c(min(x),max(x))
      } else {
        ylimi = xlim
      }

      if (any(is.na(ylim))) {
        xlimi <- c(min(at) - 1, max(at) + 1)
      } else {
        xlimi = ylim
      }

      plot.window(xlim = xlimi, ylim = ylimi)
      box(which ="plot", col = 1, lwd = 1)
      axis(2, las = las)

      if (any(labels != "none")) {
        axis(1, at = at, labels = labels, las = las)
      }

      if (all(is.na(xlab))) {
        title(ylab = deparse(substitute(x)))
      } else if (any(xlab != "none")) {
        title(ylab = xlab)
      }
    }
  }

  # par for each vioplot
  wex <- rep(wex, length.out = ngp)
  h <- rep(h, length.out = ngp)
  col <- rep(col, length.out = ngp)
  pch <- rep(pch, length.out = ngp)
  lwd <- rep(lwd, length.out = ngp)
  lty <- rep(lty, length.out = ngp)
  cex.points <- rep(cex.points, length.out = ngp)
  col.points <- rep(col.points, length.out = ngp)
  bg <- rep(bg, length.out = ngp)
  side <- rep(side, length.out = ngp)
  fill <- rep(fill, length.out = ngp)
  pts <- rep(pts, length.out = ngp)
  boxplot.fill <- rep(boxplot.fill, length.out = ngp)

  # vioplots
  hi <- NULL
  x.lettre<-NULL
  for (i in 1:ngp) {
    hi <- c(hi, vioplot(xgp[[i]], weights = wgp[[i]], wex = wex[i], h = h[i],
                        add = T, nblim = nblim, at = at[i], pch = pch[i],
                        lwd = lwd[i], col = col[i], col.points = col.points[i],
                        lty = lty[i], cex = cex.points[i], bg = bg[i],
                        boxplot = boxplot, side = side[i], fill = fill[i],
                        pts = pts[i],boxplot.fill = boxplot.fill[i],
                        horizontal = horizontal))

    x.lettre <-c(x.lettre,quart[[i]][3])
  }

  if (testdone) {
    line <- ifelse(is.na(line.Test),ifelse(add, -2.2, -1.2),line.Test)
    mtext(side = side.Test, adj = adj.Test, text = parse(text=stat), line = line, cex = cex, las = 1)

    if (type.letters=="greek") {
      if (horizontal) {
        text(x.lettre, at, parse(text = lettres), pos = pos.letters,
        cex = cex, col = col.letters)
      } else {
        text(at, x.lettre, parse(text = lettres), pos = pos.letters,
            cex = cex, col = col.letters)
      }
    } else {
      if (horizontal) {
        text(x.lettre, at, lettres, pos = pos.letters, cex = cex,
              col = col.letters)
      } else {
        text(at, x.lettre, lettres, pos = pos.letters, cex = cex,
          col = col.letters)
      }
    }
  }
  names(hi) <- levels(gp)
  par(par0)
  return(hi)
}
Rekyt/divr documentation built on April 1, 2020, 5:35 a.m.