Outsource uncertainty

CITATION.cff

@@ -8,7 +8,7 @@ message: 'To cite package "b3gbi" in publications use:'
 type: software
 license: MIT
 title: 'b3gbi: General Biodiversity Indicators for Biodiversity Data Cubes'
-version: 0.8.12
+version: 0.9.0
 abstract: Calculate general biodiversity indicators from GBIF data cubes. Includes
   many common indicators such as species richness and evenness, which can be calculated
   over time (trends) or space (maps).
@@ -78,6 +78,23 @@ references:
     orcid: https://orcid.org/0000-0003-4777-038X
   year: '2025'
   doi: 10.32614/CRAN.package.dplyr
+- type: software
+  title: dubicube
+  abstract: 'dubicube: Calculation and Interpretation of Data Cube Indicator Uncertainty'
+  notes: Imports
+  url: https://b-cubed-eu.github.io/dubicube/
+  authors:
+  - family-names: Langeraert
+    given-names: Ward
+    email: ward.langeraert@inbo.be
+    orcid: https://orcid.org/0000-0002-5900-8109
+    affiliation: Research Institute for Nature and Forest (INBO)
+  - family-names: Van Daele
+    given-names: Toon
+    email: toon.vandaele@inbo.be
+    orcid: https://orcid.org/0000-0002-1362-853X
+    affiliation: Research Institute for Nature and Forest (INBO)
+  year: '2025'
 - type: software
   title: ggplot2
   abstract: 'ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics'

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: b3gbi
 Title: General Biodiversity Indicators for Biodiversity Data Cubes
