R/get_pvalue_position.R

In rstatix: Pipe-Friendly Framework for Basic Statistical Tests

#' @include utilities.R
NULL
#'Autocompute P-value Positions For Plotting Significance
#'@description Compute p-value x and y positions for plotting significance
#'  levels.
#'
#'  See the Datanovia tutorials
#'  \href{https://www.datanovia.com/learn/data-visualization/ggpubr/p-values-from-tests}{P-values from Tests on ggplots in R (rstatix)}
#'  and
#'  \href{https://www.datanovia.com/learn/data-visualization/ggpubr/auto-p-values}{Auto P-values in ggplot with geom_pwc (ggpubr)}
#'  for worked walkthroughs.
#'@inheritParams t_test
#'@param ref.group a character string specifying the reference group. If
#'  specified, for a given grouping variable, each of the group levels will be
#'  compared to the reference group (i.e. control group).
#'@param fun summary statistics functions used to compute automatically suitable
#'  y positions of p-value labels and brackets. Possible values include:
#'  \code{"max", "mean", "mean_sd", "mean_se", "mean_ci", "median",
#'  "median_iqr", "median_mad"}.
#'
#'  For example, if \code{fun = "max"}, the y positions are guessed as follow:
#'  \itemize{ \item 1. Compute the maximum of each group (groups.maximum) \item
#'  2. Use the highest groups maximum as the first bracket y position \item 3.
#'  Add successively a step increase for remaining bracket y positions. }
#'
#'  When the main plot is a boxplot, you need the option \code{fun = "max"}, to
#'  have the p-value bracket displayed at the maximum point of the group.
#'
#'  In some situations the main plot is a line plot or a barplot showing the
#'  \code{mean+/-error bars} of the groups, where error can be SE (standard
#'  error), SD (standard deviation) or CI (confidence interval). In this case,
#'  to correctly compute the bracket y position you need the option \code{fun =
#'  "mean_se"}, etc.
#'@param step.increase numeric vector with the increase in fraction of total
#'  height for every additional comparison to minimize overlap.
#'@param y.trans a function for transforming y axis scale. Value can be
#'  \code{log2}, \code{log10} and \code{sqrt}. Can be also any custom function
#'  that can take a numeric vector as input and returns a numeric vector,
#'  example: \code{y.trans = function(x){log2(x+1)}}
#'@param test an object of class \code{rstatix_test} as returned by
#'  \code{\link{t_test}()}, \code{\link{wilcox_test}()},
#'  \code{\link{sign_test}()}, \code{\link{tukey_hsd}()},
#'  \code{\link{dunn_test}()}.
#'@param x variable on x axis.
#'@param group group variable (legend variable).
#'@param dodge dodge width for grouped ggplot/test. Default is 0.8. Used only
#'  when \code{x} specified.
#'@param stack logical. If TRUE, computes y position for a stacked plot. Useful
#'  when dealing with stacked bar plots.
#'@param scales Should scales be fixed (\code{"fixed"}, the default), free
#'  (\code{"free"}), or free in one dimension (\code{"free_y"})?. For the y
#'  position, \code{"free"} or \code{"free_y"} compute the step increase per plot
#'  panel (otherwise a global step increase is used). For the x position, when
#'  \code{scales = "free"} the x positions (\code{xmin}/\code{xmax}) are computed
#'  \strong{per facet}, so that brackets align correctly when a faceted plot
#'  (\code{facet_*(scales = "free")}) shows a different set of x-axis levels in
#'  each panel (\code{"free_y"} keeps the global x positions).
#'@section Added columns: \code{add_y_position()} adds \code{y.position};
#'  \code{add_x_position()} adds \code{x}, \code{xmin}, \code{xmax}; and
#'  \code{add_xy_position()} adds all of these. A column named \code{groups} (the
#'  pair of compared groups, used internally to position dodged grouped
#'  comparisons) is also added. \strong{\code{groups} is a reserved internal
#'  column}: do not rely on it and avoid using \code{groups} as one of your own
#'  column names in the test object, since it is overwritten here. (It may be
#'  renamed to a dotted, less collision-prone name in a future major version.)
#' @examples
#' # Data preparation
#' #::::::::::::::::::::::::::::::::::::
#' df <- ToothGrowth
#' df$dose <- as.factor(df$dose)
#' df$group <- factor(rep(c(1, 2), 30))
#' head(df)
#'
#' # Stat tests
#' #::::::::::::::::::::::::::::::::::::
#' stat.test <- df %>%
#'   t_test(len ~ dose)
#' stat.test
#'
#' # Add the test into box plots
#' #::::::::::::::::::::::::::::::::::::
#' stat.test <- stat.test %>%
#'   add_y_position()
#' \donttest{
#'  if(require("ggpubr")){
#'    ggboxplot(df, x = "dose", y = "len") +
#'      stat_pvalue_manual(stat.test, label = "p.adj.signif", tip.length = 0.01)
#'   }
#' }
#' @seealso The Datanovia tutorials: \href{https://www.datanovia.com/learn/data-visualization/ggpubr/p-values-from-tests}{P-values from Tests on ggplots in R (rstatix)}, \href{https://www.datanovia.com/learn/data-visualization/ggpubr/auto-p-values}{Auto P-values in ggplot with geom_pwc (ggpubr)}.

