R/Graph_Boxplots.R

In OceanTracer: Marine Source Unmixing Toolkit

#' Box and Whiskers Plot for Tracers
#'
#' Visualizes the distribution of continuous tracer variables across sources and 
#' mixtures. It displays summary statistics (median, hinges, whiskers) alongside individual outliers.
#'
#' @param data Data frame containing source and mixture data in the standard format.
#' @param tracers Numeric vector. Indices of the tracers to plot (1 refers to the first tracer column). Default is 1 to (ncol(data) - 2).
#' @param ncol Integer. Number of plot columns per row in the grid layout. Default is 3.
#' @param colors Character vector. Optional hex codes or color names for custom boxplot colors. Default is NULL.
#' @param seed Integer. Random seed for reproducible jittering or sampling. Default is NULL.
#'
#' @return Invisibly returns the combined plot grid object created by \code{gridExtra::grid.arrange()}.
#'
#' @examples
#' data <- data.frame(
#'    ID = paste0("S", 1:12),
#'    Group = rep(c("SourceA", "SourceB", "Mixture"), each = 4),
#'    Tracer1 = c(10, 12, 11, 13, 20, 22, 21, 23, 15, 16, 15, 17),
#'    Tracer2 = c(5, 6, 5, 7, 15, 14, 16, 15, 10, 11, 10, 12)
#' )
#' Graph_boxplots(data, tracers = 1:2, ncol = 2)
#'
#' @export
Graph_boxplots <- function(data, tracers = 1:(ncol(data) - 2), ncol = 3, colors = NULL, seed = NULL) {
  
  # ── Internal helpers ────────────────────────────────────────────────────────
  
  is_averaged_internal <- function(d) {
    if (sum(grepl("^n$", tolower(colnames(d)))) != 1) return(FALSE)
    if ((ncol(d) - 3) %% 2 != 0) return(FALSE)
    ntracer <- (ncol(d) - 3) / 2
    nms <- colnames(d)[-c(1:2)][1:(2 * ntracer)]
    v1  <- nms[grepl("^[Mm]ean_", nms)]
    v2  <- nms[grepl("^[Ss][Dd]_", nms)]
    length(v1) == ntracer && length(v2) == ntracer &&
      all(gsub("^[Mm]ean_", "", v1) == gsub("^[Ss][Dd]_", "", v2))
  }
  
  raw_dataset_internal <- function(d) {
    ntracer <- (ncol(d) - 3) / 2
    n_col   <- which(tolower(colnames(d)) == "n")
    
    d[, 2] <- factor(d[, 2], levels = unique(d[, 2]))
    mixture_label <- levels(d[, 2])[nlevels(d[, 2])]
    is_mixture    <- d[, 2] == mixture_label
    
    raw_data <- data.frame()
    
    # Sources: expand using n, mean, sd
    source_rows <- d[!is_mixture, ]
    for (i in seq_len(nrow(source_rows))) {
      n_samples <- source_rows[i, n_col]
      if (n_samples < 1) n_samples <- 1
      temp_data <- data.frame(ID    = rep(source_rows[i, 1], n_samples),
                              Group = rep(as.character(source_rows[i, 2]), n_samples))
      for (j in seq_len(ntracer)) {
        mean_val <- source_rows[i, 2 + j]
        sd_val   <- source_rows[i, 2 + ntracer + j]
        temp_data[colnames(source_rows)[2 + j]] <-
          if (n_samples > 1) stats::rnorm(n_samples, mean = mean_val, sd = sd_val)
        else               mean_val
      }
      raw_data <- rbind(raw_data, temp_data)
    }
    
    # Mixture: append directly
    mix_rows <- d[is_mixture, ]
    for (i in seq_len(nrow(mix_rows))) {
      temp_data <- data.frame(ID    = mix_rows[i, 1],
                              Group = as.character(mix_rows[i, 2]))
      for (j in seq_len(ntracer)) {
        temp_data[colnames(mix_rows)[2 + j]] <- mix_rows[i, 2 + j]
      }
      raw_data <- rbind(raw_data, temp_data)
    }
    
    colnames(raw_data) <- gsub("^[Mm]ean_", "", colnames(raw_data))
    raw_data
  }
  
  if (!is.null(seed)) set.seed(seed)
  
  # ── Pre-processing ──────────────────────────────────────────────────────────
  if (is_averaged_internal(data)) {
    data <- raw_dataset_internal(data)
  }
  
  group_col <- colnames(data)[2]
  data[, 2] <- factor(data[, 2], levels = unique(data[, 2]))
  
  # Identify target tracer columns by index
  all_tracer_cols <- colnames(data)[-c(1, 2)]
  selected_tracer_names <- all_tracer_cols[tracers]
  
  # Helper functions for stat_summary
  make_fun1 <- function(sub_df) {
    function(x) {
      y <- sub_df$value
      data.frame(
        y = min(y, na.rm = TRUE) - 0.05 * (max(y, na.rm = TRUE) - min(y, na.rm = TRUE)),
        label = paste0("n=", length(x))
      )
    }
  }
  
  make_fun2 <- function(sub_df) {
    function(x) {
      y <- sub_df$value
      data.frame(
        y = stats::median(x, na.rm = TRUE) + 0.03 * (max(y, na.rm = TRUE) - min(y, na.rm = TRUE)),
        label = signif(mean(x, na.rm = TRUE), 3)
      )
    }
  }
  
  g_legend <- function(a.gplot) {
    tmp <- ggplot2::ggplot_gtable(ggplot2::ggplot_build(a.gplot))
    leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
    if (length(leg) > 0) tmp$grobs[[leg]] else NULL
  }
  
  blank <- grid::grid.rect(gp = grid::gpar(col = "white"))
  glist <- list()
  j <- 1
  
  for (i in seq_along(selected_tracer_names)) {
    cur_var <- selected_tracer_names[i]
    
    sub_df <- data.frame(
      Group = data[, 2],
      value = data[[cur_var]]
    )
    colnames(sub_df)[1] <- group_col
    
    p <- ggplot2::ggplot(
      sub_df,
      ggplot2::aes(x = .data[[group_col]], y = .data[["value"]], color = .data[[group_col]])
    ) +
      ggplot2::geom_boxplot() +
      ggplot2::ggtitle(as.character(cur_var)) +
      ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
      ggplot2::xlab("") + ggplot2::ylab("") +
      ggplot2::stat_summary(fun.data = make_fun1(sub_df), geom = "text",
                            fontface = "bold", color = "black", alpha = 0.5, size = 4) +
      ggplot2::stat_summary(fun.data = make_fun2(sub_df), geom = "text",
                            fontface = "bold", color = "black", alpha = 0.5, size = 4) +
      ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90),
                     legend.position = "none")
    
    if (!is.null(colors)) {
      # Handle color names safely to prevent scale mismatch warnings
      plot_colors <- colors
      if (!is.null(names(plot_colors)) && !all(names(plot_colors) %in% levels(sub_df[[group_col]]))) {
        plot_colors <- unname(plot_colors)
      }
      p <- p +
        ggplot2::scale_color_manual(values = plot_colors) +
        ggplot2::scale_fill_manual(values = plot_colors)
    }
    
    glist[[j]] <- p
    
    if (i == 1) {
      mylegend <- g_legend(p)
      if (!is.null(mylegend)) { j <- j + 1; glist[[j]] <- mylegend }
    } else if (i %% ncol == 0) {
      j <- j + 1; glist[[j]] <- blank
    }
    j <- j + 1
  }
  
  res_grid <- gridExtra::grid.arrange(
    grobs = glist,
    ncol = ncol + 1,
    widths = c(rep(10, each = ncol), 3)
  )
  
  invisible(res_grid)
}

# ── Aliases ──────────────────────────────────────────────────────────────────
#' @rdname Graph_boxplots
#' @export
boxPlot <- Graph_boxplots

#' @rdname Graph_boxplots
#' @export
box_plot <- Graph_boxplots