R/bs_frag.R

In diversityArch: Computes Diversity Indices for Archaeological Data

#' @title Bootstrap estimates of the decomposed Shannon's diversity index with fragments
#'
#' @description Computes bootstrap confidence intervals of the decomposed Shannon's diversity index using evidence from fragments
#'
#' @param x Vector of dimension S (number of species) with the number of individuals observed in each species. NA values are allowed.
#' @param gx Vector of dimension S of factors indicating the group of each species. G groups.
#' @param f Vector of dimension G with the number (>0) of fragments in each group
#' @param gf Vector of dimension G of factors indicating the groups to which the fragments "f" belong.
#' @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_sh}: Confidence interval for Shannon's diversity index with fragments.
#'   \item \code{ci_bt}: Confidence interval for "between groups".
#'   \item \code{ci_wt}: Confidence interval for "within groups".
#' }
#'
#' @seealso \link{shannon}, \link{dec_shannon}, \link{shannon_frag}
#'
#' @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)
#' data(Qesem_f)
#' bs_frag(Qesem_s$HU, Qesem_s$Group, Qesem_f$HU, Qesem_f$Group)
#'
#' @export
bs_frag <- function(x, gx, f, gf, B = 1000, cl = 0.95){
  stopifnot(is.numeric(B),
            B==round(B),
            cl>0 & cl<1,
            length(x)==length(as.character(gx)),
            length(f)==length(as.character(gf))
  )
  gx <- as.factor(gx)
  ok <- !is.na(x) & x > 0 & !is.na(gx)
  x <- x[ok]
  gx <- droplevels(gx[ok])

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

  gf <- as.factor(gf)
  ok <- !is.na(f) & f > 0 & !is.na(gf)
  f <- f[ok]
  gf <- droplevels(gf[ok])

  Nf <- sum(f)
  Pf <- f / Nf

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

    p <- counts / sum(counts)
    P <- tapply(p, gx, sum)
    Nc <- tapply(counts, gx, sum)

    # fragmentary evidence
    ok <- !is.na(cf) & cf > 0 & !is.na(gf)
    cf <- cf[ok]
    gf <- droplevels(gf[ok])
    b1 <- data.frame(Nc=Nc, g=as.factor(names(Nc)))
    b2 <- data.frame(Nf=cf, g=gf)
    b <- merge(b1, b2, by = "g", all = TRUE)
    b[is.na(b)] <- 0
    Pc <- (b$Nc + b$Nf)/sum(b$Nc + b$Nf)
    names(Pc) <- names(P)

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

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

    c(sh, bt, wt)
  }

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

  # bootstrap replications
  for (b in seq_len(B)) {
    counts <- as.vector(stats::rmultinom(1, size = N, prob = p))
    cf <- as.vector(stats::rmultinom(1, size = Nf, prob = Pf))
    res[, b] <- decomp_once(counts, gx, cf, gf)
  }

  sh_valid <- res["sh", is.finite(res["sh", ])]
  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_sh = unname(stats::quantile(sh_valid, probs = probs, na.rm = TRUE)),
              ci_bt = unname(stats::quantile(bt_valid, probs = probs, na.rm = TRUE)),
              ci_wt = unname(stats::quantile(wt_valid, probs = probs, na.rm = TRUE))
              )
  )
}