R/bootstraping.R

In poppr: Genetic Analysis of Populations with Mixed Reproduction

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
#
# This software was authored by Zhian N. Kamvar and Javier F. Tabima, graduate
# students at Oregon State University; Jonah C. Brooks, undergraduate student at
# Oregon State University; and Dr. Nik Grünwald, an employee of USDA-ARS.
#
# Permission to use, copy, modify, and distribute this software and its
# documentation for educational, research and non-profit purposes, without fee,
# and without a written agreement is hereby granted, provided that the statement
# above is incorporated into the material, giving appropriate attribution to the
# authors.
#
# Permission to incorporate this software into commercial products may be
# obtained by contacting USDA ARS and OREGON STATE UNIVERSITY Office for
# Commercialization and Corporate Development.
#
# The software program and documentation are supplied "as is", without any
# accompanying services from the USDA or the University. USDA ARS or the
# University do not warrant that the operation of the program will be
# uninterrupted or error-free. The end-user understands that the program was
# developed for research purposes and is advised not to rely exclusively on the
# program for any reason.
#
# IN NO EVENT SHALL USDA ARS OR OREGON STATE UNIVERSITY BE LIABLE TO ANY PARTY
# FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING
# LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION,
# EVEN IF THE OREGON STATE UNIVERSITY HAS BEEN ADVISED OF THE POSSIBILITY OF
# SUCH DAMAGE. USDA ARS OR OREGON STATE UNIVERSITY SPECIFICALLY DISCLAIMS ANY
# WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
# MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE AND ANY STATUTORY
# WARRANTY OF NON-INFRINGEMENT. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"
# BASIS, AND USDA ARS AND OREGON STATE UNIVERSITY HAVE NO OBLIGATIONS TO PROVIDE
# MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
#==============================================================================#
#' Calculate a dendrogram with bootstrap support using any distance applicable
#' to genind or genclone objects.
#'
#' @param x a [genind-class][genind], [genpop-class][genpop],
#'   [genclone-class][genclone], [genlight-class], [snpclone-class] or
#'   \link{matrix} object.
#'
#' @param strata a formula specifying the strata to be used to convert x to a
#'   genclone object if x is a genind object. Defaults to `NULL`. See details.
#'
#' @param tree a text string or function that can calculate a tree from a
#'   distance matrix. Defaults to "upgma". Note that you must load the package
#'   with the function for it to work.
#'
#' @param distance a character or function defining the distance to be applied
#'   to x. Defaults to [nei.dist()].
#'
#' @param sample An integer representing the number of bootstrap replicates
#'   Default is 100.
#'
#' @param cutoff An integer from 0 to 100 setting the cutoff value to return the
#'   bootstrap values on the nodes. Default is 0.
#'
#' @param showtree If `TRUE` (Default), a dendrogram will be plotted. If
#'   `FALSE`, nothing will be plotted.
#'
#' @param missing any method to be used by [missingno()]: "mean"
#'   (default), "zero", "loci", "genotype", or "ignore".
#'
#' @param mcutoff a value between 0 (default) and 1 defining the percentage of
#'   tolerable missing data if the `missing` parameter is set to "loci" or
#'   "genotype". This should only be set if the distance metric can handle
#'   missing data.
#'
#' @param quiet if `FALSE` (default), a progress bar will be printed to
#'   screen.
#'
#' @param root is the tree rooted? This is a parameter passed off to
#'   [ape::boot.phylo()]. If the `tree` parameter returns a
#'   rooted tree (like UPGMA), this should be `TRUE`, otherwise (like
#'   neighbor-joining), it should be false. When set to `NULL` (default),
#'   the tree is considered rooted if [ape::is.ultrametric()] is true.
#'
#' @param ... any parameters to be passed off to the distance method.
#'
#' @return an object of class [ape::phylo()].
#'
#' @details This function automates the process of bootstrapping genetic data to
#'   create a dendrogram with bootstrap support on the nodes. It will randomly
#'   sample with replacement the loci of a `gen` (genind/genpop) object or the columns of a
#'   numeric matrix, **assuming that all loci/columns are independent**. The
#'   process of randomly sampling `gen` objects with replacement is carried out
#'   through the use of an internal class called
#'   [bootgen-class]. This is necessary due to the fact that columns in
#'   the genind matrix are defined as alleles and are thus interrelated. This
#'   function will specifically bootstrap loci so that results are biologically
#'   relevant. With this function, the user can also define a custom distance to
#'   be performed on the genind or genclone object. If you have a data frame-like
#'   object where all of the columns are independent or pairs of columns are
#'   independent, then it may be simpler to use [ape::boot.phylo()] to calculate
#'   your bootstrap support values.
#'
#'   \subsection{the strata argument}{
#'   There is an argument called `strata`. This argument is useful for when
#'   you want to bootstrap by populations from a [adegenet::genind()]
#'   object. When you specify strata, the genind object will be converted to
#'   [adegenet::genpop()] with the specified strata.
#'   }
#'
#' @note [prevosti.dist()] and [diss.dist()] are exactly the
#'   same, but [diss.dist()] scales better for large numbers of
#'   individuals (n > 125) at the cost of required memory.
#'   \subsection{missing data}{
#'   Missing data is not allowed by many of the distances. Thus, one of
#'   the first steps of this function is to treat missing data by setting it to
#'   the average allele frequency in the data set. If you are using a distance
#'   that can handle missing data (Prevosti's distance), you can set
#'   `missing = "ignore"` to allow the distance function to handle any
#'   missing data. See [missingno()] for details on missing
#'   data.}
#'   \subsection{Bruvo's Distance}{
#'   While calculation of Bruvo's distance
#'   is possible with this function, it is optimized in the function
#'   [bruvo.boot()].}
#'
#' @seealso [nei.dist()] [edwards.dist()]
#'   [rogers.dist()] [reynolds.dist()]
#'   [prevosti.dist()] [diss.dist()]
#'   [bruvo.boot()] [ape::boot.phylo()]
#'   [adegenet::dist.genpop()] [dist()]
#'   [bootgen2genind()] [bootgen-class]
#'
#' @references
#' Kamvar ZN, Brooks JC and Grünwald NJ (2015) Novel R tools for analysis of
#' genome-wide population genetic data with emphasis on clonality.
#' Frontiers in Genetics 6:208. \doi{10.3389/fgene.2015.00208}
#' @export
#' @md
#' @keywords bootstrap
#' @aliases bootstrap
#' @examples
#'
#' data(nancycats)
#' nan9 <- popsub(nancycats, 9)
#'
#' set.seed(9999)
#' # Generate a tree using nei's distance
#' neinan <- aboot(nan9, dist = nei.dist)
#'
#' set.seed(9999)
#' # Generate a tree using custom distance
#' bindist <- function(x) dist(tab(x), method = "binary")
#' binnan <- aboot(nan9, dist = bindist)
#'
#' \dontrun{
#' # Distances from other packages.
#' #
#' # Sometimes, distance functions from other packages will have the constraint
#' # that the incoming data MUST be genind. Internally, aboot uses the
#' # bootgen class ( class?bootgen ) to shuffle loci, and will throw an error
#' # The function bootgen2genind helps fix that. Here's an example of a function
#' # that expects a genind class from above
#' bindist <- function(x){
#'   stopifnot(is.genind(x))
#'   dist(tab(x), method = "binary")
#' }
#' #
#' # Fails:
#' # aboot(nan9, dist = bindist)
#' ## Error: is.genind(x) is not TRUE
#' #
#' # Add bootgen2genind to get it working!
#' # Works:
#' aboot(nan9, dist = function(x) bootgen2genind(x) %>% bindist)
#'
#' # AFLP data
#' data(Aeut)
#'
#' # Nei's distance
#' anei <- aboot(Aeut, dist = nei.dist, sample = 1000, cutoff = 50)
#'
#' # Rogers' distance
#' arog <- aboot(Aeut, dist = rogers.dist, sample = 1000, cutoff = 50)
#'
#' # This can also be run on genpop objects
#' strata(Aeut) <- other(Aeut)$population_hierarchy[-1]
#' Aeut.gc <- as.genclone(Aeut)
#' setPop(Aeut.gc) <- ~Pop/Subpop
#' Aeut.pop <- genind2genpop(Aeut.gc)
#' set.seed(5000)
#' aboot(Aeut.pop, sample = 1000) # compare to Grunwald et al. 2006
#'
#' # You can also use the strata argument to convert to genpop inside the function.
#' set.seed(5000)
#' aboot(Aeut.gc, strata = ~Pop/Subpop, sample = 1000)
#'
#' # And genlight objects
#' # From glSim:
#' ## 1,000 non structured SNPs, 100 structured SNPs
#' x <- glSim(100, 1e3, n.snp.struc=100, ploid=2)
#' aboot(x, distance = bitwise.dist)
#'
#' # Utilizing other tree methods
#'
#' library("ape")
#'
#' aboot(Aeut.pop, tree = fastme.bal, sample = 1000)
#'
#' # Utilizing options in other tree methods
#'
#' myFastME <- function(x) fastme.bal(x, nni = TRUE, spr = FALSE, tbr = TRUE)
#' aboot(Aeut.pop, tree = myFastME, sample = 1000)
#'
#' }
#==============================================================================#
aboot <- function(
  x,
  strata = NULL,
  tree = "upgma",
  distance = "nei.dist",
  sample = 100,
  cutoff = 0,
  showtree = TRUE,
  missing = "mean",
  mcutoff = 0,
  quiet = FALSE,
  root = NULL,
  ...
) {
  if (!is.null(strata)) {
    if (!is.genind(x)) {
      warning(
        "Sorry, the strata argument can only be used with genind objects."
      )
    } else {
      x <- genind2genpop(x, pop = strata, quiet = TRUE, process.other = FALSE)
    }
  }
  if (is.matrix(x)) {
    if (is.null(rownames(x))) {
      rownames(x) <- .genlab("", nrow(x))
    }
    if (is.null(colnames(x))) {
      colnames(x) <- .genlab("L", ncol(x))
    }
    xboot <- x
  } else if (!is(x, "genlight") && x@type == "PA") {
    xboot <- x@tab
    colnames(xboot) <- locNames(x)
    rownames(xboot) <- if (is.genpop(x)) popNames(x) else indNames(x)
  } else if (is(x, "gen")) {
    if (is.genind(x)) {
      if (missing %in% c("loci", "geno", "ignore")) {
        x <- missingno(x, missing, quiet = quiet, cutoff = mcutoff)
        missing <- "asis"
      }
    }
    distname <- as.character(substitute(distance))
    if (length(distname) == 1 && distname %in% "diss.dist") {
      if (missing == "mean") {
        missing <- "asis"
        warning(
          "Sorry, missing = 'mean' is incompatible with diss.dist(), setting missing to 'asis'.",
          call. = FALSE
        )
      }
      xboot <- new("bootgen", x, na = missing, freq = FALSE)
    } else {
      xboot <- new("bootgen", x, na = missing, freq = TRUE)
    }
  } else if (is(x, "genlight")) {
    xboot <- x
    my_dists <- c(
      "diss.dist",
      "nei.dist",
      "prevosti.dist",
      "edwards.dist",
      "reynolds.dist",
      "rogers.dist",
      "provesti.dist"
    )
    wrong_dist <- any(as.character(substitute(distance)) %in% my_dists)
    if (wrong_dist) {
      warning(
        "Sorry, distance from genlight objects can only be calculated by bitwise.dist."
      )
      distance <- bitwise.dist
    }
  } else {
    stop("Sorry, x must be a genind, genpop, or genlight object.")
  }
  treefunk <- tree_generator(tree, distance, ...)
  xtree <- treefunk(xboot)
  if (any(xtree$edge.len < 0)) {
    xtree <- fix_negative_branch(xtree)
    warning(negative_branch_warning())
  }
  treechar <- paste(substitute(tree), collapse = "")
  if (is.null(root)) {
    root <- ape::is.ultrametric(xtree)
  }
  nodelabs <- boot.phylo(
    xtree,
    xboot,
    treefunk,
    B = sample,
    rooted = root,
    quiet = quiet
  )
  nodelabs <- (nodelabs / sample) * 100
  nodelabs <- ifelse(nodelabs >= cutoff, nodelabs, NA)
  if (!is.genpop(x)) {
    if (is.matrix(x)) {
      xtree$tip.label <- rownames(x)
    } else if (!is.null(indNames(x))) {
      xtree$tip.label <- indNames(x)
    }
  } else {
    xtree$tip.label <- popNames(x)
  }
  xtree$node.label <- nodelabs
  if (showtree) {
    poppr.plot.phylo(xtree, treechar, root)
  }
  return(xtree)
}

