R/SAC_spatial.R

In FelixMay/MoBspatial: Spatial Simulation and Scale-Dependent Analysis of Biodiversity Changes

#' Sample species richness
#'
#' Expected species richness in a random sample of fixed size.
#'
#' @param abund_vec Species abundance distribution of the community (integer vector)
#' @param n Sample size in terms of number of individuals (integer)
#'
#' @return Expected number of species in a sample of n individuals
#'
#' @details The expected number of species is calculated after Hurlbert 1971, Equation 3.
#'
#' \code{spec_sample} is similar to the function \code{\link[vegan]{rarefy}}
#' in the R package \code{\link{vegan}}.
#'
#' @references
#' Hurlbert, S.H. 1971. The nonconcept of species diversity: a critique and +
#' alternative parameters. Ecology 52, 577-586.
#'
#' @examples
#' sad1 <- sim_sad(100, 1000)
#' spec_sample(abund_vec = sad1, n = 20)
#'
#' @export
#'
spec_sample <- function(abund_vec, n)
{
   abund_vec <- abund_vec[abund_vec > 0]
   n_total <- sum(abund_vec)
   ldiv <- lchoose(n_total, n)
   p1 <- exp(lchoose(n_total - abund_vec, n) - ldiv)
   out <- sum(1 - p1)
   names(out) <- n

   return(out)
}

#' Species rarefaction curve
#'
#' Expected species richness as a function of sample size
#'
#' @param abund_vec Species abundance distribution of the community (integer vector)
#'
#' @return Numeric Vector with expected species richness in samples of 1, 2, 3 ... n individuals
#'
#' @details This function essentially evaluates \code{\link{spec_sample}} for
#' sample sizes from 1 to \code{sum(abund_vec)}. It is similar to the function
#' \code{\link[vegan:rarefy]{vegan:rarecurve}} in the R package \code{\link{vegan}}.
#'
#' @references
#'
#' Gotelli & Colwell 2001. Quantifying biodiversity: procedures and pitfalls
#' in the measurement and comparison of species richness. Ecology Letters 4, 379--391.
#'
#' @examples
#' sad1 <- sim_sad(100, 2000, sad_type = "lnorm", sad_coef = list("meanlog" = 2,
#'                                                                 "sdlog" = 1))
#' rc1 <- rare_curve(sad1)
#' plot(rc1, type = "l", xlab = "Sample size", ylab = "Expected species richness")
#'
#' @export
#'
rare_curve <- function(abund_vec)
{
   abund_vec <- abund_vec[abund_vec > 0]
   n_total <- sum(abund_vec)
   n_vec <- 1:n_total

   rc <- sapply(n_vec, function(n_sample){spec_sample(abund_vec, n = n_sample)})
   names(rc) <- n_vec
   return(rc)
}

#' Non-spatial and spatially-explicit species sampling curves
#'
#' Expected species richness as function of sample size (no. of individuals),
#' when individuals are sampled randomly (rarefaction) or when nearest-neighbours
#' are samples (accumulation).
#'
#' @param comm Community object
#'
#' @param method Partial match to \code{accumulation} or \code{rarefaction}. Also both
#' methods can be included at the same time.
#'
#' @details
#' Non-spatial sampling corresponds to the species rarefaction curve, which only
#' depends on the species abundance distribution and can thus be also calculated
#' from abundance data (see \code{\link{rare_curve}}).
#'
#' In contrast the species-accumulation curve starts from a focal individual and
#' only samples the nearest neighbours of the focal individual. The final
#' species accumulation curves is calculated as the mean over the accumulation
#' curves starting from all individuals.
#'
#' In contrast to the rarefaction curve the accumulation curve is not only
#' influenced by the species abundance distribution, but also by the spatial
#' distribution of individuals.
#'
#' @return A dataframe with 2-3 columns. The first column indicates the sample
#' size (numbers of individuals), and the second and third column indicate
#' the expected species richness for spatial sampling (column: "spec_accum")
#' and/or random sampling (column "spec_rarefied")
#'
#' @examples
#' sim_com1 <- sim_thomas_community(s_pool = 100, n_sim = 1000)
#' sac1 <- spec_sample_curve(sim_com1, method = c("rare","acc"))
#'
#' head(sac1)
#' plot(sac1)
#'
#' @export
#' @importFrom methods is
#'
spec_sample_curve <- function(comm, method = c("accumulation" ,"rarefaction"))
{
   if (!is(comm, "community"))
      stop("spec_sample_curve requires a community object as input. See ?community.")

   out_dat <- data.frame(n = 1:nrow(comm$census))

   method <- match.arg(method, several.ok = TRUE)

   if ("accumulation" %in% method) {
      out_dat$spec_accum <- sSAC1_C(comm$census$x, comm$census$y,
                                    as.integer(comm$census$species))
   }

   if ("rarefaction" %in% method) {
      abund <- community_to_sad(comm)
      out_dat$spec_rarefied <- rare_curve(abund)

   }

   class(out_dat) <- c("spec_sample_curve", "data.frame")

   return(out_dat)
}


#' Plot species sampling curves
#'
#' @param x Species sampling curve generated by \code{\link{spec_sample_curve}}
#'
#' @param ... Additional graphical parameters used in \code{\link[graphics:plot.default]{graphics::plot}}.
#'
#' @return This function is called for its side effects and has no return value.
#'
#'
#' @examples
#' sim_com1 <- sim_thomas_community(s_pool = 100, n_sim = 1000)
#' sac1 <- spec_sample_curve(sim_com1, method = c("rare","acc"))
#' plot(sac1)
#'
#' @export
#'
plot.spec_sample_curve <- function(x, ...)
{
   graphics::plot(x[[1]], x[[2]], type = "n",
                  xlab = "No. of individuals sampled",
                  ylab = "Expected no.of species",
                  main = "Species sampling curves", ...)

   if (ncol(x) == 2) {
      if (names(x)[2] == "spec_accum") {
         legend_text <- c("Accumulation")
         line_col <- "red"
      }
      if (names(x)[2] == "spec_rarefied") {
         legend_text <- c("Rarefaction")
         line_col <- "blue"
      }
      graphics::lines(x[[1]], x[[2]], col = line_col)
   }

   if (ncol(x) == 3) {
      legend_text <- c("Accumulation","Rarefaction")
      line_col <- c("red","blue")
      graphics::lines(x[[1]], x[[2]], col = line_col[1])
      graphics::lines(x[[1]], x[[3]], col = line_col[2])
   }

   graphics::legend("bottomright", legend = legend_text, lty = 1, col = line_col)
}