R/mobsim_functions.R

In T-Engel/hammeRs:

ENS <- function(com) {
  require(vegan)
  ENS <- 1 / (2 - (rarefy(com, 2)))
  if (is.infinite(ENS) | specnumber(com) == sum(com))
    ENS <- NA
  return(as.numeric(ENS))
}

#' Simulate a sample community of a prescribed S_PIE
#'
#'A wrapper for the mobsim function simsad that internally optimises the CV parameter of the lognormal SAD to generate a community sample of a user specified S_PIE value.
#'
#' @param ENS_target S_PIE value of the sample
#' @param S species richness of the sample
#' @param n_sim number of individuals in the sample. If not specified it will be calculated from the \code{mean_abu} parameter.
#' @param mean_abu mean species abundance. Alternative argument to parameterize number of individuals \code{n_sim}. Default is 10 individuals per species.
#' @param sims number of simulations used for the internal optimization process. Default is 100. Increasing this parameter can significantly slow down the function
#'
#' @return a vector of abundances
#' @export
#'
#' @examples
#' #simulate community of S_PIE=20 with 100 species and 800 individuals
#' comm<- sim_ENS(20,100, 800)
#' mobr::calc_PIE(comm, ENS=T)
sim_ENS <- function(ENS_target,
                    S,
                    n_sim = NULL,
                    mean_abu = 10,
                    sims = 100) {
  require(mobsim)
  ENS <- function(com) {
    require(vegan)
    ENS <- 1 / (2 - (rarefy(com, 2)))
    if (is.infinite(ENS) | specnumber(com) == sum(com))
      ENS <- NA
    return(as.numeric(ENS))
  }
  if (is.null(n_sim))
    n_sim <- S * mean_abu
  seeds = sample(1:(sims * 1000), sims, replace = F)


  testsad <- function(par,
                      S = S,
                      seed = sample(1:1000)) {
    set.seed(seed)
    com <-
      sim_sad(
        s_pool = S,
        n_sim = n_sim,
        sad_type = "lnorm",
        sad_coef = list(cv_abund = par[1]),
        fix_s_sim = T
      )
    return(abs(ENS_target - ENS(com)))
  }
  pars <- sapply(seeds, function(seed) {
    pars_opt <-
      optim(
        par = c(0),
        testsad,
        S = S,
        seed = seed,
        method = "Brent",
        lower = 0,
        upper = S,
        control = list(maxit = 10000)
      )
    return(c(
      seed = seed,
      par = pars_opt$par,
      diff = pars_opt$value
    ))
  })
  pars <- as.data.frame(t(pars))
  pars <- subset(pars, diff == min(pars$diff))
  set.seed(pars$seed)
  com <-
    sim_sad(
      s_pool = S,
      n_sim = n_sim,
      sad_type = "lnorm",
      sad_coef = list(cv_abund = pars$par),
      fix_s_sim = T
    )
  return(as.numeric(com))
}