#==============================================================================#
#' Produce a table of diversity statistics
#'
#' Calculate diversity statistics on a matrix containing counts of multilocus
#' genotypes per population.
#'
#' @param z a table of integers representing counts of MLGs (columns) per
#'   population (rows)
#'
#' @param H logical whether or not to calculate Shannon's index
#' @param G logical whether or not to calculate Stoddart and Taylor's index (aka
#'   inverse Simpson's index).
#' @param lambda logical whether or not to calculate Simpson's index
#' @param E5 logical whether or not to calculate Evenness
#' @param ... any functions that can be calculated on a vector or matrix of
#'   genotype counts.
#'
#' @return a numeric matrix giving statistics (columns) for each population
#'   (rows).
#'
#' @details This function will calculate any diversity statistic for counts of
#'   multilocus genotypes per population. This does not count allelic diversity.
#'   The calculations of H, G, and lambda are all performed by
#'   [vegan::diversity()]. E5 is calculated as \deqn{E_{5} =
#'   \frac{(1/\lambda) - 1}{e^{H} - 1}}{(G - 1)/(exp(H) - 1)}.
#'
#' @export
#' @md
#' @seealso [diversity_boot()] [diversity_ci()]
#'   [poppr()]
#' @author Zhian N. Kamvar
#' @examples
#' library(poppr)
#' data(Pinf)
#' tab <- mlg.table(Pinf, plot = FALSE)
#' diversity_stats(tab)
#' \dontrun{
#' # Example using the poweRlaw package to calculate the negative slope of the
#' # Pareto distribution.
#'
#' library("poweRlaw")
#' power_law_beta <- function(x){
#'   xpow <- displ(x[x > 0])                 # Generate the distribution
#'   xpow$setPars(estimate_pars(xpow))       # Estimate the parameters
#'   xdat <- plot(xpow, draw = FALSE)        # Extract the data
#'   xlm <- lm(log(y) ~ log(x), data = xdat) # Run log-log linear model for slope
#'   return(-coef(xlm)[2])
#' }
#'
#' Beta <- function(x){
#'   x <- drop(as.matrix(x))
#'   if (length(dim(x)) > 1){
#'     res <- apply(x, 1, power_law_beta)
#'   } else {
#'     res <- power_law_beta(x)
#'   }
#'   return(res)
#' }
#'
#' diversity_stats(tab, B = Beta)
#' }
#==============================================================================#
diversity_stats <- function(
  z,
  H = TRUE,
  G = TRUE,
  lambda = TRUE,
  E5 = TRUE,
  ...
) {
  E.5 <- function(x) {
    H <- vegan::diversity(x)
    G <- vegan::diversity(x, "inv")
    (G - 1) / (exp(H) - 1)
  }
  BASIC_FUNS <- list(
    H = vegan::diversity,
    G = function(x) vegan::diversity(x, "inv"),
    lambda = function(x) vegan::diversity(x, "simp"),
    E.5 = E.5
  )
  FUNS <- c(BASIC_FUNS, list(...))
  STATS <- names(FUNS)
  STATS <- STATS[c(H, G, lambda, E5, rep(TRUE, length(list(...))))]

  # Indicator to determine if the input is from diversity_boot
  boot <- is.null(dim(z))

  nrows <- ifelse(boot, 1, nrow(z))
  if (boot) {
    dims <- list(NULL, Index = STATS)
  } else {
    dims <- list(Pop = rownames(z), Index = STATS)
  }

  mat <- matrix(nrow = nrows, ncol = length(STATS), dimnames = dims)
  for (i in STATS) {
    if (i == "E.5" && H && G) {
      mat[, i] <- (mat[, "G"] - 1) / (exp(mat[, "H"]) - 1)
    } else {
      mat[, i] <- FUNS[[i]](z)
    }
  }
  return(drop(mat))
}