#'@describeIn get_pvalue_position compute the p-value y positions
#'@export
get_y_position <- function(data, formula, fun = "max", ref.group = NULL, comparisons = NULL,
                           step.increase = 0.12, y.trans = NULL, stack = FALSE,
                           scales = c("fixed", "free", "free_y")){
  # Estimate step increase
  # 1. Get groups y scale
  outcome <- get_formula_left_hand_side(formula)
  group <- get_formula_right_hand_side(formula)
  if(.is_empty(group))
    group <- NULL
  if(is_grouped_df(data))
    group <- c(group , dplyr::group_vars(data))
  yscale <- get_y_scale(data, outcome, group, fun, stack = stack)
  # 2. Step increase
  # If fixed scales, then a global step increase is computed,
  # otherwise step.increase is estimated by panel in get_y_position_core().
  scales <- match.arg(scales)
  if(scales == "fixed"){
    step.increase <- step.increase*(yscale$max - yscale$min)
  }
  get_y_position_core(
    data = data, formula = formula, fun = fun, ref.group = ref.group,
    comparisons = comparisons, step.increase = step.increase, y.trans = y.trans,
    stack = stack, scales = scales
    )
}


get_y_position_core <- function(data, formula, fun = "max", ref.group = NULL, comparisons = NULL,
                            step.increase = 0.12, y.trans = NULL, stack = FALSE, scales = "fixed"){
  if(is_grouped_df(data)){
    results <- data %>%
      doo(
        get_y_position_core, formula = formula, fun = fun,
        ref.group = ref.group, comparisons = comparisons,
        step.increase = step.increase, y.trans = y.trans,
        stack = stack, scales = scales
        )
    return(results)
  }

  outcome <- get_formula_left_hand_side(formula)
  group <- get_formula_right_hand_side(formula)
  # Possible comparisons between groups
  if(is.null(comparisons)){
    if(.is_empty(group)){
      comparisons <- list(c("1", "null model"))
    }
    else{
      comparisons <- data %>%
        get_comparisons(!!group, !!ref.group)
    }
  }
  ncomparisons <- length(comparisons)
  group1 <- comparisons %>% get_group(1)
  group2 <- comparisons %>% get_group(2)

  # Estimate y axis scale
  yscale <- get_y_scale(data, y = outcome, group = group, fun = fun, stack = stack)
  if(is.null(step.increase)) step.increase <- yscale$max/20
  else if(scales %in% c("free", "free_y")){
    step.increase <- step.increase*(yscale$max - yscale$min)
  }
  ystart <- yscale$max + step.increase
  # The seq() below spreads ncomparisons brackets over [ystart, yend], so the
  # gap between consecutive brackets is (yend - ystart)/(ncomparisons - 1).
  # Using (ncomparisons - 1) here makes that gap exactly step.increase (#201).
  yend <- ystart + (step.increase*(ncomparisons-1))

  if(is.null(ref.group)) ref.group <- ""
  if(ref.group %in% c("all", ".all.")){
    y.position <- yscale$y + step.increase
  }
  else{
    y.position <- seq(
      from = ystart, to = yend,
      length.out = ncomparisons
    )
  }
  if(!is.null(y.trans)) y.position <- y.trans(y.position)
  results <- tibble(group1, group2, y.position) %>%
    mutate(groups = combine_this(group1, group2))
  results
}



