R/bs.R

In diversityArch: Computes Diversity Indices for Archaeological Data

#' @title Bootstrap estimates of the diversity indices
#'
#' @description Computes bootstrap standard error and confidence interval of the diversity indices
#'
#' @param x Vector of dimension S (number of species) with the number of individuals observed in each species. NA values are allowed.
#' @param ind Index to be computed. The default value is "shannon". Other possible values are: "simpson_E", "simpson_D", "menhinick", "margalef", "evenness", "equitability" and "dominance".
#' @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{s}: Boostrap standard error.
#'   \item \code{ci}: Confidence interval.
#' }
#'
#' @seealso \link{shannon}, \link{dec_shannon}, \link{dominance}, \link{equitability}, \link{evenness}, \link{margalef}, \link{menhinick}, \link{simpson_D}, \link{simpson_E}
#'
#' @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(Qesem_s$HU)
#'
#' @export
bs <- function(x, ind = "shannon", B = 1000, cl = 0.95){

  x <- check_x(x)
  stopifnot(is.numeric(B),
            B==round(B),
            cl>0 & cl<1,
            ind %in% c("shannon", "simpson_E", "simpson_D", "menhinick",
                       "margalef", "evenness", "equitability", "dominance")
            )

  n <- sum(x)
  p <- x/n
  m <- length(x)
  y <- numeric(B)
  for (i in 1:B){
    xx <- sample(1:m, n, replace = T, prob = p)
    xb <- as.vector(table(xx))
    y[i] <- switch(ind,
                   shannon = shannon(xb),
                   simpson_E = simpson_E(xb),
                   simpson_D = simpson_D(xb),
                   menhinick = menhinick(xb),
                   margalef = margalef(xb),
                   evenness = evenness(xb),
                   equitability = equitability(xb),
                   dominance = dominance(xb)
                   )
  }

  alpha <- 1-cl
  q <- unname(stats::quantile(y, probs = sort(c(alpha/2, 1-alpha/2))))

  return(list(s = stats::sd(y), ci = q)          )
}