#==============================================================================#
#' Perform a bootstrap analysis on diversity statistics
#'
#' This function is intended to perform bootstrap statistics on a matrix of
#' multilocus genotype counts in different populations. Results from this
#' function should be interpreted carefully as the default statistics are known
#' to have a downward bias. See the details for more information.
#'
#'
#' @param tab a table produced from the \pkg{poppr} function
#'   [mlg.table()]. MLGs in columns and populations in rows
#' @param n an integer > 0 specifying the number of bootstrap replicates to
#'   perform (corresponds to `R` in the function [boot::boot()].
#' @param n.boot an integer specifying the number of samples to be drawn in each
#'   bootstrap replicate. If `n.boot` < 2 (default), the number of samples
#'   drawn for each bootstrap replicate will be equal to the number of samples in
#'   the data set.
#' @param n.rare a sample size at which all resamplings should be performed.
#'   This should be no larger than the smallest sample size. Defaults to
#'   `NULL`, indicating that each population will be sampled at its own
#'   size.
#' @inheritParams diversity_stats
#' @param ... other parameters passed on to [boot::boot()] and
#'   [diversity_stats()].
#'
#' @return a list of objects of class "boot".
#' @seealso [diversity_stats()] for basic statistic calculation,
#'    [diversity_ci()] for confidence intervals and plotting, and
#'    [poppr()]. For bootstrap sampling:
#'    [stats::rmultinom()] [boot::boot()]
#'
#' @details
#'   Bootstrapping is performed in three ways:
#'   \itemize{
#'     \item if `n.rare` is a number greater than zero, then bootstrapping
#'     is performed by randomly sampling without replacement *n.rare*
#'     samples from the data.
#'
#'     \item if `n.boot` is greater than 1, bootstrapping is performed by
#'     sampling n.boot samples from a multinomial distribution weighted by the
#'     proportion of each MLG in the data.
#'
#'     \item if `n.boot` is less than 2, bootstrapping is performed by
#'     sampling N samples from a multinomial distribution weighted by the
#'     proportion of each MLG in the data.
#'   }
#'
#'   \subsection{Downward Bias}{
#'     When sampling with replacement, the diversity statistics here present a
#'     downward bias partially due to the small number of samples in the data.
#'     The result is that the mean of the bootstrapped samples will often be
#'     much lower than the observed value. Alternatively, you can increase the
#'     sample size of the bootstrap by increasing the size of `n.boot`. Both
#'     of these methods should be taken with caution in interpretation. There
#'     are several R packages freely available that will calculate and perform
#'     bootstrap estimates of Shannon and Simpson diversity metrics (eg.
#'     \pkg{entropart}, \pkg{entropy}, \pkg{simboot}, and
#'     \pkg{EntropyEstimation}. These packages also offer unbiased estimators of
#'     Shannon and Simpson diversity. Please take care when attempting to
#'     interpret the results of this function.
#'   }
#' @export
#' @md
#' @author Zhian N. Kamvar
#' @examples
#' library(poppr)
#' data(Pinf)
#' tab <- mlg.table(Pinf, plot = FALSE)
#' diversity_boot(tab, 10L)
#' \dontrun{
#' # This can be done in a parallel fashion (OSX uses "multicore", Windows uses "snow")
#' system.time(diversity_boot(tab, 10000L, parallel = "multicore", ncpus = 4L))
#' system.time(diversity_boot(tab, 10000L))
#' }
#' @importFrom boot boot boot.ci norm.ci
#==============================================================================#
diversity_boot <- function(
  tab,
  n,
  n.boot = 1L,
  n.rare = NULL,
  H = TRUE,
  G = TRUE,
  lambda = TRUE,
  E5 = TRUE,
  ...
) {
  if (!is.null(n.rare)) {
    FUN <- rare_sim_boot
    mle <- n.rare
  } else {
    FUN <- multinom_boot
    mle <- n.boot
  }
  res <- apply(
    tab,
    1,
    boot_per_pop,
    rg = FUN,
    n = n,
    mle = mle,
    H = H,
    G = G,
    lambda = lambda,
    E5 = E5,
    ...
  )
  return(res)
}