#' @describeIn get_pvalue_position add p-value y positions to an object of class \code{rstatix_test}
#' @export
add_y_position <- function(test, fun = "max", step.increase = 0.12,
                           data = NULL, formula = NULL,  ref.group = NULL, comparisons = NULL,
                           y.trans = NULL, stack = FALSE, scales = c("fixed", "free", "free_y"))
  {
  scales <- match.arg(scales)
  asserttat_group_columns_exists(test)
  .attributes <- get_test_attributes(test)
  args  <- get_test_arguments(test)
  test <- keep_only_tbl_df_classes(test)
  if(!is.null(args)){
    if(missing(data)) data <- args$data
    if(missing(formula)) formula <- args$formula
    if(missing(ref.group)) ref.group <- args$ref.group
    if(missing(comparisons)) comparisons <- args$comparisons
  }
  if(is.null(data) | is.null(formula)){
    stop("data and formula arguments should be specified.")
  }
  # Repeated-measures pairwise tests (e.g. friedman_conover_test) carry a
  # 'outcome ~ within | subject' formula. For bracket positioning only the
  # 'outcome ~ within' part is relevant, so drop the '| subject' block to let the
  # grouping (x) variable be parsed correctly. This is a no-op for the usual
  # 'outcome ~ group' formulas (which contain no '|'), so existing behavior is
  # unchanged.
  formula <- drop_formula_block_term(formula)
  # When no explicit comparisons were used, derive them from the test's own
  # group1/group2 rows so that y positions are computed for exactly the
  # comparisons present (e.g. after filtering the test) instead of for the full
  # comparison set, which would otherwise produce uneven bracket spacing (#197).
  # Restricted to the standard ungrouped pairwise case; grouped tests,
  # ref.group = "all" and one-sample tests keep their existing behavior.
  if(is.null(comparisons) && !is_grouped_df(data) && nrow(test) > 0 &&
     all(c("group1", "group2") %in% colnames(test))){
    is.special <- any(test$group1 %in% c("all", ".all.")) ||
      any(test$group2 %in% c("null model"))
    if(!is.special){
      comparisons <- purrr::map2(
        as.character(test$group1), as.character(test$group2), c
      )
      # name like get_comparisons() (V1, V2, ...) so the derived `groups`
      # list-column is identical to the unfiltered/default path
      names(comparisons) <- paste0("V", seq_along(comparisons))
    }
  }
  positions <- get_y_position(
    data = data, formula = formula, fun = fun,
    ref.group = ref.group, comparisons = comparisons,
    step.increase = step.increase, y.trans = y.trans,
    stack = stack, scales = scales
  )
  if(nrow(test) == nrow(positions)){
    test$y.position <- positions$y.position
    test$groups <- positions$groups
  }
  else{
    # this occurs when tests are grouped by two variables (for ex),
    # but y positions are grouped by one variable
    # merging positions and test data frame
    if("y.position" %in% colnames(test)){
      test <- test %>% select(-any_of("y.position"))
    }
    if("groups" %in% colnames(test)){
      test <- test %>% select(-any_of("groups"))
    }
    common.columns <- intersect(colnames(test), colnames(positions))
    test <- test %>% dplyr::left_join(positions, by = common.columns)
  }

  test %>%
    set_test_attributes(.attributes)
}

# Compute y scale depending on fun
get_y_scale <- function(data, y, group, fun = "max", stack = FALSE){
  if(!.is_empty(group)){
    desc.stat <- data %>%
      group_by(!!!syms(group)) %>%
      get_summary_stats(!!y, type = fun)
  }
  else {
    desc.stat <- data %>%
      get_summary_stats(!!y, type = fun)
  }
  # Add error bars positions if any
  fun.splitted <- unlist(strsplit(fun, "_", fixed = TRUE))
  .center <- fun.splitted[1]
  .error <- ifelse(length(fun.splitted) == 2, fun.splitted[2], 0)
  .center <- desc.stat %>% pull(tidyselect::all_of(.center))
  if(.error != 0) .error <- desc.stat %>% pull(tidyselect::all_of(.error))
  if(stack){
    .center <- rep(sum(.center), length(.center))
  }
  y <- .center + .error
  ymax <- max(y, na.rm = TRUE)
  ymin <- min(y, na.rm = TRUE)
  list(y = y, max = ymax, min = ymin)
}


# Return group1 and group2 values from possible pairs
get_group <- function(possible.pairs, index){
  possible.pairs %>%
    map( function(x){as.character(x[index])}) %>%
    unlist()
}
# Combine vectors: c(a, b) & c(c, d) --> list(c(a, c), c(b, d))
# ... two or more character vectors
combine_this <- function(...){
  params <- list(...)
  purrr::pmap(params, c)
}


