R/para_explo_abund.R

In parasiteR: A Theorical-Practical Approach to Parasitological Data Analysis

#' Exploratory plots of parasite abundance distributions
#'
#' Generates exploratory visualizations of parasite abundance distributions across taxa and optional grouping variables. The function produces histograms combined with kernel density curves to facilitate the assessment of distributional patterns, including skewness, dispersion and zero inflation.
#'
#' The function reshapes the input dataset into a long format, where parasite taxa are treated as a single variable and their abundances as observations. For each parasite taxon and combination of grouping variables (if provided), the function generates:
#' \itemize{
#'   \item A \strong{histogram} representing the distribution of parasite abundance values
#'   \item A \strong{kernel density curve} (when sufficient data are available), providing a smoothed approximation of the underlying distribution.
#' }
#' Both elements are scaled to represent density, allowing direct comparison between distributions. These plots are intended for exploratory purposes and should not be used as formal inference tools. Faceting is applied to display each taxon and grouping combination in separate panels.
#' Special cases are handled as follows:
#' \itemize{
#'   \item When all abundance values for a given combination are zero, no histogram or density curve is drawn and a message is displayed indicating that the parasite was not recorded for that combination.
#'   \item When the number of observations is insufficient (less than 2), a message is displayed indicating that there is not enough data to compute a meaningful distribution.
#'   \item Density curves are only computed when there are more than two observations and non-zero variation.
#'}
#' All plots use independent scales (free scales) to better represent the variability within each facet.
#'
#' @usage
#' para_explo_abund(dataset, sp_cols, group_vars = NULL,
#'  bins = 30, n_col = NULL, verbose = FALSE)
#'
#' @param dataset Data frame containing parasite data.
#' @param sp_cols Vector with the names of the columns containing parasite abundance (taxa) to be plotted.
#' @param group_vars Vector with the names of categorical variables used to define groups (e.g., "Sex", "Site"). Default = \code{NULL}.
#' @param bins Integer specifying the number of bins used in the histogram. Higher values provide finer resolution but may introduce noise, while lower values produce smoother but less detailed distributions. Default = \code{30}.
#' @param n_col Integer specifying the number of columns in the faceted plot layout. If  \code{NULL}, the number of columns is determined automatically by ggplot2.Default = \code{NULL}.
#' @param verbose A logical value indicating if progress messages should be given. Default = \code{FALSE}.
#'
#' @return A ggplot2 object containing the generated faceted plots. This object can be further customized using standard ggplot2 functions.
#'
#' @examples
#'
#'#Species 1 and 2
#'
#'para_explo_abund (para_data$dataset,
#'                  sp_cols = c("Sp1", "Sp2"),
#'                  group_vars = c("Site", "Sp_host"),
#'                  bins = 30,
#'                  n_col = 4,
#'                  verbose = TRUE)
#'
#'#Species 3 and 4
#'
#'para_explo_abund (para_data$dataset,
#'                  sp_cols = c("Sp3", "Sp4"),
#'                  group_vars = c("Site", "Sp_host"),
#'                  bins = 30,
#'                  n_col = 4,
#'                  verbose = TRUE)
#'
#' @author Juan Manuel Cabrera, Exequiel Furlan and Elisa Helman
#'
#' @export


para_explo_abund <- function(dataset, sp_cols, group_vars = NULL, bins = 30, n_col = NULL, verbose = FALSE) {

  Abund<-NULL
  n_obs<-NULL
  all_zero<-NULL
  density<-NULL
  label<-NULL
  all_na<-NULL

  if (verbose) message("Checking function arguments...")

  # ---------------------------
  # Validaciones
  # ---------------------------
  if (is.null(sp_cols) || length(sp_cols) == 0) {
    stop("The species columns must be specified (sp_cols).")
  }

  if (!all(sp_cols %in% colnames(dataset))) {
    stop("Some of the specified species columns do not exist in the dataset.")
  }

  if (!is.null(group_vars) && !all(group_vars %in% colnames(dataset))) {
    stop("Some of the specified categorical variables do not exist in the dataset.")
  }

  # ---------------------------

  if (verbose) message("Drawing abundance exploratory plots...")
  data_long <- dataset %>%
    tidyr::pivot_longer(cols = dplyr::all_of(sp_cols),
                 names_to = "Sp",
                 values_to = "Abund") %>%
    dplyr::mutate(Abund = as.numeric(as.character(Abund)))

  # Completar combinaciones
  data_long <- data_long %>%
    tidyr::complete(!!!dplyr::syms(c(group_vars, "Sp")))

  # Resumen por grupo
  summary_df <- data_long %>%
    dplyr::group_by(dplyr::across(dplyr::all_of(c(group_vars, "Sp")))) %>%
    dplyr::summarise(
      n_obs = sum(!is.na(Abund)),
      all_na = n_obs == 0,
      all_zero = n_obs > 0 & all(Abund == 0, na.rm = TRUE),
      .groups = "drop"
    )

  # Unir info al dataset
  data_long <- dplyr::left_join(data_long, summary_df,
                         by = c(group_vars, "Sp"))

  # Base
  p <- ggplot2::ggplot() +
    ggplot2::theme_minimal()

  # Histogramas (solo si hay variación real)
  p <- p +
    ggplot2::geom_histogram(
      data = data_long %>% dplyr::filter(n_obs > 1 & !all_zero & !all_na),
      ggplot2::aes(x = Abund, y = ggplot2::after_stat(density)),
      bins = bins,
      fill = "grey70",
      color = "black"
    )

  # Densidad (solo si tiene sentido)
  p <- p +
    ggplot2::geom_density(
      data = data_long %>% dplyr::filter(n_obs > 2 & !all_zero & !all_na),
      ggplot2::aes(x = Abund),
      color = "red",
      linewidth = 1
    )

  # Mensajes para casos problemáticos
  message_df <- summary_df %>%
    dplyr::filter(all_na | all_zero | n_obs == 1) %>%
    dplyr::mutate(
      label = dplyr::case_when(
        all_na ~ "Not analyzed",
        n_obs == 1 ~ "One host\nanalyzed",
        all_zero ~ "Parasite\nabundance\nzero"
      )
    )

  p <- p +
    ggplot2::geom_text(
      data = message_df,
      ggplot2::aes(x = 0, y = 0, label = label),
      color = "black",
      size = 3,
      hjust = 0.5,
      vjust = 0.5
    )

  # Facetas
  if (!is.null(group_vars)) {
    p <- p + ggplot2::facet_wrap(
      stats::as.formula(paste("~", paste(c(group_vars, "Sp"), collapse = " + "))),
      scales = "free",
      ncol = n_col
    )
  } else {
    p <- p + ggplot2::facet_wrap(~Sp, scales = "free", ncol = n_col)
  }

  p <- p + ggplot2::labs(
    x = "Parasite abundance",
    y = "Density"
  )
  if (verbose) message("Calculation completed")
  return(p)
}