R/coverage_functions.R

In T-Engel/hammeRs:

#

#' Calculate sample coverage
#'
#' Sample coverage is calculated according to Chao & Jost (2012).
#'
#' @param m site by species matrix
#'
#' @return a vector
#' @export
#'
#' @examples
#' library(vegan)
#' data(BCI)
#' coverage(BCI)
coverage<-function(m){
  if (is.null(dim(m))) m= matrix(m, ncol=length(m))
  n<-apply(m,1, sum)
  f1<-apply(m,1,function(x)length(x[x==1]))
  f2<-apply(m,1,function(x)length(x[x==2]))
  cover=1-(f1/n)*(((n-1)*f1)/((n-1)*f1+2*f2))
  return(cover)
}

#' Calculate maximum recommendable sample coverage
#'
#' Extrapolations of diversity curves should not surpass 2N. This function returns the expected coverage corresponding to a an extrapolated sample size of 2N.
#' I use iNEXT's internal function Chat.Ind() to calculate this.
#'
#'
#' @param m site by species matrix
#'
#' @return a vector
#' @export
#'
#' @seealso [recommended_C()]
#'
#' @examples
#' library(vegan)
#' data(BCI)
#' C_max(BCI)

C_max<-function(m){
  require(iNEXT)
  if (is.null(dim(m))) m= matrix(m, ncol=length(m))
  C_max<-apply(m,1, function(x)iNEXT:::Chat.Ind(x=x ,m = 2*sum(x)))
  return(C_max)
}

#' Recommended coverage
#'
#' Calculates recommended coverage for rarefaction. This is the coverage corresponding to 2N of the sample in m with the lowest observed coverage.
#' I use iNEXT's internal function Chat.Ind() to calculate this.
#'
#'
#' @param m site by species matrix
#' @param na.rm bolean. Passed on to min() function. Shall missing values be removed?
#'
#' @return a vector
#' @export
#'
#' @seealso [C_max()]
#' @examples
#' library(vegan)
#' data(BCI)
#' C_max(BCI)

recommended_C<-function(m, na.rm=T) min(C_max(m),na.rm = na.rm)


#' Calculate rarefied diversity estimates along coverage (e.g. S_cov)
#'
#' This function is similar to iNEXT's estimateD(). It allows interpolation and extrapolation of diversities (of order q) along coverage.
#' It is much faster than estimateD() for individual point estimates and uses data inputs like vegan and mobr (i.e site by species matrix).
#' It uses iNext's internal function invChat.Ind() to do the interpolation/ extrapolation. If no C is NULL,
#' it will default to maximum recommendable sample coverage (using recommended_C()).
#'
#' @param m site by species matrix or single vector
#' @param C Numeric value between 0 and 1 specifying desired coverage level.
#' @param q integer vector of diversity orders
#'
#' @return a vector if q is a single value; a matrix if q is a vector.
#' @export
#' @seealso [recommended_C()], [C_max()]
#' @examples
#' library(vegan)
#' data(BCI)
#' rarefy_coverage(BCI)
rarefy_coverage<- function(m, C=NULL, q=0){
  require(iNEXT)
  require(pbapply)
  m=as.matrix(m)
  if(is.null(C)){
    C=recommended_C(m)
    print(paste("C was not specified by user. Defaulting to maximum recommendable sample coverage = ",
                C,sep=""))
  }
  if (is.null(dim(m))) m= matrix(m, ncol=length(m))
  if(length(C)>1)stop("C must be a single numeric value between 0 and 1. If you need more then 1 coverage levels, use iNext().")
  out<-pbsapply(1:dim(m)[1],function(x) {
    comm<-m[x,]
    comm<- comm[comm>0]
    if(sum(comm)>0){
      a<-iNEXT:::invChat.Ind(comm,C=C)[,paste("q = ",q,sep=""),drop=F]
      }else{
        a<-NA
      }
    names<-colnames(a)
    a<-as.numeric(a)
    names(a)<-names
    return(a)
  }
  )
  if (is.null(dim(out))==F){
    colnames(out)=rownames(m)
    out=t(out)} else{
      names(out)=rownames(m)
    }

  return(out)

}


#' Calculate sample coverage of diversity order 2
#'
#' I define second order coverage as the complement of the slope of the diversity curve for order 2 at the given sample size N.
#' Practically, it is computed as the slope  of the second order diversity curve between k=N and k= N+1.
#'
#'
#' @param m site by species matrix
#'
#' @return a vector
#' @export
#'
#' @examples
#' library(vegan)
#' data(BCI)
#' coverage2(BCI)
coverage2<-function(m){
  require(vegan)
  if (is.null(dim(m))) m= matrix(m, ncol=length(m))
  N= rowSums(m)
  simpson=diversity(m, index = "simpson")
  cov2= ifelse(N>1, 1+1/(1-simpson)-1/(1-N^2/(N^2-1)*simpson),NA)
  return(cov2)
}