Add beta_C code

NAMESPACE

@@ -4,7 +4,10 @@ S3method(plot,mob_out)
 S3method(plot,mob_stats)
 S3method(print,mob_in)
 S3method(subset,mob_in)
+export(C_target)
+export(Chat)
 export(avg_nn_dist)
+export(beta_C)
 export(calc_PIE)
 export(calc_SPIE)
 export(calc_chao1)
@@ -14,6 +17,7 @@ export(compare_samp_rarefaction)
 export(get_delta_stats)
 export(get_mob_stats)
 export(get_null_comm)
+export(invChat)
 export(kNCN_average)
 export(make_mob_in)
 export(plot_N)
@@ -23,6 +27,7 @@ export(rarefaction)
 import(dplyr)
 import(ggplot2)
 import(purrr)
+import(stats)
 importFrom(egg,ggarrange)
 importFrom(geosphere,centroid)
 importFrom(geosphere,midPoint)

R/beta_C.R

@@ -0,0 +1,208 @@
+#' Calculate expected sample coverage C_hat
+#'
+#' Returns expected sample coverage  of a sample x for a smaller than observed sample size m (Chao & Jost, 2012).
+#' This code was copied from INEXT's internal function Chat.Ind() (Hsieh et al 2016).
+#'
+#' @param x integer vector (species abundances)
+#' @param m integer. (smaller than observed sample size)
+#'
+#' @return a numeric value.
+#' @export
+#'
+#' @examples
+#' \donttest{
+#' library(vegan)
+#' data(BCI)
+#'
+#' # What is the expected coverage corresponding to a sample size of 50 at the gamma scale?
+#' Chat(colSums(BCI), 50)
+#' }
+Chat <- function (x, m)
+{
+    x <- x[x > 0]
+    n <- sum(x)
+    f1 <- sum(x == 1)
+    f2 <- sum(x == 2)
+    f0.hat <- ifelse(f2 == 0, (n - 1) / n * f1 * (f1 - 1) / 2, (n -
+                                                                    1) / n * f1 ^
+                         2 / 2 / f2)
+    A <- ifelse(f1 > 0, n * f0.hat / (n * f0.hat + f1), 1)
+    Sub <- function(m) {
+        if (m < n) {
+            xx <- x[(n - x) >= m]
+            out <- 1 - sum(xx / n * exp(
+                lgamma(n - xx + 1) - lgamma(n -
+                                                xx - m + 1) - lgamma(n) + lgamma(n - m)
+            ))
+        }
+        if (m == n)
+            out <- 1 - f1 / n * A
+        if (m > n)
+            out <- 1 - f1 / n * A ^ (m - n + 1)
+        out
+    }
+    sapply(m, Sub)
+}
+
+#' Number of individuals corresponding to a desired coverage (inverse C_hat)
+#'
+#' If you wanted to resample a vector to a certain expected sample coverage, how many individuals would you have to draw?
+#' This is C_hat solved for the number of individuals. This code is a modification of INEXT's internal function invChat.Ind() (Hsieh et al 2016).
+#'
+#' @param x integer vector.
+#' @param C numeric. between 0 and 1
+#'
+#' @return a numeric value
+#' @export
+#' @import stats
+#' @examples
+#' \donttest{
+#' library(vegan)
+#' data(BCI)
+#'
+#' # What sample size corresponds to an expected sample coverage of 55%?
+#' invChat(colSums(BCI), 0.55)
+#' }
+#'
+invChat <- function (x, C)
+{
+    m <- NULL
+    n <- sum(x)
+    refC <- Chat(x, n)
+    f <- function(m, C)
+        abs(Chat(x, m) - C)
+    # for interpolation
+    if (refC > C) {
+        opt <- optimize(f,
+                        C = C,
+                        lower = 0,
+                        upper = sum(x))
+        mm <- opt$minimum
+    }
+    # for extrapolation
+    if (refC <= C) {
+        f1 <- sum(x == 1)
+        f2 <- sum(x == 2)
+        if (f1 > 0 & f2 > 0) {
+            A <- (n - 1) * f1 / ((n - 1) * f1 + 2 * f2)
+        }
+        if (f1 > 1 & f2 == 0) {
+            A <- (n - 1) * (f1 - 1) / ((n - 1) * (f1 - 1) + 2)
+        }
+        if (f1 == 1 & f2 == 0) {
+            A <- 1
+        }
+        if (f1 == 0 & f2 == 0) {
+            A <- 1
+        }
+        mm <- (log(n / f1) + log(1 - C)) / log(A) - 1
+        mm <- n + mm
+
+    }
+    if (mm > 2 * n)
+        warning(
+            "The maximum size of the extrapolation exceeds double reference sample size, the results for q = 0 may be subject to large prediction bias."
+        )
+    return(mm)
+}
+
+
+#' Calculate beta_C
+#'
+#' Beta_C uses coverage-based rarefaction to standardize beta-diversity. Calculates the ratio between 
+#' gamma and alpha scale IBR curve for a target gamma-scale sample coverage (i.e. a measure of sample completeness).
+#'
+#' @param x a site by species matrix
+#' @param C target coverage. value between 0 and 1.
+#' @param extrapolation logical. should extrapolation be used?
+#' @param interrupt logical. SHould the function throw an error when C exceeds the maximum recommendable coverage?
+#'
+#' @return a numeric value
+#' @export
+#'
+#' @examples
+#' \donttest{
+#' library(vegan)
+#' data(BCI)
+#'
+#' # What is beta_C for a coverage value of 60%?
+#' beta_C(BCI,C = 0.6)
+#' }
+beta_C <- function(x,
+                   C,
+                   extrapolation = T,
+                   interrupt = T) {
+    x <- as.matrix(x)
+    total <- colSums(x)
+    N <- round(invChat(total, C))
+    C_max = C_target(x, factor = ifelse(extrapolation, 2, 1))
+    if (C > C_max & interrupt == T) {
+        if (extrapolation == F) {
+            stop(
+                paste0(
+                    "Coverage exceeds the maximum possible value for interpolation (i.e. C_target = ",
+                    round(C_max, 4),
+                    "). Use extrapolation or reduce the value of C."
+                )
+            )
+        } else{
+            stop(
+                paste0(
+                    "Coverage exceeds the maximum possible value recommendable for extrapolation (i.e. C_target = ",
+                    round(C_max, 4),
+                    "). Reduce the value of C."
+                )
+            )
+        }
+    }
+    if (N > 1) {
+        gamma_value = rarefaction(x = total,
+                                  method = "IBR",
+                                  effort = N)
+        alpha_value = mean(apply(
+            x,
+            1,
+            rarefaction,
+            method = "IBR",
+            effort = N
+        ))
+        beta = gamma_value / alpha_value
+    } else {
+        beta = NA
+    }
+    attr(beta, "C") = C
+    attr(beta, "N") = N
+    return(beta)
+
+}
+
+#' Calculate the recommended maximum coverage value for the computation of beta_C 
+#'
+#' Returns the estimated gamma-scale coverage that corresponds to the largest allowable sample size 
+#' (i.e. the smallest observed sample size at the alpha scale multiplied by an extrapolation factor).
+#' The default (factor = 2) allows for extrapolation up to 2 times the observed sample size of the 
+#' smallest alpha sample. For factor= 1, only interpolation is applied.
+#' Its not recommendable to use factors larger than 2.
+#'
+#' @param x a site by specie matrix
+#' @param factor numeric. how far do you want to extrapolate? 
+#'
+#' @return numeric value
+#' @export
+#'
+#' @examples
+#' \donttest{
+#' library(vegan)
+#' data(BCI)
+#'
+#' # What is the largest possible C that I can use to calculate beta_C for my site by species matrix?
+#' C_target(BCI)
+#' }
+#'
+
+C_target <- function(x, factor = 2) {
+    x <- as.matrix(x)
+    n = min(factor * rowSums(x))
+    out <- Chat(colSums(x), n)
+    return(out)
+}