#==============================================================================#
#' Perform bootstrap statistics, calculate, and plot confidence intervals.
#'
#' This function is for calculating bootstrap statistics and their confidence
#' intervals. It is important to note that the calculation of confidence
#' intervals is not perfect (See Details). Please be cautious when interpreting
#' the results.
#'
#' @param tab a [genind()], [genclone()],
#'   [snpclone()], OR a matrix produced from
#'   [poppr::mlg.table()].
#' @param n an integer defining the number of bootstrap replicates (defaults to
#'   1000).
#' @param n.boot an integer specifying the number of samples to be drawn in each
#'   bootstrap replicate. If `n.boot` < 2 (default), the number of samples
#'   drawn for each bootstrap replicate will be equal to the number of samples
#'   in the data set. See Details.
#' @param ci the percent for confidence interval.
#' @param total argument to be passed on to [poppr::mlg.table()] if
#'   `tab` is a genind object.
#' @param rarefy if `TRUE`, bootstrapping will be performed on the smallest
#'   population size or the value of `n.rare`, whichever is larger.
#'   Defaults to `FALSE`, indicating that bootstrapping will be performed
#'   respective to each population size.
#' @param n.rare an integer specifying the smallest size at which to resample
#'   data. This is only used if `rarefy = TRUE`.
#' @param plot If `TRUE` (default), boxplots will be produced for each
#'   population, grouped by statistic. Colored dots will indicate the observed
#'   value.This plot can be retrieved by using `p <- last_plot()` from the
#'   \pkg{ggplot2} package.
#' @param raw if `TRUE` (default) a list containing three elements will be
#'   returned
#' @param center if `TRUE` (default), the confidence interval will be
#'   centered around the observed statistic. Otherwise, if `FALSE`, the
#'   confidence interval will be bias-corrected normal CI as reported from
#'   [boot::boot.ci()]
#' @param ... parameters to be passed on to [boot::boot()] and
#'   [diversity_stats()]
#'
#' @return \subsection{raw = TRUE}{
#'
#'    - **obs** a matrix with observed statistics in columns,
#'   populations in rows
#'    - **est** a matrix with estimated statistics in columns,
#'   populations in rows
#'    - **CI** an array of 3 dimensions giving the lower and upper
#'   bound, the index measured, and the population.
#'    - **boot** a list containing the output of
#'   [boot::boot()] for each population.
#'   }
#'
#' \subsection{raw = FALSE}{ a data frame with the statistic observations,
#' estimates, and confidence intervals in columns, and populations in rows. Note
#' that the confidence intervals are converted to characters and rounded to
#' three decimal places. }
#'
#' @details
#'   \subsection{Bootstrapping}{
#'   For details on the bootstrapping procedures, see
#'   [diversity_boot()]. Default bootstrapping is performed by
#'   sampling \strong{N} samples from a multinomial distribution weighted by the
#'   relative multilocus genotype abundance per population where \strong{N} is
#'   equal to the number of samples in the data set. If \strong{n.boot} > 2,
#'   then \strong{n.boot} samples are taken at each bootstrap replicate. When
#'   `rarefy = TRUE`, then samples are taken at the smallest population
#'   size without replacement. This will provide confidence intervals for all
#'   but the smallest population.
#'   }
#'   \subsection{Confidence intervals}{
#'   Confidence intervals are derived from the function
#'   [boot::norm.ci()]. This function will attempt to correct for bias
#'   between the observed value and the bootstrapped estimate. When `center
#'   = TRUE` (default), the confidence interval is calculated from the
#'   bootstrapped distribution and centered around the bias-corrected estimate
#'   as prescribed in Marcon (2012). This method can lead to undesirable
#'   properties, such as the confidence interval lying outside of the maximum
#'   possible value. For rarefaction, the confidence interval is simply
#'   determined by calculating the percentiles from the bootstrapped
#'   distribution. If you want to calculate your own confidence intervals, you
#'   can use the results of the permutations stored in the `$boot` element
#'   of the output.
#'   }
#'   \subsection{Rarefaction}{
#'   Rarefaction in the sense of this function is simply sampling a subset of
#'   the data at size **n.rare**. The estimates derived from this method
#'   have straightforward interpretations and allow you to compare diversity
#'   across populations since you are controlling for sample size.
#'   }
#'   \subsection{Plotting}{ Results are plotted as boxplots with point
#'   estimates. If there is no rarefaction applied, confidence intervals are
#'   displayed around the point estimates. The boxplots represent the actual
#'   values from the bootstrapping and will often appear below the estimates and
#'   confidence intervals.
#'   }
#' @note
#'   \subsection{Confidence interval calculation}{ Almost all of the statistics
#'   supplied here have a maximum when all genotypes are equally represented.
#'   This means that bootstrapping the samples will always be downwardly biased.
#'   In many cases, the confidence intervals from the bootstrapped distribution
#'   will fall outside of the observed statistic. The reported confidence
#'   intervals here are reported by assuming the variance of the bootstrapped
#'   distribution is the same as the variance around the observed statistic. As
#'   different statistics have different properties, there will not always be
#'   one clear method for calculating confidence intervals. A suggestion for
#'   correction in Shannon's index is to center the CI around the observed
#'   statistic (Marcon, 2012), but there are theoretical limitations to this.
#'   For details, see <https://stats.stackexchange.com/q/156235/49413>.
#'   }
#'
#'   \subsection{User-defined functions}{
#'   While it is possible to use custom functions with this, there are three
#'   important things to remember when using these functions:
#'
#'     1. The function must return a single value.
#'     2. The function must allow for both matrix and vector inputs
#'     3. The function name cannot match or partially match any arguments
#'     from [boot::boot()]
#'
#'   Anonymous functions are okay \cr(e.g. `function(x)
#'   vegan::rarefy(t(as.matrix(x)), 10)`).
#'   }
#' @export
#' @md
#' @seealso [diversity_boot()] [diversity_stats()]
#'   [poppr()] [boot::boot()] [boot::norm.ci()]
#'   [boot::boot.ci()]
#' @author Zhian N. Kamvar
#' @references
#' Marcon, E., Herault, B., Baraloto, C. and Lang, G. (2012). The Decomposition
#' of Shannon’s Entropy and a Confidence Interval for Beta Diversity.
#' *Oikos* 121(4): 516-522.
#'
#' @examples
#' library(poppr)
#' data(Pinf)
#' diversity_ci(Pinf, n = 100L)
#' \dontrun{
#' # With pretty results
#' diversity_ci(Pinf, n = 100L, raw = FALSE)
#'
#' # This can be done in a parallel fasion (OSX uses "multicore", Windows uses "snow")
#' system.time(diversity_ci(Pinf, 10000L, parallel = "multicore", ncpus = 4L))
#' system.time(diversity_ci(Pinf, 10000L))
#'
#' # We often get many requests for a clonal fraction statistic. As this is
#' # simply the number of observed MLGs over the number of samples, we
#' # recommended that people calculate it themselves. With this function, you
#' # can add it in:
#'
#' CF <- function(x){
#'  x <- drop(as.matrix(x))
#'  if (length(dim(x)) > 1){
#'    res <- rowSums(x > 0)/rowSums(x)
#'  } else {
#'    res <- sum(x > 0)/sum(x)
#'  }
#'  return(res)
#' }
#' # Show pretty results
#'
#' diversity_ci(Pinf, 1000L, CF = CF, center = TRUE, raw = FALSE)
#' diversity_ci(Pinf, 1000L, CF = CF, rarefy = TRUE, raw = FALSE)
#' }
#'
#==============================================================================#
diversity_ci <- function(
  tab,
  n = 1000,
  n.boot = 1L,
  ci = 95,
  total = TRUE,
  rarefy = FALSE,
  n.rare = 10,
  plot = TRUE,
  raw = TRUE,
  center = TRUE,
  ...
) {
  allowed_objects <- c("genind", "genclone", "snpclone")
  if (!is.matrix(tab) & inherits(tab, allowed_objects)) {
    tab <- mlg.table(tab, total = total, plot = FALSE)
  }
  rareval <- NULL
  if (rarefy) {
    rareval <- max(min(rowSums(tab)), n.rare)
    msg <- paste("\nSamples for rarefaction:", rareval)
    message(msg)
  } else {
    if (center == TRUE) {
      msg <- paste(
        "\nConfidence Intervals have been centered around observed",
        "statistic.\nPlease see ?diversity_ci for details.\n"
      )
      message(msg)
    } else if (n.boot > 2) {
      msg <- paste(
        "\nTake caution in interpreting the results.\nWhile sampling",
        n.boot,
        "samples will produce a confidence interval,\nit",
        "might be unclear what the correct interpretation should be.\n"
      )
      message(msg)
    }
  }
  res <- diversity_boot(tab, n, n.rare = rareval, n.boot = n.boot, ...)
  orig <- get_boot_stats(res)
  statnames <- colnames(orig)
  CI <- get_all_ci(
    res,
    ci = ci,
    index_names = statnames,
    center = if (rarefy) FALSE else center,
    btype = if (rarefy) "percent" else "normal",
    bci = if (rarefy) FALSE else TRUE
  )
  est <- get_boot_se(res, "mean")
  out <- list(obs = orig, est = est, CI = CI, boot = res)
  if (plot) {
    boot_plot(res = res, orig = orig, CI = if (rarefy) NULL else CI)
  }
  out <- if (raw) out else do.call("pretty_info", out)
  return(out)
}