R/beta_functions.R

In T-Engel/hammeRs:

#' Calculate beta diveristies of S_cov, S_N and S_PIE
#'
#'
#' Beta diversities quantify spatial species turnover of samples belonging to the sama gamma scale/ group. beta_S_PIE characterises spatial turnover of dominant species (i.e. the base of the individual based rarefaction curve).
#' Beta diversities of S_N and S_cov reflect spatial turnover of rare species.
#'
#' @param mob_in A mob_in object that contains the commnity matrix and the env dataframe eith a grouping variable specified in group_var.
#' @param group_var A character string referring to the column in mob_in$env that groups the alpha-scale samples into gammas.
#' @param effort_samples The number of individuals N used for beta_S_N. If not specified, it will automatically use the recommended sample size. That is the number of individuals in the smallest sample if extrapolate = FALSE or twice as many if extrapolate = TRUE.
#' @param stand_cov The coverage used for coverage based rarefaction. If not specified, it will automatically use the recommended sample coverage. That is the lowest coverage value of the alpha samples if extrapolate = FALSE or the expected coverage of that sample if it was extrapolated to 2N. See recommended_C() for details.
#' @param extrapolate A boolean vecotor of length 1 or 2 specifing whether extrapolation should be used for individial- and coverage-based richness standardisation. The first value in the vector corresponds to individual-based standardisation, the second to coverage-based standardisation. If a single value is supplied, it is adopted for both.
#' @param na.rm Boolean. Passed on to the function that aggregates the alpha diversities to a mean.
#'
#' @export
#'
#' @examples
#' library(mobr)
#' data("inv_comm")
#' data("inv_plot_attr")
#' inv_mob_in<-make_mob_in(inv_comm, inv_plot_attr)
#' betas<-beta
beta_diversity = function(mob_in, group_var,
                         effort_samples = NULL, stand_cov=NULL,
                         extrapolate = c("N"=TRUE, "cov"=TRUE),
                         na.rm= T){

  require(mobr)
  require(tidyverse)
  if (length(extrapolate)==1) extrapolate<- rep(extrapolate,2)
  names(extrapolate)<-c("N", "cov")
  N<-rowSums(mob_in$comm)
  if(is.null(effort_samples))effort_samples= as.numeric(ifelse(extrapolate["N"]==T,2*min(N[N>0]),min(N[N>0])))

  if(effort_samples> min(N) & extrapolate["N"]==F) stop(paste("effort_samples must not exceed the  smallest sample size (here: ",
                                                              min(N), ") if extrapolate is FALSE!", sep = ""))

  if(is.null(stand_cov)) stand_cov= as.numeric(ifelse(extrapolate["cov"]==T, recommended_C(mob_in$comm), min(coverage(mob_in$comm),na.rm = T)))

  if(stand_cov> min(coverage(mob_in$comm),na.rm = T) & extrapolate["cov"]==F) stop(paste("stand_cov must not exceed the coverage of the smallest sample (here: ",
                                                                                         min(coverage(mob_in$comm),na.rm = T), ") if extrapolate is FALSE!", sep = ""))
  if(effort_samples==2) warning("Effort_samples is 2. This is too low to draw any meaningful conclusions in contrasting beta_S_N and beta_S_PIE")

  # alphas
  alphas=data.frame(
    group=factor(mob_in$env[, group_var]),
    N=N,
    S=vegan::specnumber(mob_in$comm),
    S_n= apply(mob_in$comm, 1, mobr::rarefaction, method = "indiv",effort = effort_samples,extrapolate = extrapolate[1],quiet=T),
    S_PIE= calc_PIE(mob_in$comm,ENS = T),
    S_cov=rarefy_coverage(m=mob_in$comm,C = stand_cov,q = 0))

  alphas <- alphas %>% group_by(group) %>% summarise_all(mean, na.rm= na.rm)

  # Abundance distribution pooled in groups
  group_id  = factor(mob_in$env[, group_var])
  abund_group = aggregate(mob_in$comm, by = list(Group=group_id), FUN = "sum")
  rownames(abund_group)<- abund_group$Group
  abund_group= as.matrix(abund_group[,-1])

  # gammas
  gammas=data.frame(
    group=factor(rownames(abund_group)),
    N=rowSums(abund_group),
    S=vegan::specnumber(abund_group),
    S_n= apply(abund_group, 1, mobr::rarefaction, method = "indiv",effort = effort_samples,extrapolate = extrapolate[1],quiet=T),
    S_PIE= calc_PIE(abund_group,ENS = T),
    S_cov=rarefy_coverage(m=abund_group,C = stand_cov,q = 0))

  # merge
  out<-alphas %>% left_join(gammas, by = "group", suffix= c("_alpha", "_gamma"))

  # betas
  out<-out %>%mutate(beta_S= S_gamma/S_alpha,
                     beta_S_n= S_n_gamma/S_n_alpha,
                     beta_S_PIE= S_PIE_gamma/S_PIE_alpha,
                     beta_S_cov= S_cov_gamma/S_cov_alpha,
                     N_stand= effort_samples,
                     C_stand= stand_cov
                     )
  return(out)
}