R/para_aggregation_CI.R

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

#' Aggregation indices and confidence intervals for parasite distributions
#'
#' This function calculates point estimates and bootstrap confidence intervals (CIs) for six commonly used measures of parasite aggregation among hosts. Following the recommendations of Morrill et al. (2023) and Reiczigel et al. (2019), CIs are estimated via non-parametric bootstrap using the bias-corrected and accelerated (BCa) method, except for the parameter k of the negative binomial distribution, for which maximum likelihood estimation with likelihood profile confidence intervals is used. Parasite aggregation refers to the nearly universal phenomenon where most hosts harbor few or no parasites, while a minority of hosts harbor heavy infections. This pattern has important consequences for host health, parasite transmission, and population dynamics. The six aggregation measures included in this function capture different aspects of this phenomenon and can be grouped into three distinct categories based on their properties and interpretations.
#'
#' Parasite aggregation is quantified through six indices that capture different attributes of how parasites are distributed among hosts. Following Morrill et al. (2023), these measures can be grouped into three categories based on their mathematical properties and biological interpretations:
#'
#' \itemize{
#'  \item \strong{Abundance-driven measures (sensitive to mean abundance):}
#'    \itemize{
#'            \item VMR (Variance-to-Mean Ratio): Calculated as the variance divided by the mean abundance. Values > 1 indicate aggregation, values = 1 indicate a random (Poisson) distribution, and values < 1 indicate a more even distribution. This measure is strongly correlated with mean abundance, making it unsuitable for comparing aggregation across samples with different mean abundances.
#'            \item Mean crowding: Represents the average number of conspecific parasites sharing a host from the perspective of an individual parasite. Higher values indicate greater crowding experienced by parasites within hosts. This measure is strongly positively correlated with mean abundance, making it unsuitable for comparing aggregation across samples with different mean abundances.
#'            }
#'  \item \strong{Distribution-driven measures (independent of mean abundance):}
#'    \itemize{
#'            \item \emph{k} of the negative binomial distribution: An inverse measure of aggregation where smaller values indicate higher aggregation. When k tends towards infinity, the distribution approaches a Poisson (random) pattern. This measure has a clear biological interpretation in terms of mean crowding: 1/k represents the proportion of the mean abundance by which mean crowding exceeds the mean abundance. Importantly, k is not necessarily correlated with mean abundance, contrary to some expectations in the literature.
#'            \item Patchiness: Calculated as mean crowding divided by the mean abundance. This measure describes how many times more "crowded" an average parasite is compared to a hypothetical random distribution. Like k, patchiness is not correlated with mean abundance and does not require that the parasite distribution fits a negative binomial.
#'            }
#'  \item \strong{Inequality-driven measures (based on Lorenz curves):}
#'    \itemize{
#'            \item Poulin's D: Ranges from 0 (perfectly even distribution) to 1 (maximum possible aggregation, all parasites in one host). This measure is based on the Lorenz curve and is particularly suitable for comparing aggregation across samples or different host-parasite associations. It shows a strong negative correlation with prevalence.
#'            \item Hoover's index: Also ranges from 0 to 1 and represents the proportion of parasites that would need to be redistributed to achieve a perfectly even distribution among hosts. For example, a value of 0.7 indicates that 70\% of the parasites would need to be moved to achieve evenness. This measure has a clear biological interpretation but can be constrained at low prevalence values (often equaling 1 - prevalence when all infected hosts have infrapopulations at least as large as the mean abundance).
#'            }
#' }
#'
#' Confidence intervals are estimated using a non-parametric bootstrap approach. Specifically, BCa bootstrap intervals are computed by resampling the observed abundance values with replacement a specified number of times (perm). This method adjusts for both bias and skewness in the bootstrap distribution and is recommended for aggregation measures by Morrill et al. (2023) and Reiczigel et al. (2019). For the parameter k of the negative binomial distribution, confidence intervals are estimated using maximum likelihood estimation with likelihood profile confidence intervals, following the specific recommendation of Morrill et al. (2023). For all other indices (VMR, mean crowding, patchiness, Poulin's D, and Hoover's index), BCa bootstrap intervals are used.
#' 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. When sample size is small (fewer than 3 hosts), CIs may not be estimable and results should be interpreted with caution. The interpretation of results remains the responsibility of the user.
#'
#' @usage
#' para_aggregation_CI(dataset, sp_cols, group_vars = NULL, perm = 2000,
#'  decimal_places = 2, ci_method= "bca", combine_ci = FALSE,
#'  conf_level = 0.95, save_csv = FALSE, verbose = FALSE)
#'
#'
#'
#' @param dataset Data frame with parasitological 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 is \code{NULL}.
#' @param perm Number of permutations to perform for CIs estimation. Default is \code{2000}.
#' @param decimal_places Number of decimal places to include in the calculation. Default is \code{2}.
#' @param ci_method Method for bootstrap CIs. Options are \code{"bca"} (bias-corrected and accelerated) or \code{"perc"} (percentile). Default = \code{"bca"}.
#' @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\% CIs).
#' @param save_csv Logical. If \code{TRUE}, the resulting table is automatically exported as a .CSV file to the current working directory. Default = \code{FALSE}.
#' @param verbose Logical. If \code{TRUE}, progress messages should be given. Default = \code{FALSE}.
#'
#' @return A data frame containing aggregation estimates and CIs 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{VMR}: Variance-to-Mean Ratio
#'   \item \code{mean_crowding}: Mean crowding index
#'   \item \code{patchiness}: Patchiness index
#'   \item \code{poulin_D}: Poulin's D index
#'   \item \code{hoover}: Hoover's index
#'   \item \code{k}: Parameter of the negative binomial distribution
#'   \item \code{Lower_<ci_method>}: Lower bound of the CI estimated using the method specified in \code{ci_method} (\code{"bca"} or \code{"perc"}).
#'   \item \code{Upper_<ci_method>}: Upper bound of the CI estimated using the method specified in \code{ci_method} (\code{"bca"} or \code{"perc"}).
#'   \item \code{CI_<ci_method>}: If \code{combine_ci = TRUE}, confidence interval stored as a single column (\code{Lower_<ci_method> - Upper_<ci_method>}).
#'   \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 aggregation indices with BCa confidence intervals
#'agg_CI <- para_aggregation_CI(para_data$dataset,
#'                                sp_cols =  c("Sp1","Sp2"),
#'                                group_vars = c("Site","Sp_host"),
#'                                perm = 2000,
#'                                decimal_places = 2,
#'                                ci_method = "bca",
#'                                combine_ci = FALSE,
#'                                conf_level = 0.95,
#'                                save_csv = FALSE,
#'                                verbose = TRUE)
#'agg_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.
#'
#' Morrill, A., Poulin, R., Forbes, M., 2023. Interrelationships and properties of parasite aggregation measures:
#' a user’s guide. Int. \emph{Journal of Parasitology} 53, 763–776.
#'
#' @author Juan Manuel Cabrera, Exequiel Furlan and Elisa Helman
#'
#' @export

