R/bs_dec.R

In diversityArch: Computes Diversity Indices for Archaeological Data

#' @title Bootstrap estimates of the decomposed Shannon's diversity index
#'
#' @description Computes bootstrap confidence intervals of the Shannon's diversity index
#'
#' @param x Vector of dimension S (number of species) with the number of individuals observed in each species. NA values are allowed.
#' @param groups Vector of dimension S of factors indicating the groups of each species.
#' @param B Number of bootstrap samples. The default is 1000.
#' @param cl Confidence level. A value between 0 and 1. The default is 0.95.
#'
#' @return
#' \itemize{
#'   \item \code{ci_bt}: Confidence interval for "between groups".
#'   \item \code{ci_wt}: Confidence interval for "within groups".
#' }
#'
#' @seealso \link{shannon}, \link{dec_shannon}
#'
#' @references
#' Arnaud Barat, Andreu Sansó, Maite Arilla-Osuna, Ruth Blasco, Iñaki Pérez-Fernández, Gabriel Cifuentes-Alcobenda, Rubén Llorente, Daniel Vivar-Ríos, Ella Assaf, Ran Barkai, Avi Gopher, & Jordi Rosell-Ardèvol (2026) <doi:10.1007/s10816-026-09802-3>.
#'
#' @examples
#' data(Qesem_s)
#' bs_dec(Qesem_s$HU, Qesem_s$Group)
#'
#' @export
bs_dec <- function(x, groups, B = 1000, cl = 0.95){
  stopifnot(is.numeric(B),
            B==round(B),
            cl>0 & cl<1,
            length(x)==length(as.character(groups))
            )
  groups <- as.factor(groups)
  ok <- !is.na(x) & x > 0 & !is.na(groups)
  x <- x[ok]
  groups <- droplevels(groups[ok])

  n <- sum(x)
  p <- x / n

  decomp_once <- function(counts, groups) {
    ok <- !is.na(counts) & counts > 0 & !is.na(groups)
    counts <- counts[ok]
    groups <- droplevels(groups[ok])

    p <- counts / sum(counts)
    P <- tapply(p, groups, sum)

    # between-group component
    bt <- -sum(P[P > 0] * log(P[P > 0]))

    # within-group component
    wt <- 0
    levs <- names(P)
    for (lv in levs) {
      idx <- groups == lv
      Pg <- P[lv]
      if (Pg > 0) {
        q <- p[idx] / Pg
        wt <- wt - Pg * sum(q[q > 0] * log(q[q > 0]))
      }
    }

    c(bt, wt)
  }

  res <- matrix(NA_real_, nrow = 2, ncol = B)
  rownames(res) <- c("bt", "wt")

  # bootstrap replications
  for (b in seq_len(B)) {
    counts_b <- as.vector(stats::rmultinom(1, size = n, prob = p))
    res[, b] <- decomp_once(counts_b, groups)
  }

  bt_valid <- res["bt", is.finite(res["bt", ])]
  wt_valid <- res["wt", is.finite(res["wt", ])]

  alpha <- 1 - cl
  probs <- sort(c(alpha / 2, 1 - alpha / 2))

  return(list(ci_bt = unname(stats::quantile(bt_valid, probs = probs, na.rm = TRUE)),
              ci_wt = unname(stats::quantile(wt_valid, probs = probs, na.rm = TRUE))
              )
         )
}