R/shannon_frag.R

In diversityArch: Computes Diversity Indices for Archaeological Data

#' @title Shannon diversity decomposition
#'
#' @description Computes Shannon diversity and its decomposition
#'
#' @param x Vector of dimension S (number of species) with the number of individuals observed in each species. NA values are allowed. 0 values are converted to NA.
#' @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.
#'
#' @return
#' \itemize{
#'   \item \code{shannon}: Shannon's total Entropy.
#'   \item \code{within}: Within groups entropy.
#'   \item \code{between}: Between groups entropy.
#'   \item \code{groups}: A data frame with information about each group: relative frequency, internal entropy and number of species.
#' }
#'
#' @seealso \link{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)
#' data(Qesem_f)
#' shannon_frag(Qesem_s$HU, Qesem_s$Group, Qesem_f$HU, Qesem_f$Group)
#'
#' @export
shannon_frag <- function(x, gx, f, gf){
  stopifnot(length(x)==length(as.character(gx)),
            length(f)==length(levels(gx)),
            length(f)==length(as.character(gf)),
            is.factor(gx),
            is.factor(gf),
            is.numeric(f),
            all(f >= 0),
            isTRUE(all.equal(f, round(f)))
  )
  x <- ifelse(x==0, NA, x)
  xx <- stats::na.omit(data.frame(x, gx))
  stopifnot(
    is.numeric(xx$x),
    all(xx$x > 0),
    isTRUE(all.equal(xx$x, round(xx$x)))
  )

  d <- dec_shannon(xx$x, xx$gx)
  g <- d$groups
  l <- levels(xx[,2])
  k <- length(l)
  y <- numeric(k)
  Hg <- numeric(k)
  for (i in 1:k){
    y[i] <- as.numeric(sum(xx[xx[,2]==l[i],1]) + f[gf==l[i]])
    Hg[i] <- g[g[,1]==l[i],3]
  }
  P <- y/sum(y)
  b <- -sum(P*log(P))
  w <- sum(P*Hg)
  gg <- data.frame(group=l, Prop=P, Entropy=Hg, S=g$S)
  rownames(gg) <- NULL

  return(list(shannon = w + b,
              within = w,
              between = b,
              groups = gg))
}