#' @describeIn get_pvalue_position compute and add p-value x positions.
#' @export
add_x_position <- function(test, x = NULL, group = NULL, dodge = 0.8,
                           scales = c("fixed", "free", "free_y")){
  scales <- match.arg(scales)
  # Checking
  asserttat_group_columns_exists(test)
  .attributes <- get_test_attributes(test)
  if(any(test$group1 %in% c("all", ".all."))) {
    # case when ref.group = "all"
    test$group1 <- test$group2
  }
  groups <- c(as.character(test$group1), as.character(test$group2)) %>%
    unique() %>%
    setdiff(c("all", ".all."))   # case when ref.group = "all"

  is_rstatix_test <- inherits(test, "rstatix_test")
  is.null.model <- ("null model" %in% test$group2) & all(test$group1 %in% 1)
  is.grouped.by.legend <- test %>% is_stat_test_grouped_by(group)
  is.grouped.by.x <- test %>% is_stat_test_grouped_by(x)
  is.basic <- !is.grouped.by.x & !is.grouped.by.legend

  # Data preparation
  if(is_rstatix_test) {
    data <- attr(test, "args")$data
    if(is.basic & is.null(x))
      x <- get_formula_right_hand_side(.attributes$args$formula)
    else if(is.grouped.by.x & is.null(group))
      group <- get_formula_right_hand_side(.attributes$args$formula)
  }
  else if(is.basic){
    data <- data.frame(x = groups, stringsAsFactors = FALSE)
    x <- "x"
  }
  else{
    if(is.grouped.by.x)
      data <- expand.grid(x = unique(test[[x]]), group = groups)
    else if(is.grouped.by.legend)
      data <- expand.grid(x = groups, group = test[[group]])
    colnames(data) <- c(x, group)
  }
  if(is.null.model) {
    data$group <- rep(1, nrow(data))
    group <- "group"
  }

  # Add xmin and x max
  if(is.basic){
    x_coords <- as_numeric_group(data[[x]])
    xmin_id <- as.character(test$group1)
    xmax_id <- as.character(test$group2)
  }
  else{
    x_coords <- get_grouped_x_position(data, x = x, group = group, dodge = dodge)
    if(is.grouped.by.legend){
      # Add x position to stat test when the test is grouped by the legend variable
      # Case when you group by legend and pairwise compare between x-axis groups
      xmin_id <- paste(test$group1, test[[group]], sep = "_")
      xmax_id <- paste(test$group2, test[[group]], sep = "_")
    }
    else if(is.grouped.by.x){
      # Add x position to stat test when the test is grouped by the x variable
      # Case when you pairwise compare legend groups at each x-axis position,
      # so the data is grouped by x position
      xmin_id <- paste(test[[x]], test$group1, sep = "_")
      xmax_id <- paste(test[[x]], test$group2, sep = "_")
      test$x <- unname(as_numeric_group(test[[x]]))
    }
  }
  test$xmin <- unname(x_coords[xmin_id])
  test$xmax <- unname(x_coords[xmax_id])
  if(is.null.model) test$xmax <- test$xmin
  # Free x scales in facets (#203): for a faceted test (grouped by one or more
  # variables) plotted with facet_*(scales = "free"), each facet shows only its
  # own x levels at positions 1..m, but x_coords above are global. Remap
  # xmin/xmax to each facet's positions so the brackets line up within the panel.
  # The panel's x levels are taken from the facet's actual DATA (in factor order),
  # not from the comparison endpoints, so a level that is plotted but not part of
  # any (retained) comparison still occupies its position. Only the basic
  # (non-dodged) case is affected; scales = "fixed" (default) / "free_y" leave the
  # global positions unchanged.
  facet.vars <- if(is.basic && scales == "free") .get_facet_vars(test) else character(0)
  facet.data <- attr(test, "args")$data
  if(length(facet.vars) > 0 && !is.null(facet.data) && !is.null(x) &&
     all(facet.vars %in% colnames(facet.data)) && x %in% colnames(facet.data)){
    x.levels <- levels(as.factor(facet.data[[x]]))   # global factor order
    test.fid <- do.call(paste, c(as.list(test[facet.vars]), sep = "\r"))
    data.fid <- do.call(paste, c(as.list(facet.data[facet.vars]), sep = "\r"))
    for(f in unique(test.fid)){
      present <- x.levels[x.levels %in% as.character(facet.data[[x]][data.fid == f])]
      rows <- test.fid == f
      test$xmin[rows] <- match(as.character(test$group1[rows]), present)
      test$xmax[rows] <- match(as.character(test$group2[rows]), present)
    }
  }
  test %>% set_test_attributes(.attributes)
}

# Leading grouping/facet variables of a (grouped) test result: the columns that
# appear before the first comparison/outcome column (.y., term or group1).
.get_facet_vars <- function(test){
  stop.cols <- c(".y.", "term", "group1")
  idx <- match(stop.cols, colnames(test))
  idx <- idx[!is.na(idx)]
  if(length(idx) == 0) return(character(0))
  first <- min(idx)
  if(first > 1) colnames(test)[seq_len(first - 1)] else character(0)
}