-Version: 0.8.12
+Version: 0.9.0
 Authors@R: c(
     person("Shawn", "Dove", email = "Shawn.Dove@allzool.bio.uni-giessen.de",
            role = c("aut", "cre"), comment = c(ORCID = "0000-0001-9465-5638")),
@@ -19,6 +19,7 @@ Depends:
 Imports:
     boot,
     dplyr,
+    dubicube,
     ggplot2,
     iNEXT,
     labeling,
@@ -60,7 +61,8 @@ Remotes:
     github::hrbrmstr/mgrs,
     github::ropensci/bold,
     github::ropensci/taxize,
-    github::ropensci/rnaturalearthhires
+    github::ropensci/rnaturalearthhires,
+    github::b-cubed-eu/dubicube
 BugReports: https://github.com/b-cubed-eu/b3gbi/issues
 Funding: This project receives funding from the European Union’s Horizon Europe 
     Research and Innovation Programme (ID No 101059592).

NAMESPACE

@@ -67,6 +67,7 @@ S3method(print,sim_cube)
 export("%>%")
 export(ab_rarity_map)
 export(ab_rarity_ts)
+export(add_ci)
 export(area_rarity_map)
 export(area_rarity_ts)
 export(calc_ci)

NEWS.md

@@ -1,3 +1,9 @@
+# b3gbi 0.9.0 - Major update:
+
+* Confidence intervals are no longer calculated in the core indicator workflow, except for Hill diversity (Hill diversity is handled by the iNEXT package, which calculates confidence intervals internally). They are now calculated using a separaed function, add_ci(), which can be applied to any indicator_ts object for which reasonable confidence intervals can be calculated. This gives the user more freedom. After adding CIs, you can recalculate with different parameters by setting 'replace = TRUE' when you call add_ci().
+* Bootstrapping for confidence intervals is now done across the entire cube by default. This uses the dubicube package, which is now added as a dependency. The option to calculate at indicator level is still available by setting 'level = "indicator"' when calling add_ci(). Indicator level bootstrapping is a faster but less robust method.
+* The default number of bootstrap replicates when calculating confidence intervals is now 1000. This improves robustness at the sake of speed. This can still be changed using the 'num_bootstrap' parameter when calling add_ci().
+
 # b3gbi 0.8.12 - Minor update:
 
 * The cellCode column from the cube is now retained in indicator map outputs. This allows users to trace back grid cells to their original codes in the cube.

R/add_ci.R

@@ -0,0 +1,207 @@
+#' Add Confidence Intervals to an Indicator Object
+#'
+#' @description
+#' This function calculates bootstrap confidence intervals for an existing
+#' `indicator_ts` object. It supports both cube-level bootstrapping (resampling
+#' occurrence records) and indicator-level bootstrapping (resampling
+#' calculated values), allowing for advanced transformations during the
+#' CI calculation process.
+#'
+#' @param indicator An object of class `indicator_ts` to which confidence
+#'   intervals should be added.
+#' @param num_bootstrap (Optional) Number of bootstrap replicates to perform.
+#'   (Default: 1000)
+#' @param bootstrap_level (Optional) Level at which to perform bootstrapping:
+#'   * `cube` (default): Bootstrapping is done by resampling the
+#'     occurrence records in the cube. This is mathematically more robust as it
+#'     captures the underlying sampling uncertainty.
+#'   * `indicator`: Bootstrapping is done by resampling indicator
+#'     values. This is faster for large cubes but less robust.
+#'
+#' @param ci_type (Optional) Type of bootstrap confidence intervals to
+#'   calculate. (Default: `"norm"`). Supported options are:
+#'   * `norm`: Normal approximation intervals.
+#'   * `basic`: Basic bootstrap intervals.
+#'   * `perc`: Percentile intervals.
+#'   * `bca`: Bias-corrected and accelerated intervals.
+#'   * `none`: No confidence intervals calculated.
+#'
+#' @param trans (Optional) A function for transforming the indicator values
+#'   before calculating confidence intervals (e.g., `log`).
+#'   (Default: identity function)
+#' @param inv_trans (Optional) The inverse of the transformation function
+#'   `trans` (e.g., `exp`). Used to back-transform the intervals
+#'   to the original scale. (Default: identity function)
+#' @param confidence_level (Optional) The confidence level for the calculated
+#'   intervals (e.g., 0.95 for 95% CIs). (Default: 0.95)
+#' @param overwrite (Optional) Logical. If the indicator already contains
+#'   confidence intervals (`ll` and `ul` columns), should they
+#'   be replaced? (Default: TRUE)
+#' @param ... Additional arguments passed to the underlying functions.
+#'
+#' @details
+#' The function acts as a bridge to the `dubicube` package for statistical
+#' heavy lifting. For certain indicators (e.g., Hill numbers), confidence
+#' intervals cannot be added post-hoc as they are calculated internally by
+#' the `iNext` package during the initial calculation. In such cases,
+#' a warning is issued and the original object is returned.
+#'
+#' @return An updated object of class `indicator_ts` containing the
+#'   original data with the following additional columns:
+#'   * `ll`: Lower limit of the confidence interval.
+#'   * `ul`: Upper limit of the confidence interval.
+#'   * `est_boot`: The bootstrap estimate of the indicator value.
+#'   * `se_boot`: The bootstrap standard error.
+#'   \item `bias_boot`: The bootstrap estimate of bias.
+#'   \item `int_type`: The type of interval calculated (e.g., 'perc').
+#'   \item `conf`: The confidence level used.
+#'
+#' @seealso [dubicube::bootstrap_cube()], [dubicube::calculate_bootstrap_ci()]
+#'
+#' @export
+add_ci <- function(indicator,
+                   num_bootstrap = 1000,
+                   bootstrap_level = c("cube",
+                                       "indicator"),
+                   ci_type = c("norm",
+                               "basic",
+                               "perc",
+                               "bca",
+                               "none"),
+                   trans = function(t) t,
+                   inv_trans = function(t) t,
+                   confidence_level = 0.95,
+                   overwrite = TRUE,
+                   ...) {
+
+  # Check for correct object class
+  if (!inherits(indicator, "indicator_ts")) {
+    stop("indicator must be an indicator_ts object.")
+  }
+
+  ll <- ul <- year <- est_original <- NULL
+
+  # List of indicators for which bootstrapped confidence intervals should not
+  # be calculated
+  noci_list <- c("obs_richness",
+                 "cum_richness",
+                 "occ_turnover",
+                 "tax_distinct",
+                 "hill0",
+                 "hill1",
+                 "hill2")
+
+  # Match ci_type argument
+  ci_type <- match.arg(ci_type)
+  bootstrap_level <- match.arg(bootstrap_level)
+
+  # If indicator is in noci_list, return indicator without calculating CIs
+  if (indicator$div_type %in% noci_list) {
+    if (indicator$div_type %in% c("hill0", "hill1", "hill2")) {
+      warning(
+        paste0(
+          "Confidence intervals cannot calculated for ",
+          indicator$div_type,
+          " as they are handled by the iNext package when calculating
+          your indicator. Returning indicator without adding CIs."
+        )
+      )
+    } else
+      warning(
+        paste0(
+          "Cannot calculate sensible confidence intervals for ",
+          indicator$div_type, ". Returning indicator without CIs."
+        )
+      )
+    return(indicator)
+  }
+
+  # Extract data from indicator object
+  x <- indicator$data
+  raw_data <- indicator$raw_data
+
+  if(any(c("ll", "ul") %in% names(x)) & !overwrite) {
+    warning(
+      paste0(
+        "Indicator already contains confidence intervals. Returning indicator
+        without adding CIs. Use 'replace = TRUE' argument to recalculate CIs."
+      )
+    )
+    return(indicator)
+  }
+
+  # Remove existing confidence intervals if overwrite = TRUE
+  if (overwrite &
+      all(c("ll", "ul") %in% names(x))) {
+    x <- x %>%
+      dplyr::select(-ll, -ul)
+  }
+
+  # Calculate confidence intervals
+  if (bootstrap_level == "indicator") {
+
+    # Send data to calc_ci for indicator level bootstrapping
+    indicator <- calc_ci(raw_data,
+                         indicator = indicator,
+                         num_bootstrap = num_bootstrap,
+                         ci_type = ci_type,
+                         ...)
+
+  } else if (bootstrap_level == "cube") {
+
+    # Identify indicators that require species-level grouping
+    species_level_indicators <- c("spec_occ", "spec_range")
+
+    if (indicator$div_type %in% species_level_indicators) {
+      group_cols <- c("year", "taxonKey")
+    } else {
+      group_cols <- "year"
+    }
+
+    # Bootstrap cube data
+    bootstrap_results <- dubicube::bootstrap_cube(
+      data_cube = raw_data,
+      fun = calc_ts,
+      grouping_var = group_cols,
+      samples = num_bootstrap,
+      seed = 123,
+      progress = TRUE,
+      processed_cube = FALSE,
+      #method = "whole_cube",
+    )
+
+    # Calculate confidence intervals from bootstrap results
+    ci_df <- dubicube::calculate_bootstrap_ci(
+      bootstrap_samples_df = bootstrap_results,
+      grouping_var = group_cols,
+      type = ci_type,
+      h = trans,
+      hinv = inv_trans,
+      conf = confidence_level,
+      data_cube = raw_data,
+      fun = calc_ts
+    ) %>%
+      dplyr::select(-est_original)
+
+    # Join confidence intervals to indicator object
+    if (nrow(ci_df) > 0) {
+      # Convert negative values to zero as rarity cannot be less than zero
+      ci_df$ll <- ifelse(ci_df$ll > 0, ci_df$ll, 0)
+      # Join confidence intervals to indicator values by year
+      x <- x %>%
+        dplyr::full_join(ci_df,
+                         by = group_cols)
+      indicator$data <- x
+      return(indicator)
+    } else {
+      warning(
+        paste0(
+          "Unable to calculate confidence intervals. There may be ",
+          "insufficient data."
+        )
+      )
+    }
+  } else {
+    stop("Invalid bootstrap_level. Choose 'cube' or 'indicator'.")
+  }
+}

R/calc_map_evenness_core.R

@@ -7,7 +7,7 @@ calc_map_evenness_core <- function(x,
                   internal function, not meant to be called directly.",
                   inherits(x, c("data.frame", "sf")))
 
-  num_occ <- obs <- cellid <- taxonKey <- cellCode <- . <- NULL
+  num_occ <- obs <- cellid <- taxonKey <- diversity_val <- cellCode <- . <- NULL
 
   type <- match.arg(type, names(available_indicators))