R/para_abundance_CI.R

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

#' Mean or median abundance estimation and confidence intervals
#'
#' This function calculates point estimates and confidence intervals (CIs) for parasite abundance, using either the mean or the median as a measure of central tendency. Confidence intervals are estimated via a non-parametric bootstrap approach based on resampling (permutations) of the observed data. Specifically, the function implements bias-corrected and accelerated (BCa) bootstrap intervals, which adjust for both bias and skewness in the bootstrap distribution. This approach does not assume a specific underlying distribution and is particularly robust for overdispersed and zero-inflated parasitological data.
#'
#' Parasite abundance is defined as the number of individuals of a given parasite taxon per host. For each taxon, abundance metrics are calculated based on the observed counts across hosts. The function reshapes the dataset into long format and computes abundance statistics for each parasite taxon and grouping combination (if specified).
#' The following are estimated:
#'
#' \itemize{
#'  \item \code{A} is the total parasite abundance
#'  \item \code{nH} is the number of hosts analyzed
#'  \item \code{nH_inf} is the number of infected hosts
#' }
#'
#' Depending on the argument \code{c_median}, the function calculates:
#'
#' \itemize{
#'  \item Mean abundance \code{MeanA}: average number of parasites per host
#'  \item Median abundance \code{MedA}: median number of parasites per host
#' }
#'
#' Confidence intervals are estimated using a non-parametric bootstrap approach. Specifically, bias-corrected and accelerated (BCa) bootstrap intervals are computed by resampling the observed abundance values with replacement a specified number of times \code{perm}. This method adjusts for both bias and skewness in the bootstrap distribution.
#' Statistical considerations: parasite abundance data are typically overdispersed and zero-inflated, making parametric assumptions inappropriate in many cases. The use of bootstrap methods allows robust estimation of confidence intervals without assuming normality. Mean abundance is sensitive to extreme values, whereas median abundance provides a more robust measure under highly skewed distributions. When sample size is small, bootstrap confidence intervals may be unstable or wide, and results should be interpreted with caution. The interpretation of results remains the responsibility of the user.
#' @usage
#' para_abundance_CI(dataset, c_median = TRUE,
#'  sp_cols, group_vars = NULL,  perm = 2000, decimal_places = 2,
#'  combine_ci = FALSE,  conf_level = 0.95, verbose = FALSE)
#'
#' @param dataset Data frame with parasitological data.
#' @param c_median Logical. If \code{TRUE}, the results will include the median as a central tendency of measure; if \code{FALSE}, the results will include the mean of the data.
#' @param sp_cols Vector with the names of the columns containing abundance of parasites (taxa) to calculate the parasitological descriptors.
#' @param group_vars Vector with the names of categorical variables used to define groups (e.g., "Sex", "Site"). Default = \code{NULL}.
#' @param perm Number of permutations to perform for confidence interval estimation. Default = \code{2000}.
#' @param decimal_places Number of decimal places to include in the calculation. Default = \code{2}.
#' @param combine_ci Logical. If \code{TRUE}, the interval is expressed as a single column (min - max). If \code{FALSE}, the interval is split into separate lower and upper limit columns.
#' @param conf_level Confidence level for the interval estimation (e.g., \code{0.95} for 95\% CI).
#' @param verbose A logical value indicating if progress messages should be given. Default = \code{FALSE}.
#'
#' @return A data frame containing abundance estimates and confidence intervals for each parasite taxon, either globally or by group. The following variables are returned:
#' \itemize{
#'   \item \code{nH}: Number of hosts analyzed
#'   \item \code{nH_inf}: Number of infected hosts
#'   \item \code{A}: Total parasite abundance
#'   \item \code{MeanA}:Mean parasite abundance
#'   \item \code{MedA}: Median parasite abundance
#'   \item \code{Lower_CI}: Lower bound of the bootstrap confidence interval
#'   \item \code{Upper_CI}: Upper bound of the bootstrap confidence interval
#'   \item \code{CI}: If \code{combine_ci = TRUE}, confidence interval expressed are store as a single column (Lower_CI - Upper_CI)
#'   \item \code{Observation}: Categorical description of the data context:
#'        \itemize{
#'          \item \code{"Not analyzed"}: No valid observations are available for the given combination (all values are missing or the combination is absent in the dataset); therefore, no estimates can be computed.
#'          \item \code{"One host analyzed"}: Only a single host analyzed is available for the given combination; thus, no population-level inference is possible and statistical summary measures are not estimated.
#'          \item \code{"No hosts infested"}: Hosts are present for the given combination, but none are infested (abundance = 0 for all observations); consequently, no statistical summary measures of abundance or intensity can be estimated.
#'          \item \code{"One host infested"}: Only a single infested host is recorded for the given combination; therefore, no sample-based estimation of intensity or related summary measures is possible.
#'          \item \code{"Multiple hosts infested"}: More than one infested host is recorded for the given combination, allowing the estimation of summary measures.
#' }
#' }
#'
#' @examples
#' #Calculate the CI for the median abundance
#' med_abun_CI <- para_abundance_CI(para_data$dataset,
#'                                 c_median = TRUE,
#'                                 sp_cols =  c("Sp1"),
#'                                 group_vars = c("Site"),
#'                                 decimal_places = 2,
#'                                 conf_level = 0.95,
#'                                 combine_ci = TRUE,
#'                                 verbose = TRUE)
#' med_abun_CI
#' #Calculate the CI for the mean abundance
#' mean_abun_CI <- para_abundance_CI(para_data$dataset,
#'                                  c_median = FALSE,
#'                                  sp_cols =  c("Sp1"),
#'                                  group_vars = c("Site"),
#'                                  decimal_places = 2,
#'                                  conf_level = 0.95,
#'                                  combine_ci = TRUE,
#'                                  verbose = TRUE)
#' mean_abun_CI
#'
#' @references
#' Bush, A.O., Lafferty, K.D., Lotz, J.M., Shostak, A.W. (1997). Parasitology meets ecology on its own terms:
#' Margolis revisited. \emph{Journal of Parasitology}, 83(4), 575–583.
#'
#' Reiczigel, J., Marozzi, M., Fabian, I., Rózsa, L. (2019). Biostatistics for parasitologists – a primer to
#' quantitative parasitology. \emph{Trends in Parasitology}, 35(4), 277–281.
#'
#' @author Juan Manuel Cabrera, Exequiel Furlan and Elisa Helman
#'
#' @export