# Compute grouped x positions or coordinates
# data is a dataframe containing the x and the group columns
get_grouped_x_position<- function(data, x, group, dodge = 0.8){
  data <- data.frame(x = data[[x]], group = data[[group]]) %>%
    dplyr::distinct(.data$x, .data$group) %>%
    dplyr::arrange(.data$x, .data$group)
  data$x.position <- as_numeric_group(data$x)
  # Add group.ranks and ngroups at x position
  data <- data %>%
    rstatix::df_nest_by(vars = "x") %>%
    mutate(
      data = map(.data$data, function(data){data$group.ranks = 1:nrow(data); data}),
      n = unlist(map(.data$data, nrow))
    ) %>%
    tidyr::unnest(cols = "data")
  # Compute x coords
  d <- data
  x_coords <-  (((dodge - dodge*d$n) / (2*d$n)) + ((d$group.ranks - 1) * (dodge / d$n))) + d$x.position
  names(x_coords) <- paste(d$x, d$group, sep = "_")
  x_coords
}

# Check if a stat test is grouped by a given variable
is_stat_test_grouped_by <- function(test, x = NULL){
  answer <- FALSE
  if(!is.null(x)){
    if(x %in% colnames(test)){
      answer <- TRUE
    }
  }
  answer
}

# Return a numeric named vector
# c("a", "b", "a") -> c(a = 1, b = 2, a = 1)
as_numeric_group <- function(x){
  grp <- x %>% as.factor() %>% as.numeric()
  names(grp) <- x
  grp
}




#' @describeIn get_pvalue_position compute and add both x and y positions.
#' @export
add_xy_position <- function(test, x = NULL,  group = NULL, dodge = 0.8, stack = FALSE,
                            fun = "max", step.increase = 0.12,
                            scales = c("fixed", "free", "free_y"), ...){
  if(missing(dodge)){
    if(stack) dodge <- 0
  }
  test %>%
    add_y_position(
      fun = fun, step.increase = step.increase,
      stack = stack, scales = scales, ...
      ) %>%
    add_x_position(x = x, group = group, dodge = dodge, scales = scales)
}


# Helper functions
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
asserttat_group_columns_exists <- function(data){
  groups.exist <- all(c("group1", "group2") %in% colnames(data))
  if(!groups.exist){
    stop("data should contain group1 and group2 columns")
  }
}

.is_grouped_test <- function(test, x = NULL){
  answer <- FALSE
  if(!is.null(x)){
    if(x %in% colnames(test)){
      answer <- TRUE
    }
  }
  answer
}

.contains_selected_comparisons <- function(test){
  answer <- FALSE
  if(is_rstatix_test(test)){
    comparisons <- attr(test, "args")$comparisons
    answer <- !is.null(comparisons)
  }
  answer
}


# To be removed
add_x_position0 <- function(test, x = NULL, dodge = 0.8){

  asserttat_group_columns_exists(test)
  .attributes <- get_test_attributes(test)
  group1 <- set_diff(test$group1, c("all", ".all."), keep.dup = TRUE)
  group2 <- set_diff(test$group2, c("all", ".all."), keep.dup = TRUE)
  groups <- unique(c(group1, group2))
  if(.contains_selected_comparisons(test)){
    # get all data groups
    group.var <- get_formula_right_hand_side(.attributes$args$formula)
    groups <- .attributes$args$data %>%
      convert_as_factor(vars = group.var) %>%
      pull(tidyselect::all_of(group.var)) %>%
      levels()
  }
  group1.coords <- group1.ranks <-  match(group1, groups)
  group2.coords <- group2.ranks <- match(group2, groups)
  n <- length(groups)

  if(.is_grouped_test(test, x)){
    test <- test %>%
      .as_factor(x) %>%
      mutate(x = as.numeric(!!sym(x))) # x levels order
    xpos <- test$x
    group1.coords <- (((dodge - dodge*n) / (2*n)) + ((group1.ranks - 1) * (dodge / n))) + xpos
    group2.coords <- (((dodge - dodge*n) / (2*n)) + ((group2.ranks - 1) * (dodge / n))) + xpos
  }

  if(.is_empty(group1)) {
    # case when ref.group = "all"
    group1.coords <- group2.coords
  }
  test %>%
    mutate(xmin = group1.coords, xmax = group2.coords) %>%
    set_test_attributes(.attributes)
}