R/BoxPlot.R

In CorReg: Linear Regression Based on Linear Structure Between Variables

#' Boxplot with confidence interval and ANOVA on the plot.
#' 
#' @description classical boxplot function improved with integrated confidence level on the mean for each group ploted on the graph and also ANOVA with p-value and its interpretation given in the legend.
#' 
#' @param num a numeric vector to plot boxplot(num~grp). Represents the value that will be compared between the groups.
#' @param grp a qualitative vector (factor) to plot boxplot(num~grp). Represents the groups we will compare.
#' @param data a data.frame (or list) from which the variables in formula should be taken.
#' @param AnoVa boolean to compute or not anova (when multiple groups) to see if they differ in mean.If false the Kruskal-Wallis Rank Sum test is computed instead.
#' @param risk the risk value used for confidence intervals. 
#' @param lang lingustic parameter to specify the language of the legend
#' @param verbose boolean to make a test and print the result in the subtitle
#' @param ... Other graphical parameters 
#' 
#' @examples
#' repart = c(20, 40, 40)
#' X = data.frame(num = c(rnorm(repart[1], 10, 1), rnorm(repart[2], 11, 1), rnorm(repart[3], 10, 1)),
#'                grp = rep(c("A", "B", "C"), repart))
#' 
#' 
#' BoxPlot(X$num, X$grp, data = X, ylab = "num", main = "boxplot with confidence intervals")
#' # Confidence interval in red with mean in blue.
#' 
#' @export
BoxPlot <- function(num, grp = NULL, data = NULL, AnoVa = TRUE, risk = 0.05, lang = c("en", "fr"), verbose = TRUE, ...) {
   
   if (!is.null(grp)) {
      formula = num ~ grp
      if (AnoVa) { # if we want anova: we compute it and stock the p-value in the subtitle
         resanov = aov(formula, data = data)
         pval = summary(resanov)[[1]][["Pr(>F)"]][[1]]
      } else { # suspicion of non-normalityand/or heteroscedasticity so kruskal-Wallis instead.
         grp = as.factor(grp) # characters do not work
         pval = kruskal.test(num, grp, data = data)$p.value
      }
      phrase = ""
      if (pval < risk) {
         if (lang[1] == "fr") {
            phrase = "Les groupes semblent distincts en moyenne."
         } else {
            phrase = "Groups seem to differ in mean"
         }
      } else {
         if (lang[1] == "fr") {
            phrase = "Les groupes ne semblent pas distincts en moyenne"
         } else {
            phrase = "groups do not seem to differ in mean"
         }
      }
      test = ifelse(AnoVa, "(AnoVa)", "(Kruskal-Wallis)")
      if (verbose) {sub = paste("p-value=", round(pval, 4), " ", phrase, test)}
      else {sub = ""}
   } else { # 1 seule classe : grp=NULL
      sub = ""
      formula = num
      pval = NA
   }
   
   res = boxplot(formula, sub = sub, ...)
   
   res$pval = pval
   
   
   y0 = c()
   y1 = c()
   y2 = c()
   if (is.null(grp)) {
      grp = 1
      groupes = 1
      x0 = 1
   } else {
      groupes = sort(unique(grp))
      
      x0 = seq(from = 1, length.out = length(groupes), by = 1)
      
   }
   for (k in groupes) { # k n'est pas forcement un indice, c'est le nom du groupe en cours
      Xloc = num[grp == k]
      y0 = c(y0, mean(Xloc) - qnorm(1 - risk / 2) * sd(Xloc) / sqrt(length(Xloc)))
      y1 = c(y1, mean(Xloc) + qnorm(1 - risk / 2) * sd(Xloc) / sqrt(length(Xloc)))
      y2 = c(y2, mean(Xloc))
   }
   arrows(x0 = x0, y0 = y0, x1 = x0, y1 = y2, angle = 90, code = 3, lwd = 2, col = "blue")
   arrows(x0 = x0, y0 = y0, x1 = x0, y1 = y1, angle = 90, code = 3, lwd = 2, col = "red")
   res$confint = rbind(y1, y2, y0)
   row.names(res$confint) <- c("upper", "mean", "lower")
   res$risk = risk
   
   return(res)
}