para_abundance_CI <- function(dataset, c_median = TRUE, sp_cols, group_vars = NULL,
						   perm = 2000, decimal_places = 2, combine_ci = FALSE,
						   conf_level=0.95, verbose = FALSE) {

  Abund<-NA
  A<-NA
  CI<-NA
  MeanA<-NA
  MedA<-NA
  nH_inf<-NA
  nH<-NA
  Sp<-NA
  boot_result<-NA
  bca_conf_interval<-NA
  Lower_CI<-NA
  Upper_CI<-NA
  Observation<-NA
  M<-NA


  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.")
  }

  # ---------------------------
  # Función bootstrap
  # ---------------------------
  if (c_median) {
    ab_function <- function(dataset, indices) {
      sample_data <- dataset[indices]
      return(stats::median(sample_data, na.rm = TRUE))
    }
  } else {
    ab_function <- function(dataset, indices) {
      sample_data <- dataset[indices]
      return(mean(sample_data, na.rm = TRUE))
    }
  }
  # ---------------------------
  # Formato largo
  # ---------------------------
  data_long <- dataset %>%
    tidyr::pivot_longer(
      cols = dplyr::all_of(sp_cols),
      names_to = "Sp",
      values_to = "Abund"
    )
  if (verbose) message("Calculating abundance...")
  # ---------------------------
  # Cálculo + bootstrap
  # ---------------------------
  abundance_results <- data_long %>%
    dplyr::group_by(dplyr::across(dplyr::all_of(group_vars)), Sp) %>%
    dplyr::summarise(
      nH = sum(!is.na(Abund)),
      nH_inf = sum(Abund > 0, na.rm = TRUE),
      A = dplyr::if_else(nH > 1, sum(Abund, na.rm = TRUE), NA_real_),
      MeanA = dplyr::if_else(
        nH > 1,
        round(mean(Abund, na.rm = TRUE), decimal_places),
        NA_real_
      ),
      MedA = dplyr::if_else(
        nH > 1,
        round(stats::median(Abund, na.rm = TRUE), decimal_places),
        NA_real_
      ),
      boot_result = dplyr::if_else(MeanA || MedA && nH > 1,
                                   list(boot::boot(data = Abund, statistic = ab_function, R = perm)),
                                   list(boot::boot(data = 0, statistic = ab_function, R = perm))
      ),
      .groups = "drop"
    ) %>%
    dplyr::rowwise() %>%
    dplyr::mutate(
      bca_conf_interval = list(
        tryCatch(suppressMessages(
          suppressWarnings(boot::boot.ci(boot_result, type = "bca", conf_level = conf_level))),
          error = function(e) {return (NULL)})),
      Lower_CI = ifelse(
        !is.null(bca_conf_interval),
        round(bca_conf_interval$bca[4], decimal_places),
        NA_real_),
      Upper_CI = ifelse(
        !is.null(bca_conf_interval),
        round(bca_conf_interval$bca[5], decimal_places),
        NA_real_)
    ) %>% dplyr::ungroup()

  # ---------------------------
  # Observaciones
  # ---------------------------
  abundance_results <- abundance_results %>%
    dplyr::mutate(
      Observation = dplyr::case_when(
        nH == 0 ~ "Not analyzed",
        nH == 1 ~ "One host examined",
        nH_inf == 0 ~ "Not hosts infested",
        nH_inf == 1 ~ "One host infested",
        TRUE ~ "Multiple hosts infested"
      )
    )

  # ---------------------------
  # WARNINGS
  # ---------------------------

  #CI no estimables
  failed_ci <- any(is.na(abundance_results$Lower_CI) | is.na(abundance_results$Upper_CI))

  if (failed_ci) {
    warning("Confidence intervals could not be estimated due to insufficient variability in the data.")
  }

  #CI poco confiables
  unreliable_ci <- any(abundance_results$nH_inf > 0 & abundance_results$nH <= 5)

  if (unreliable_ci) {
    warning("Some confidence intervals may be unreliable due to low variability or small sample size.")
  }

  # ---------------------------
  # Combinar CI
  # ---------------------------
  if (combine_ci) {
    abundance_results <- abundance_results %>%
      dplyr::mutate(
        CI = ifelse(
          nH_inf > 0,
          paste(Lower_CI, Upper_CI, sep = " - "),
          NA_character_
        )
      ) %>%
      dplyr::select(
        dplyr::all_of(group_vars),
        Sp, nH, nH_inf, A,
        MeanA, MedA,
        CI, Observation
      )
  } else {
    abundance_results <- abundance_results %>%
      dplyr::select(
        dplyr::all_of(group_vars),
        Sp, nH, nH_inf, A,
        MeanA, MedA,
        Lower_CI, Upper_CI,
        Observation
      )
  }

  # ---------------------------
  # Selección final
  # ---------------------------
  if (c_median) {
    abundance_results <- abundance_results %>%
      dplyr::select(-MeanA)
  } else {
    abundance_results <- abundance_results %>%
      dplyr::select(-MedA)
  }
  if (verbose) message("Calculation completed")
  return(abundance_results)
}