para_aggregation_CI <- function(dataset, sp_cols, group_vars = NULL,
                              perm = 2000, decimal_places = 2, ci_method = "bca",
                              combine_ci = FALSE, conf_level = 0.95, save_csv = FALSE,
                              verbose = FALSE){

  Abund<- NA
  CI <- NA
  MeanI <- NA
  MedI <- NA
  nH_inf <- NA
  nH <- NA
  VMR <- NA
  mean_crowding <- NA
  patchiness <- NA
  poulin_D <- NA
  hoover <- NA
  k <- NA
  VMR_CI_lower <- NA
  VMR_CI_upper <- NA
  MC_CI_lower <- NA
  MC_CI_upper <- NA
  patch_CI_lower <- NA
  patch_CI_upper <- NA
  D_CI_lower <- NA
  D_CI_upper <- NA
  H_CI_lower <- NA
  H_CI_upper <- NA
  k_CI_lower <- NA
  k_CI_upper <- NA
  Observation <- NA
  Sp <- NA
  VMR_CI <- NA
  MC_CI <- NA
  patch_CI <- NA
  D_CI <- NA
  H_CI <- NA
  k_CI <- NA
  k_CI_ML <- NA
  k_Lower_ML <- NA
  k_Upper_ML <- NA
  current_date<-NA
  file_name<-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.")
  }
  if (!(ci_method %in% c("bca", "perc"))) {
    stop("Invalid ci_method select. Choose between \"bca\" or \"perc\"")
  }

  calc_aggregation <- function(x) {
    x <- x[!is.na(x)]

    if (length(x) < 2) {
      return(list(
        VMR = NA, mean_crowding = NA, patchiness = NA,
        poulin_D = NA, hoover = NA, k = NA
      ))
    }

    mean_x <- mean(x)
    var_x <- stats::var(x)

    # VMR
    vmr <- var_x / mean_x

    # Mean crowding
    mc <- (sum(x^2) / sum(x)) - 1

    # Patchiness
    patch <- mc / mean_x

    # Poulin D
    x_sorted <- sort(x)
    N <- length(x)
    cum <- cumsum(x_sorted)
    D <- 1 - (2 * sum(cum)) / (mean_x * N * (N + 1))

    # Hoover
    H <- sum(abs(x / mean_x - 1)) / (2 * N)

    # k (MLE)
    k <- tryCatch({
      suppressWarnings(MASS::fitdistr(x, "negative binomial")$estimate[1])
    }, error = function(e) NA_real_)

    return(list(
      VMR = vmr,
      mean_crowding = mc,
      patchiness = patch,
      poulin_D = D,
      hoover = H,
      k = k
    ))
  }

  # Funciones para bootstrap (BCa para todos excepto k)
  vmr_func <- function(z) {
    z <- z[!is.na(z)]
    if (length(z) < 2) return(NA)
    stats::var(z)/mean(z)
  }

  mc_func <- function(z) {
    z <- z[!is.na(z)]
    if (length(z) < 2) return(NA)
    (sum(z^2)/sum(z)) - 1
  }

  patch_func <- function(z) {
    z <- z[!is.na(z)]
    if (length(z) < 2) return(NA)
    mc <- (sum(z^2)/sum(z)) - 1
    mc / mean(z)
  }

  d_func <- function(z) {
    z <- z[!is.na(z)]
    if (length(z) < 2) return(NA)
    z_sorted <- sort(z)
    N <- length(z)
    mean_z <- mean(z)
    cum <- cumsum(z_sorted)
    1 - (2 * sum(cum)) / (mean_z * N * (N + 1))
  }

  hoover_func <- function(z) {
    z <- z[!is.na(z)]
    if (length(z) < 2) return(NA)
    mean_z <- mean(z)
    N <- length(z)
    sum(abs(z / mean_z - 1)) / (2 * N)
  }


  get_k_profile_ci <- function(x, conf = conf_level) {
    x <- x[!is.na(x)]

    if (length(x) < 3) {
      return(list(k = NA, lower = NA, upper = NA))
    }

    # Estimar k por MLE
    fit <- tryCatch({
      suppressWarnings(MASS::fitdistr(x, "negative binomial"))
    }, error = function(e) NULL)

    if (is.null(fit)) {
      return(list(k = NA, lower = NA, upper = NA))
    }

    k_est <- fit$estimate[1]
    mu_est <- fit$estimate[2]

    # Intentar perfil de verosimilitud con bbmle
    if (requireNamespace("bbmle", quietly = TRUE)) {
      tryCatch({
        # Función de verosimilitud negativa
        nll <- function(k) {
          -sum(stats::dnbinom(x, size = k, mu = mu_est, log = TRUE))
        }

        # Perfil de verosimilitud
        mle_fit <- bbmle::mle2(nll, start = list(k = k_est))
        prof <- bbmle::profile(mle_fit)
        ci <- bbmle::confint(prof, level = conf)

        return(list(
          k = k_est,
          lower = as.numeric(ci[1]),
          upper = as.numeric(ci[2])
        ))
      }, error = function(e) {
        # Si falla, usar bootstrap como alternativa
        return(get_k_bootstrap_ci(x, conf))
      })
    } else {
      # Si bbmle no está instalado, usar bootstrap
      if (verbose) message("  'bbmle' not installed. Using bootstrap for k CI.")
      return(get_k_bootstrap_ci(x, conf))
    }
  }

  # Función alternativa: bootstrap para k (cuando el perfil falla)
  get_k_bootstrap_ci <- function(x, conf = conf_level) {
    x <- x[!is.na(x)]

    if (length(x) < 3) {
      return(list(k = NA, lower = NA, upper = NA))
    }

    # Estimar k por MLE
    fit <- tryCatch({
      suppressWarnings(MASS::fitdistr(x, "negative binomial"))
    }, error = function(e) NULL)

    if (is.null(fit)) {
      return(list(k = NA, lower = NA, upper = NA))
    }

    k_est <- fit$estimate[1]

    # Bootstrap para k
    k_func <- function(z) {
      z <- z[!is.na(z)]
      if (length(z) < 2) return(NA)
      tryCatch({
        suppressWarnings(MASS::fitdistr(z, "negative binomial")$estimate[1])
      }, error = function(e) NA_real_)
    }

    tryCatch({
      boot_k <- boot::boot(x,
                           statistic = function(d, i) k_func(d[i]),
                           R = min(200, 200))  # Menos iteraciones por ser lento

      ci_k <- boot::boot.ci(boot_k, type = "bca", conf = conf)

      if (!is.null(ci_k)) {
        return(list(
          k = k_est,
          lower = ci_k$bca[4],
          upper = ci_k$bca[5]
        ))
      }
    }, error = function(e) {})

    return(list(k = k_est, lower = NA, upper = NA))
  }


  if (verbose) message("Calculating aggregation indices with confidence intervals...")

  # Formato largo
  data_long <- dataset %>%
    tidyr::pivot_longer(cols = dplyr::all_of(sp_cols),
                        names_to = "Sp",
                        values_to = "Abund")

  group_cols <- c(group_vars, "Sp")

  # Cálculo
  agg_results <- data_long %>%
    dplyr::group_by(dplyr::across(dplyr::all_of(group_cols))) %>%
    dplyr::summarise(
      nH = sum(!is.na(Abund)),
      nH_inf = sum(Abund > 0, na.rm = TRUE),
      .groups = 'drop'
    ) %>%
    dplyr::filter(nH >= 2)

    if (nrow(agg_results) == 0) {
    stop("No groups with sufficient data (minimum 2 hosts per group).")
  }

  # Inicializar lista de resultados
  results_list <- list()

  # Loop por cada grupo
  for (i in 1:nrow(agg_results)) {
    # Obtener el grupo actual
    current_group <- agg_results[i, ]

    # Filtrar datos del grupo
    group_data <- data_long |>
      dplyr::inner_join(current_group, by = dplyr::all_of(group_cols))

    # Extraer datos de abundancia
    x <- group_data$Abund
    x <- x[!is.na(x)]

    # Calcular índices puntuales
    indices <- calc_aggregation(x)

    # Inicializar CIs como NA
    vmr_ci_lower <- NA
    vmr_ci_upper <- NA
    mc_ci_lower <- NA
    mc_ci_upper <- NA
    patch_ci_lower <- NA
    patch_ci_upper <- NA
    d_ci_lower <- NA
    d_ci_upper <- NA
    h_ci_lower <- NA
    h_ci_upper <- NA
    k_ci_lower <- NA
    k_ci_upper <- NA

    # Calcular bootstrap BCa solo si hay suficientes datos (>= 3 hosts)
    # para VMR, Mean crowding, Patchiness, Poulin D, Hoover
    if (length(x) >= 3) {

      # VMR (BCa)
      tryCatch({
        boot_vmr <- boot::boot(x, statistic = function(d, i) vmr_func(d[i]), R = perm)
        ci_vmr <- boot::boot.ci(boot_vmr, type = ci_method, conf = conf_level)
        if (!is.null(ci_vmr)) {
          if (ci_method == "bca") {
            vmr_ci_lower <- ci_vmr$bca[4]
            vmr_ci_upper <- ci_vmr$bca[5]
          } else {
            if (!is.null(ci_vmr$percent)){
            vmr_ci_lower <- ci_vmr$percent[4]
            vmr_ci_upper <- ci_vmr$percent[5]
            } else {
              vmr_ci_lower <- NA_real_
              vmr_ci_upper <- NA_real_
            }
          }
        }
      }, error = function(e) {})

      # Mean crowding (BCa)
      tryCatch({
        boot_mc <- boot::boot(x, statistic = function(d, i) mc_func(d[i]), R = perm)
        ci_mc <- boot::boot.ci(boot_mc, type = ci_method, conf = conf_level)
        if (!is.null(ci_mc)) {
          if (ci_method == "bca") {
            mc_ci_lower <- ci_mc$bca[4]
            mc_ci_upper <- ci_mc$bca[5]
          } else {
            if (!is.null(ci_mc$percent)){
            mc_ci_lower <- ci_mc$percent[4]
            mc_ci_upper <- ci_mc$percent[5]
            } else {
              mc_ci_lower <- NA_real_
              mc_ci_upper <- NA_real_
          }
          }
        }
      }, error = function(e) {})

      # Patchiness (BCa)
      tryCatch({
        boot_patch <- boot::boot(x, statistic = function(d, i) patch_func(d[i]), R = perm)
        ci_patch <- boot::boot.ci(boot_patch, type = ci_method, conf = conf_level)
        if (!is.null(ci_patch)) {
          if (ci_method == "bca") {
            patch_ci_lower <- ci_patch$bca[4]
            patch_ci_upper <- ci_patch$bca[5]
          } else {
            if (!is.null(ci_patch$percent)){
            patch_ci_lower <- ci_patch$percent[4]
            patch_ci_upper <- ci_patch$percent[5]
            } else{
              patch_ci_lower <- NA_real_
              patch_ci_upper <- NA_real_
            }
          }
        }
      }, error = function(e) {})

      # Poulin D (BCa)
      tryCatch({
        boot_D <- boot::boot(x, statistic = function(d, i) d_func(d[i]), R = perm)
        ci_D <- boot::boot.ci(boot_D, type = ci_method, conf = conf_level)
        if (!is.null(ci_D)) {
          if (ci_method == "bca") {
            d_ci_lower <- ci_D$bca[4]
            d_ci_upper <- ci_D$bca[5]
          } else {
            if (!is.null(ci_D$percent)){
            d_ci_lower <-  ci_D$percent[4]
            d_ci_upper <- ci_D$percent[5]
            }else{
              d_ci_lower <- NA_real_
              d_ci_upper <- NA_real_
            }
          }
        }
      }, error = function(e) {})

      # Hoover (BCa)
      tryCatch({
        boot_H <- boot::boot(x, statistic = function(d, i) hoover_func(d[i]), R = perm)
        ci_H <- boot::boot.ci(boot_H, type = ci_method, conf = conf_level)
        if (!is.null(ci_H)) {
          if (ci_method == "bca") {
            h_ci_lower <- ci_H$bca[4]
            h_ci_upper <- ci_H$bca[5]
          } else {
            if (!is.null(ci_H$percent)){
            h_ci_lower <- ci_H$percent[4]
            h_ci_upper <- ci_H$percent[5]
            }else{
              h_ci_lower <- NA_real_
              h_ci_upper <- NA_real_
            }
          }
        }
      }, error = function(e) {})

      # k: PERFIL DE VEROSIMILITUD (RECOMENDACIÓN DE MORRILL)
      k_result <- get_k_profile_ci(x, conf_level)
      k_ci_lower <- k_result$lower
      k_ci_upper <- k_result$upper
    }

    # Crear fila de resultados
    result_row <- current_group |>
      dplyr::mutate(
        VMR = indices$VMR,
        mean_crowding = indices$mean_crowding,
        patchiness = indices$patchiness,
        poulin_D = indices$poulin_D,
        hoover = indices$hoover,
        k = indices$k,
        vmr_ci_lower,
        vmr_ci_upper,
        mc_ci_lower,
        mc_ci_upper,
        patch_ci_lower,
        patch_ci_upper,
        d_ci_lower,
        d_ci_upper,
        h_ci_lower,
        h_ci_upper,
        k_ci_lower,
        k_ci_upper,
        Observation = dplyr::case_when(
          nH == 0 ~ "Not analyzed",
          nH == 1 ~ "One host examined",
          nH_inf == 0 ~ "No hosts infected",
          nH_inf == 1 ~ "One host infected",
          TRUE ~ "Multiple hosts infected"
        )
      )
    results_list[[i]] <- result_row
  }

  # Combinar resultados
  aggregation_results <- dplyr::bind_rows(results_list)

  # ---------------------------
  # REDONDEAR RESULTADOS
  # ---------------------------
  aggregation_results <- aggregation_results |>
    dplyr::mutate(
      VMR = round(VMR, decimal_places),
      mean_crowding = round(mean_crowding, decimal_places),
      patchiness = round(patchiness, decimal_places),
      poulin_D = round(poulin_D, decimal_places),
      hoover = round(hoover, decimal_places),
      k = round(k, 3),
      vmr_ci_lower = round(vmr_ci_lower, decimal_places),
      vmr_ci_upper = round(vmr_ci_upper, decimal_places),
      mc_ci_lower = round(mc_ci_lower, decimal_places),
      mc_ci_upper = round(mc_ci_upper, decimal_places),
      patch_ci_lower = round(patch_ci_lower, decimal_places),
      patch_ci_upper = round(patch_ci_upper, decimal_places),
      d_ci_lower = round(d_ci_lower, decimal_places),
      d_ci_upper = round(d_ci_upper, decimal_places),
      h_ci_lower = round(h_ci_lower, decimal_places),
      h_ci_upper = round(h_ci_upper, decimal_places),
      k_ci_lower = round(k_ci_lower, 3),
      k_ci_upper = round(k_ci_upper, 3)
    )

  # ---------------------------
  # WARNINGS
  # ---------------------------
  failed_ci <- any(
    is.na(aggregation_results$vmr_ci_lower) |
      is.na(aggregation_results$vmr_ci_upper)
  )

  if (failed_ci) {
    warning("Some confidence intervals could not be estimated. Check if you have at least 3 hosts per group.")
  }

  # ---------------------------
  # COMBINAR CI (opcional)
  # ---------------------------
  if (combine_ci) {
    aggregation_results <- aggregation_results |>
      dplyr::mutate(
        VMR_CI = ifelse(
          nH_inf > 0 & !is.na(vmr_ci_lower),
          paste(vmr_ci_lower, vmr_ci_upper, sep = " - "),
          NA_character_
        ),
        MC_CI = ifelse(
          nH_inf > 0 & !is.na(mc_ci_lower),
          paste(mc_ci_lower, mc_ci_upper, sep = " - "),
          NA_character_
        ),
        patch_CI = ifelse(
          nH_inf > 0 & !is.na(patch_ci_lower),
          paste(patch_ci_lower, patch_ci_upper, sep = " - "),
          NA_character_
        ),
        D_CI = ifelse(
          nH_inf > 0 & !is.na(d_ci_lower),
          paste(d_ci_lower, d_ci_upper, sep = " - "),
          NA_character_
        ),
        H_CI = ifelse(
          nH_inf > 0 & !is.na(h_ci_lower),
          paste(h_ci_lower, h_ci_upper, sep = " - "),
          NA_character_
        ),
        k_CI = ifelse(
          nH_inf > 0 & !is.na(k_ci_lower),
          paste(k_ci_lower, k_ci_upper, sep = " - "),
          NA_character_
        )
      ) |>
      dplyr::select(
        dplyr::all_of(group_vars),
        Sp,
        nH,
        nH_inf,
        VMR,
        VMR_CI,
        mean_crowding,
        MC_CI,
        patchiness,
        patch_CI,
        poulin_D,
        D_CI,
        hoover,
        H_CI,
        k,
        k_CI,
        Observation
      )
    aggregation_results<-   aggregation_results %>%
      dplyr::rename(!! paste("vmr_CI", ci_method, sep="_") := VMR_CI,
                    !! paste("mc_CI", ci_method, sep="_") := MC_CI,
                    !! paste("patch_CI", ci_method, sep="_") := patch_CI,
                    !! paste("d_CI", ci_method, sep="_") := D_CI,
                    !! paste("h_CI", ci_method, sep="_") := H_CI,
                    k_CI_ML := k_CI)
  } else {
    aggregation_results <- aggregation_results |>
      dplyr::select(
        dplyr::all_of(group_vars),
        Sp,
        nH,
        nH_inf,
        VMR,
        vmr_ci_lower,
        vmr_ci_upper,
        mean_crowding,
        mc_ci_lower,
        mc_ci_upper,
        patchiness,
        patch_ci_lower,
        patch_ci_upper,
        poulin_D,
        d_ci_lower,
        d_ci_upper,
        hoover,
        h_ci_lower,
        h_ci_upper,
        k,
        k_ci_lower,
        k_ci_upper,
        Observation
      )
    aggregation_results<-   aggregation_results %>%
      dplyr::rename(!! paste("vmr_Lower", ci_method, sep="_") := vmr_ci_lower,
                    !! paste("vmr_Upper", ci_method, sep="_") := vmr_ci_upper,
                    !! paste("mc_Lower", ci_method, sep="_") := mc_ci_lower,
                    !! paste("mc_Upper", ci_method, sep="_") := mc_ci_upper,
                    !! paste("patch_Lower", ci_method, sep="_") := patch_ci_lower,
                    !! paste("patch_Upper", ci_method, sep="_") := patch_ci_upper,
                    !! paste("d_Lower", ci_method, sep="_") := d_ci_lower,
                    !! paste("d_Upper", ci_method, sep="_") := d_ci_upper,
                    !! paste("h_Lower", ci_method, sep="_") := h_ci_lower,
                    !! paste("h_Upper", ci_method, sep="_") := h_ci_upper,
                    k_Lower_ML := k_ci_lower,
                    k_Upper_ML := k_ci_upper)
  }

  if (verbose) message("Calculation completed")
  if (save_csv){
    current_date <- Sys.Date()
    file_name <- paste0("aggregation_results_", current_date, ".csv")
    readr::write_csv(aggregation_results, file_name)
    if (verbose) message(".CSV file created in current working directory")
  }
  return(aggregation_results)
}