An Introduction to iNEXT.beta3D via Examples

The package iNEXT.beta3D (iNterpolation and EXTrapolation with beta diversity for three dimensions of biodiversity) is a sequel to iNEXT. The three dimensions (3D) of biodiversity include taxonomic diversity (TD), phylogenetic diversity (PD) and functional diversity (FD). This document provides an introduction demonstrating how to run iNEXT.beta3D. An online version iNEXT.beta3D Online is also available for users without an R background.

A unified framework based on Hill numbers (for TD) and their generalizations, Hill-Chao numbers (for PD and FD), is adopted in the package. TD quantifies the effective number of species, meanPD (PD divided by reference time) quantifies the effective number of lineages, and FD quantifies the effective number of virtual functional groups (or functional “species”). Thus, TD, meanPD, and FD are all expressed in comparable units of species/lineage equivalents, allowing meaningful comparisons across the three dimensions of biodiversity; see Chao et al. (2021) for a review of the unified framework.

For each of the three dimensions, iNEXT.beta3D focuses on the multiplicative decomposition of diversity (beta = gamma/alpha, expressed in units of effective number of assemblages) for orders q = 0, 1 and 2, based on sampling data. The resulting 3D beta diversity quantifies the extent of among-assemblage differentiation, reflecting changes in species, lineages, or functional-groups composition and their abundances across assemblages.

iNEXT.beta3D provides standardized 3D estimates at a common sample size (for alpha and gamma diversity only, because sample-size-based beta diversity is not generally valid) or at a common sample coverage (for alpha, beta and gamma diversity). iNEXT.beta3D also provides standardized estimates for four classes of coverage-based dissimilarity measures.

Based on rarefaction and extrapolation (R/E) theory for Hill-numbers-based TD of orders q = 0, 1 and 2, Chao et al. (2023b) developed the corresponding R/E theory for taxonomic beta diversity, with applications to real-world spatial, temporal and spatio-temporal data. An application to Gentry’s global forest data together with a concise description of the theory, is provided in Chao et al. (2023a). The extension to phylogenetic and functional beta diversity follows a generally parallel framework. iNEXT.beta3D features two types of R/E sampling curves:

Sample-size-based (or size-based) R/E sampling curves: These curves plot standardized 3D gamma and alpha diversity as functions of sample size. Note that the size-based beta diversity is not generally a statistically valid measure (Chao et al. 2023b); therefore, the corresponding beta diversity sampling curve is not provided.
Sample-coverage-based (or coverage-based) R/E sampling curves: These curves plot standardized 3D gamma, alpha, and beta diversity, as well as four classes of dissimilarity measures, as functions of sample coverage, an objective measure of sample completeness.

Note that sufficient data are required to run iNEXT.beta3D. If a dataset contains only a few species, along with their abundances, phylogenies, or traits, it may not contain sufficient information for reliable iNEXT.beta3D analyses.

HOW TO CITE iNEXT.beta3D

If you publish work based on results obtained from iNEXT.beta3D, please cite at least one of the following methodological papers (2023a, b) and also the iNEXT.beta3D package:

Chao, A., Chiu, C.-H., Hu, K.-H., and Zelený, D. (2023a). Revisiting Alwyn H. Gentry’s forest transect data: a statistical sampling-model-based approach. Japanese Journal of Statistics and Data Science, 6, 861-884. (https://doi.org/10.1007/s42081-023-00214-1)
Chao, A., Thorn, S., Chiu, C.-H., Moyes, F., Hu, K.-H., Chazdon, R. L., Wu, J., Magnago, L. F. S., Dornelas, M., Zelený, D., Colwell, R. K., and Magurran, A. E. (2023b). Rarefaction and extrapolation with beta diversity under a framework of Hill numbers: the iNEXT.beta3D standardization. Ecological Monographs e1588.(https://doi.org/10.1002/ecm.1588)
Chao, A., Lee, Y.-T., Lin, K.-L., and Chuang, P.-Y. (2026). The iNEXT.beta3D package: interpolation and extrapolation with beta diversity for three dimensions of biodiversity. R package version 1.1.0.

SOFTWARE NEEDED TO RUN iNEXT.beta3D IN R

Required: R
Suggested: RStudio IDE

HOW TO RUN iNEXT.beta3D:

The iNEXT.beta3D package is available from CRAN and can be downloaded from Anne Chao’s GitHub iNEXT.beta3D_github using the following commands. For a first-time installation, additional visualization extension package (ggplot2 from CRAN) and relevant package (iNEXT.3D from CRAN) must be installed and loaded.


## install iNEXT.beta3D package from CRAN
install.packages("iNEXT.beta3D")

## install the latest version from GitHub
install.packages('devtools')
library(devtools)
install_github('AnneChao/iNEXT.beta3D')

## import packages
library(iNEXT.beta3D)

There are three main functions in this package:

iNEXTbeta3D(): Computes standardized 3D estimates at common sample sizes (for gamma and alpha diversity) or common sample coverage levels (for gamma, alpha, and beta diversity), as well as four classes of coverage-based dissimilarity measures. Users can obtain estimates at default or user-specified sample sizes or coverage levels.
ggiNEXTbeta3D(): Visualizes the output from iNEXTbeta3D().
DataInfobeta3D(): Provides basic information about (1) the reference sample from each individual assemblage, (2) the gamma reference sample from the pooled assemblage, and (3) the alpha reference sample from the joint assemblage.

DATA INPUT FORMAT

To assess taxonomic, phylogenetic, and functional beta diversity among assemblages, information on species identity (or a unique identification code) and assemblage affiliation is required. For each dataset, species should be arranged in rows, with species names (or identification codes) used as row names. The columns represent assemblages for abundance data or sampling units for incidence data, as described below. Two types of species abundance/incidence data are supported:

Individual-based abundance data (datatype = "abundance"): The input data for a single dataset with N assemblages should be arranged as a species-by-assemblage abundance matrix/data.frame, with species in rows and assemblages in columns. Species names (or unique identification codes) should be used as row names, and assemblage names as column names. Users can input multiple datasets collected from different localities, regions, plots, time periods, or other sampling contexts. The input for multiple datasets is a list of matrices/data.frames, with each element representing one dataset. Different datasets may contain different numbers of assemblages. For each dataset, iNEXTbeta3D() computes gamma, alpha, and beta diversity, as well as dissimilarity among all N assemblages and other relevant statistics. Users can also set by_pair = TRUE to obtain the corresponding results for every pair of assemblages.
Sampling-unit-based incidence raw data (datatype = "incidence_raw"): For a dataset with N assemblages, the input data consist of a list of N matrices/data.frames, with each matrix representing a species-by-sampling-unit incidence matrix for one assemblage. Species are arranged in rows and sampling units in columns. Each entry must be either 0 (non-detection) or 1 (detection). Users can input multiple datasets, with each dataset represented by a list of species-by-sampling-unit incidence matrices/data.frames. The number of sampling units may vary among datasets. However, within each dataset, all assemblages must have the same number of sampling units. For each dataset, iNEXTbeta3D() computes gamma, alpha, and beta diversity, as well as dissimilarity among all N assemblages and other relevant statistics, based on incidence frequency counts obtained from the sampling units. Users can also set by_pair = TRUE to obtain the corresponding results for every pair of assemblages.

Species abundance data format

We use the tree species abundance data collected from two rainforest fragments/localities (“Marim” and “Rebio2”) in Brazil to assess beta diversity between the Edge and Interior habitats (assemblages) within each fragment; see Chao et al. (2023b) for details of the analysis. The demo data differ slightly from those analyzed in Chao et al. (2023b), as seven species were removed from the original pooled data because phylogenetic information was unavailable for these species. The data, named Brazil_rainforests, consist of a list of two species-by-assemblage abundance matrices, one for each fragment (“Marim” and “Rebio2”). Each matrix contains two assemblages (columns), “Edge” and “Interior”, with species in rows. Run the following code to view the data (only the first 15 rows of each matrix are shown).


data(Brazil_rainforests)
Brazil_rainforests

#> $Marim
#>                            Edge Interior
#> Acosmium_lentiscifolium       1        0
#> Allophylus_petiolulatus       5        0
#> Alseis_involuta               2        0
#> Ampelocera_glabra             1        0
#> Andira_legalis                0        1
#> Andira_ormosioides            0        1
#> Apuleia_leiocarpa             1        0
#> Aspidosperma_illustre         0        3
#> Astrocaryum_aculeatissimum    1        0
#> Astronium_concinnum           4        1
#> Barnebydendron_riedelii       0        2
#> Bauhinia_forficata            1        0
#> Brosimum_glaucum              4        0
#> Calyptranthes_lucida          0        4
#> Campomanesia_lineatifolia     1        0
#> 
#> $Rebio2
#>                             Edge Interior
#> Albizia_polycephala            1        0
#> Allophylus_petiolulatus        3        3
#> Alseis_involuta                1        0
#> Amaioua_intermedia             0        1
#> Ampelocera_glabra              0        3
#> Anaxagorea_silvatica           0        6
#> Annona_dolabripetala           1        0
#> Aspidosperma_cylindrocarpon    2        0
#> Astrocaryum_aculeatissimum     7        1
#> Astronium_concinnum           12        1
#> Astronium_graveolens          13        1
#> Beilschmiedia_linharensis      1        0
#> Brosimum_glaucum               2        2
#> Brosimum_sp1                   0        1
#> Calyptranthes_lucida           2        1

Species incidence raw data format

We use tree species data collected from two second-growth rainforests, Cuatro Rios (CR) and Juan Enriquez (JE), in Costa Rica, as demo data to assess temporal beta diversity among three years (2005, 2011, and 2017) within each forest, with each year treated as an assemblage. In each forest, data were collected from a 1-ha (50 m × 200 m) forest plot.

Because individuals of some tree species may exhibit intraspecific spatial aggregation within a 1-ha plot, individual trees may not be appropriately treated as independent sampling units. In such cases, it is statistically preferable to convert species abundance data within each forest into occurrence or incidence (detection/non-detection) data for subplots or quadrats. See Chao et al. (2023b) for details of the analysis.

Each 1-ha forest was divided into 100 subplots (0.01 ha each), and only species incidence (detection/non-detection) records within each subplot were used to compute the incidence frequency of each species (i.e., the number of subplots in which the species was detected). By treating species incidence frequencies across subplots as proxies for species abundances, iNEXT.beta3D standardization can be adapted to accommodate spatially aggregated data and reduce the effects of intraspecific aggregation.

The data, named Second_growth_forests, consist of a list containing two forests ("CR_2005_2011_2017" and "JE_2005_2011_2017"). Each forest is represented by a list of three species-by-sampling-unit incidence matrices corresponding to 2005, 2011, and 2017. Each matrix contains species in rows and 100 subplots (sampling units) in columns.

Run the following code to view the incidence raw data (only the first 10 rows and 6 columns of each matrix are shown).


data(Second_growth_forests)
Second_growth_forests

#> $CR_2005_2011_2017
#> $CR_2005_2011_2017$Year_2005
#>        Subplot_1 Subplot_2 Subplot_3 Subplot_4 Subplot_5 Subplot_6
#> Abaade         0         0         0         0         0         0
#> Alcflo         0         0         0         0         0         0
#> Alclat         0         1         0         0         0         0
#> Aliatl         0         0         0         0         0         0
#> Ampmac         0         0         0         0         0         0
#> Anacra         0         1         0         0         0         1
#> Annama         0         1         0         0         0         0
#> Annpap         0         0         0         0         0         0
#> Apemem         0         0         0         0         0         0
#> Ardfim         0         0         0         0         0         0
#> 
#> $CR_2005_2011_2017$Year_2011
#>        Subplot_1 Subplot_2 Subplot_3 Subplot_4 Subplot_5 Subplot_6
#> Abaade         0         0         0         0         0         0
#> Alcflo         0         0         0         0         0         0
#> Alclat         0         1         0         0         0         0
#> Aliatl         0         0         0         0         0         0
#> Ampmac         0         0         0         0         0         0
#> Anacra         0         1         0         0         0         1
#> Annama         0         0         0         0         0         0
#> Annpap         0         0         0         0         0         0
#> Apemem         0         0         0         0         0         0
#> Ardfim         0         0         0         0         0         0
#> 
#> $CR_2005_2011_2017$Year_2017
#>        Subplot_1 Subplot_2 Subplot_3 Subplot_4 Subplot_5 Subplot_6
#> Abaade         0         0         0         0         0         0
#> Alcflo         0         0         0         0         0         0
#> Alclat         0         1         0         0         0         0
#> Aliatl         0         0         0         0         0         0
#> Ampmac         0         0         0         0         0         0
#> Anacra         0         1         1         0         1         1
#> Annama         0         0         0         0         0         0
#> Annpap         0         0         0         0         0         0
#> Apemem         0         0         0         0         0         0
#> Ardfim         0         0         0         0         0         0
#> 
#> 
#> $JE_2005_2011_2017
#> $JE_2005_2011_2017$Year_2005
#>        Subplot_1 Subplot_2 Subplot_3 Subplot_4 Subplot_5 Subplot_6
#> Alccos         0         0         0         0         0         0
#> Alcflo         0         0         0         0         0         0
#> Alclat         0         0         0         0         0         0
#> Annpap         0         0         0         0         0         0
#> Apemem         0         0         0         0         0         0
#> Astcon         0         0         0         0         0         0
#> Bacgas         0         0         0         0         0         0
#> Brogui         0         0         0         0         0         0
#> Brolac         0         0         0         0         0         0
#> Byrcra         0         0         0         0         1         0
#> 
#> $JE_2005_2011_2017$Year_2011
#>        Subplot_1 Subplot_2 Subplot_3 Subplot_4 Subplot_5 Subplot_6
#> Alccos         0         0         0         0         0         0
#> Alcflo         0         0         0         0         0         0
#> Alclat         0         0         0         0         0         0
#> Annpap         0         0         0         0         0         0
#> Apemem         0         0         0         0         0         0
#> Astcon         0         0         0         0         0         0
#> Bacgas         0         0         0         0         0         0
#> Brogui         0         0         0         0         0         0
#> Brolac         0         0         0         0         0         0
#> Byrcra         0         0         0         0         1         0
#> 
#> $JE_2005_2011_2017$Year_2017
#>        Subplot_1 Subplot_2 Subplot_3 Subplot_4 Subplot_5 Subplot_6
#> Alccos         0         0         0         0         0         0
#> Alcflo         0         0         0         0         0         0
#> Alclat         0         0         0         0         0         0
#> Annpap         0         0         0         0         0         0
#> Apemem         0         0         0         0         0         0
#> Astcon         0         0         0         0         0         0
#> Bacgas         0         0         0         0         0         0
#> Brogui         0         0         0         0         0         0
#> Brolac         0         0         0         0         0         0
#> Byrcra         0         0         0         0         0         0

Phylogenetic tree format for PD

To perform phylogenetic diversity (PD) analysis, a phylogenetic tree (in Newick format) spanning all species observed across the datasets is required. For example, the phylogenetic tree spanning all species observed in the Marim and Rebio2 fragments is stored in the R dataset Brazil_tree. This dataset is used to demonstrate the computation of PD. A partial list of tip labels and node labels is shown below.


data(Brazil_tree)
Brazil_tree
#> 
#> Phylogenetic tree with 185 tips and 117 internal nodes.
#> 
#> Tip labels:
#>   Carpotroche_brasiliensis, Casearia_ulmifolia, Casearia_sp2, Casearia_oblongifolia, Casearia_commersoniana, Rinorea_bahiensis, ...
#> Node labels:
#>   magnoliales_to_asterales, poales_to_asterales, , , , , ...
#> 
#> Rooted; includes branch length(s).

Species pairwise distance matrix format for FD

To perform functional diversity (FD) analysis, a pairwise distance matrix among all species observed across the datasets is required. The distance matrix can be computed from species traits using Gower distance. For example, the pairwise distance matrix for all species observed in the Marim and Rebio2 fragments is available in the R dataset "Brazil_distM". This dataset is used to demonstrate the computation of FD. The first three rows and three columns of the distance matrix are shown below.


data(Brazil_distM)
Brazil_distM

#>                          Carpotroche_brasiliensis Astronium_concinnum Astronium_graveolens
#> Carpotroche_brasiliensis                    0.000               0.522                0.522
#> Astronium_concinnum                         0.522               0.000                0.000
#> Astronium_graveolens                        0.522               0.000                0.000

MAIN FUNCTION: iNEXTbeta3D()

We first describe the main function iNEXTbeta3D() with default arguments:


iNEXTbeta3D(data, diversity = "TD", q = c(0, 1, 2), datatype = "abundance",
            base = "coverage", level = NULL, nboot = 10, conf = 0.95,
            PDtree = NULL, PDreftime = NULL, PDtype = 'AUC',  PDcut_number = 30,
            FDdistM = NULL, FDtype = "AUC", FDtau = NULL, FDcut_number = 30,
            by_pair = FALSE)

The arguments of this function are briefly described below and explained in more detail through illustrative examples in subsequent sections. By default (with base = “coverage”), the function computes coverage-based standardized 3D gamma, alpha, and beta diversity, as well as four dissimilarity indices, using coverage levels from 0.5 to 1. For q = 0, estimates corresponding to sample sizes greater than twice the reference sample size are excluded. When base = “size”, the function computes size-based standardized 3D gamma and alpha diversity estimates for sample sizes up to twice the reference sample size in each dataset. The function also computes standardized 3D estimates at user-specified sample sizes or coverage levels.

Argument	Description
`data`	For `datatype = “abundance”`, data for a single dataset can be provided as a species-by-assemblage abundance `matrix/data.frame`, with species in rows and assemblages in columns. Data for multiple datasets should be provided as a `list` of such `matrices/data.frames`, with each element representing one dataset. For `datatype = “incidence_raw”`, data for a single dataset should be provided as a `list` of species-by-sampling-unit incidence `matrices/data.frames`, with each element representing one assemblage. Data for multiple datasets should be provided as a list of such lists, with each outer-list element representing one dataset.
`diversity`	Selection of the diversity type: `diversity = “TD”` for Taxonomic diversity, `diversity = “PD”` for Phylogenetic diversity, and `diversity = “FD”` for Functional diversity.
`q`	A numerical vector specifying the diversity orders. Default is `c(0, 1, 2)`.
`datatype`	Type of input data: individual-based abundance data (`datatype = “abundance”`) or sampling-unit-based incidence data (`datatype = “incidence_raw”`). For incidence raw data, each species-by-sampling-unit matrix contains only 0 (non-detection) or 1 (detection).
`base`	Standardization base: coverage-based rarefaction and extrapolation for gamma, alpha, beta diversity, as well as four classes of dissimilarity indices (`base = “coverage”`), or size-based rarefaction and extrapolation for gamma and alpha diversity (`base = “size”`). Default is `base = “coverage”`.
`level`	A numeric vector specifying the target sample coverage levels (when `base = “coverage”`) or sample sizes (when `base = “size”`) used to compute standardized diversity and dissimilarity estimates. For coverage-based standardization, the specified coverage levels must be between 0 and 1. Setting `level = 1` requests estimates at complete coverage (i.e., asymptotic estimates), when supported by the estimator. By default, when `base = “coverage”`, the function initially uses coverage levels ranging from 0.5 to 1 in increments of 0.025, together with the estimated gamma and alpha coverage levels at the reference sample size and twice the reference sample size. For q = 0, estimates corresponding to sample sizes greater than twice the reference sample size are excluded. Thus, the extrapolation limits correspond to SC(2n, gamma) for gamma diversity and SC(2n, alpha) for alpha and beta diversity, as well as dissimilarity measures. For q = 1 and 2, estimates can be extended to complete coverage when supported by the data. The reference coverage for beta diversity and dissimilarity is defined as that for alpha diversity. When `base = “size”`, the function computes standardized 3D gamma and alpha diversity estimates at sample-size knots spanning rarefaction, the reference sample size, and extrapolation up to twice the reference sample size.
`nboot`	A non-negative integer specifying the number of bootstrap replications used to assess sampling uncertainty and construct confidence intervals. Bootstrap replications can be time-consuming. Set `nboot = 0` to skip the bootstrap procedures. Default is `nboot = 10`. For more accurate uncertainty estimates, a larger number of replications (e.g., 100 or 200) is recommended.
`conf`	A number between 0 and 1 specifying the confidence level. Default is `conf = 0.95`.
`PDtree`	(required when `diversity = “PD”`), a phylogenetic tree (an R `phylo` object) spanning all observed species across the input datasets. The tree can be imported from a file in Newick format.
`PDreftime`	(argument only for `diversity = “PD”` and `PDtype = “meanPD”`), a numerical value specifying a fixed time reference point for phylogenetic diversity. Default is `PDreftime = NULL`, in which case the age of the root of the phylogenetic tree spanned by the pooled observed species is used.
`PDtype`	(argument only for `diversity = “PD”`), specifies the phylogenetic diversity measure: `PDtype = “meanPD”` for the effective number of equally divergent lineages at a specified reference time, or `PDtype = “AUC”` for an overall measure obtained by integrating and averaging mean phylogenetic diversity over all reference times from zero to the root age of the phylogenetic tree. Default is `PDtype = “AUC”`.
`PDcut_number`	(argument only for `diversity = “PD”` and `PDtype = “AUC”`), an integer specifying the number of equally spaced integration points used to compute the AUC over the interval from zero to the root age of the phylogenetic tree. Default is `PDcut_number = 30`, and the minimum is 3.
`FDdistM`	(required when `diversity = “FD”`), a species pairwise distance matrix covering all species observed across the input datasets.
`FDtype`	(argument only for `diversity = “FD”`), specifies the functional diversity measure: `FDtype = “tau_value”` for FD at a specified threshold of functional distinctiveness, or `FDtype = “AUC”` for an overall FD measure obtained by integrating the tau-profile over threshold values from 0 to 1. Default is `FDtype = “AUC”`.
`FDtau`	(argument only for `diversity = “FD”` and `FDtype = “tau_value”`), a numeric value between 0 and 1 specifying the threshold of functional distinctiveness used to compute FD. If `FDtau = NULL` (default), the threshold is set to the mean pairwise functional distance between two individuals randomly selected from the pooled data across all datasets (i.e., Rao’s quadratic entropy).
`FDcut_number`	(argument only for `diversity = “FD”` and `FDtype = “AUC”`), an integer specifying the number of equally spaced integration points used to compute the AUC over the interval from 0 to 1. Default is `FDcut_number = 30`, and the minimum is 3.
`by_pair`	A logical variable specifying whether diversity decomposition is performed for all assemblages as a whole or for every pair of assemblages. If `by_pair = FALSE`, gamma, alpha, and beta diversity, as well as dissimilarity, are computed for all assemblages within each dataset considered as a whole. If `by_pair = TRUE`, the corresponding estimates are computed for every pair of assemblages. Default is `by_pair = FALSE`.

This function returns an iNEXTbeta3D object that can be visualized using ggiNEXTbeta3D(), as described below. When by_pair = FALSE, ggiNEXTbeta3D() returns a single ggplot object. When by_pair = TRUE, it returns a named list of ggplot objects, with one plot for each dataset. Within each plot, different pairs of assemblages are distinguished by color.

Output of the main function iNEXTbeta3D()

By default (with base = "coverage"), for each of the three diversity dimensions (TD, PD, and FD), iNEXTbeta3D() returns an "iNEXTbeta3D" object containing seven data frames for each dataset:

gamma (standardized gamma diversity)
alpha (standardized alpha diversity)
beta (standardized beta diversity)
1-C (standardized Sørensen-type non-overlap index)
1-U (standardized Jaccard-type non-overlap index)
1-V (standardized Sørensen-type turnover index)
1-S (standardized Jaccard-type turnover index)

For PDtype = "meanPD", the gamma output includes Gamma_PD and Gamma_meanPD, whereas the alpha output includes Alpha_PD and Alpha_meanPD. These quantities satisfy Gamma_PD = Gamma_meanPD × Reftime and Alpha_PD = Alpha_meanPD × Reftime. The PD-based and meanPD-based beta diversity and dissimilarity estimates are identical after multiplicative decomposition; therefore, only one set of estimates is reported.

When base = "size", for each of the three diversity dimensions (TD, PD, and FD), iNEXTbeta3D() returns an "iNEXTbeta3D" object containing two data frames for each dataset:

gamma (size-based standardized gamma diversity)
alpha (size-based standardized alpha diversity)

Size-based beta diversity and dissimilarity indices are not statistically valid measures and are therefore not provided.

GRAPHIC DISPLAYS: FUNCTION ggiNEXTbeta3D()

The function ggiNEXTbeta3D() and its default arguments are introduced below. It supports output generated by iNEXTbeta3D() with both by_pair = FALSE and by_pair = TRUE.


ggiNEXTbeta3D(output, type = "B")

Argument	Description
`output`	An object returned by `iNEXTbeta3D()`.
`type`	(argument only for `base = “coverage”`), `type = “B”` generates R/E sampling curves for gamma, alpha, and beta diversity, whereas `type = “D”` generates R/E sampling curves for four dissimilarity indices. The `type` argument is not required when plotting size-based R/E sampling curves for gamma and alpha diversity.

type = “B” generates R/E sampling curves for gamma, alpha, and beta diversity, whereas type = “D” generates R/E sampling curves for four dissimilarity indices. The type argument is not required when plotting size-based R/E sampling curves for gamma and alpha diversity.

generates R/E sampling curves for gamma, alpha, and beta diversity, whereas generates R/E sampling curves for four dissimilarity indices. The argument is not required when plotting size-based R/E sampling curves for gamma and alpha diversity.

The ggiNEXTbeta3D() function is a wrapper around the ggplot2 package that creates rarefaction and extrapolation (R/E) sampling curves with a single line of code. When by_pair = FALSE, it returns a single "ggplot" object that can be further customized using ggplot2 functions. When by_pair = TRUE, it returns a named list of "ggplot" objects, with one plot for each dataset.

TAXONOMIC DIVERSITY (TD): RAREFACTION/EXTRAPOLATION VIA EXAMPLES

EXAMPLE 1: Abundance data with default sample sizes or coverage values

First, we use the iNEXTbeta3D() function with the Brazil_rainforests abundance data to compute coverage-based standardized taxonomic gamma, alpha, and beta diversity, as well as four dissimilarity indices, by running the following code:


## Coverage-based R/E Analysis with taxonomic diversity for abundance data
data(Brazil_rainforests)

output_TDc_abun = iNEXTbeta3D(data = Brazil_rainforests, diversity = 'TD', 
                              datatype = 'abundance', base = "coverage", nboot = 10)
output_TDc_abun

The output contains seven data frames: gamma, alpha, beta, 1-C, 1-U, 1-V, and 1-S. Each data frame includes the dataset name (Dataset), the diversity order (Order.q), and the target standardized coverage value (SC). When by_pair = TRUE, the output also includes the assemblage pair (Pair). The corresponding sample size (Size for abundance data or mT for incidence data) is reported for gamma and alpha diversity, but not for beta diversity or dissimilarity indices. The estimated diversity or dissimilarity is reported in the corresponding Gamma, Alpha, Beta, or Dissimilarity column. The Method column indicates Rarefaction, Observed, or Extrapolation, depending on whether the target coverage is less than, equal to, or greater than the coverage of the reference sample. For beta diversity and dissimilarity measures, the reference coverage is defined as that for alpha diversity. The output also includes the standard error (s.e.) and the bootstrap lower and upper confidence limits (LCL and UCL). These estimates and confidence intervals are used to construct rarefaction and extrapolation (R/E) sampling curves. For PDtype = "meanPD", gamma diversity is reported as Gamma_PD and Gamma_meanPD, whereas alpha diversity is reported as Alpha_PD and Alpha_meanPD, with Gamma_PD = Gamma_meanPD × Reftime and Alpha_PD = Alpha_meanPD × Reftime.

By default, for coverage-based standardization, the function initially considers standardized coverage levels ranging from 0.5 to 1 in increments of 0.025. For q = 0, estimates corresponding to sample sizes greater than twice the reference sample size are excluded. Thus, the extrapolation limits correspond to SC(2n, gamma) for gamma diversity and SC(2n, alpha) for alpha and beta diversity, as well as dissimilarity measures. For q = 1 and 2, estimates can be extended to complete coverage (SC = 1) when supported by the data. In addition, the following four coverage values are included in the initial set of coverage levels: SC(n, alpha), SC(2n, alpha), SC(n, gamma), and SC(2n, gamma). The corresponding results are retained when they fall within the applicable estimation range. Here, SC(n, alpha) denotes the estimated coverage for the alpha reference sample size n, whereas SC(2n, alpha) denotes the estimated coverage for an extrapolated alpha sample of size 2n in the joint assemblage. These values are reported as SC(n) and SC(2n) for "Joint assemblage (for alpha)" in the Assemblage column of the output from DataInfobeta3D(). Similarly, SC(n, gamma) and SC(2n, gamma) denote the corresponding coverage values for the gamma reference sample and its extrapolated sample of size 2n in the pooled assemblage. These values are reported as SC(n) and SC(2n) for "Pooled assemblage (for gamma)" in the output from DataInfobeta3D(). For beta diversity and dissimilarity measures, the observed sample coverage and extrapolation limit are defined as those for alpha diversity. For incidence data, the corresponding coverage values are denoted by SC(T, alpha), SC(2T, alpha), SC(T, gamma), and SC(2T, gamma), respectively.

Because the default standardized coverage range starts at 0.5, the default setting (level = NULL) is not applicable when the observed sample coverage of either the alpha or gamma reference sample is below 50%. In such cases, users should specify appropriate coverage values through the level argument. Extrapolation beyond twice the reference sample size for q = 0 is not recommended because the resulting estimates may be subject to increased bias. Below, we show the output for taxonomic beta diversity between the “Edge” and “Interior” habitats in the “Marim” fragment.

#>    Dataset Order.q    SC Beta                Method  s.e.   LCL  UCL
#> 1    Marim       0 0.500 1.11           Rarefaction 0.061 0.990 1.23
#> 2    Marim       0 0.525 1.11           Rarefaction 0.064 0.983 1.23
#> 3    Marim       0 0.550 1.10           Rarefaction 0.066 0.976 1.23
#> 4    Marim       0 0.575 1.10           Rarefaction 0.069 0.967 1.24
#> 5    Marim       0 0.600 1.10           Rarefaction 0.072 0.958 1.24
#> 6    Marim       0 0.625 1.09           Rarefaction 0.075 0.949 1.24
#> 7    Marim       0 0.650 1.09           Rarefaction 0.078 0.939 1.25
#> 8    Marim       0 0.675 1.09           Rarefaction 0.082 0.929 1.25
#> 9    Marim       0 0.696 1.09 Observed_SC(n, alpha) 0.085 0.920 1.25
#> 10   Marim       0 0.700 1.09         Extrapolation 0.085 0.919 1.25
#> 11   Marim       0 0.725 1.08         Extrapolation 0.088 0.912 1.26
#> 12   Marim       0 0.750 1.08         Extrapolation 0.089 0.909 1.26
#> 13   Marim       0 0.775 1.08         Extrapolation 0.089 0.911 1.26
#> 14   Marim       0 0.800 1.09         Extrapolation 0.087 0.917 1.26
#> 15   Marim       0 0.825 1.09         Extrapolation 0.085 0.924 1.26
#> 16   Marim       0 0.850 1.09         Extrapolation 0.082 0.932 1.25
#> 17   Marim       0 0.855 1.09 Observed_SC(n, gamma) 0.082 0.933 1.25
#> 18   Marim       0 0.875 1.09         Extrapolation 0.080 0.938 1.25
#> 19   Marim       0 0.876 1.09  Extrap_SC(2n, alpha) 0.080 0.938 1.25
#> 20   Marim       1 0.500 1.11           Rarefaction 0.056 1.000 1.22
#> 21   Marim       1 0.525 1.11           Rarefaction 0.057 0.995 1.22
#> 22   Marim       1 0.550 1.11           Rarefaction 0.059 0.991 1.22
#> 23   Marim       1 0.575 1.10           Rarefaction 0.060 0.986 1.22
#> 24   Marim       1 0.600 1.10           Rarefaction 0.061 0.981 1.22
#> 25   Marim       1 0.625 1.10           Rarefaction 0.063 0.975 1.22
#> 26   Marim       1 0.650 1.10           Rarefaction 0.065 0.970 1.22
#> 27   Marim       1 0.675 1.09           Rarefaction 0.066 0.965 1.23
#> 28   Marim       1 0.696 1.09 Observed_SC(n, alpha) 0.068 0.961 1.23
#> 29   Marim       1 0.700 1.09         Extrapolation 0.068 0.961 1.23
#> 30   Marim       1 0.725 1.09         Extrapolation 0.069 0.957 1.23
#> 31   Marim       1 0.750 1.09         Extrapolation 0.070 0.953 1.23
#> 32   Marim       1 0.775 1.09         Extrapolation 0.070 0.950 1.23
#> 33   Marim       1 0.800 1.08         Extrapolation 0.070 0.947 1.22
#> 34   Marim       1 0.825 1.08         Extrapolation 0.069 0.945 1.21
#> 35   Marim       1 0.850 1.07         Extrapolation 0.067 0.943 1.21
#> 36   Marim       1 0.855 1.07 Observed_SC(n, gamma) 0.066 0.943 1.20
#> 37   Marim       1 0.875 1.07         Extrapolation 0.064 0.942 1.20
#> 38   Marim       1 0.876 1.07  Extrap_SC(2n, alpha) 0.064 0.942 1.20
#> 39   Marim       1 0.900 1.06         Extrapolation 0.063 0.942 1.19
#> 40   Marim       1 0.925 1.06         Extrapolation 0.061 0.942 1.18
#> 41   Marim       1 0.950 1.06         Extrapolation 0.061 0.943 1.18
#> 42   Marim       1 0.969 1.07  Extrap_SC(2n, gamma) 0.061 0.947 1.19
#> 43   Marim       1 0.975 1.07         Extrapolation 0.061 0.949 1.19
#> 44   Marim       1 1.000 1.10         Extrapolation 0.057 0.991 1.22
#> 45   Marim       2 0.500 1.10           Rarefaction 0.051 1.000 1.20
#> 46   Marim       2 0.525 1.10           Rarefaction 0.052 0.997 1.20
#> 47   Marim       2 0.550 1.10           Rarefaction 0.052 0.994 1.20
#> 48   Marim       2 0.575 1.09           Rarefaction 0.052 0.991 1.20
#> 49   Marim       2 0.600 1.09           Rarefaction 0.053 0.989 1.20
#> 50   Marim       2 0.625 1.09           Rarefaction 0.053 0.986 1.19
#> 51   Marim       2 0.650 1.09           Rarefaction 0.053 0.983 1.19
#> 52   Marim       2 0.675 1.09           Rarefaction 0.054 0.981 1.19
#> 53   Marim       2 0.696 1.08 Observed_SC(n, alpha) 0.054 0.980 1.19
#> 54   Marim       2 0.700 1.08         Extrapolation 0.054 0.979 1.19
#> 55   Marim       2 0.725 1.08         Extrapolation 0.054 0.978 1.19
#> 56   Marim       2 0.750 1.08         Extrapolation 0.054 0.979 1.19
#> 57   Marim       2 0.775 1.09         Extrapolation 0.054 0.981 1.19
#> 58   Marim       2 0.800 1.09         Extrapolation 0.053 0.984 1.19
#> 59   Marim       2 0.825 1.09         Extrapolation 0.053 0.986 1.19
#> 60   Marim       2 0.850 1.09         Extrapolation 0.053 0.987 1.19
#> 61   Marim       2 0.855 1.09 Observed_SC(n, gamma) 0.053 0.987 1.19
#> 62   Marim       2 0.875 1.09         Extrapolation 0.053 0.988 1.20
#> 63   Marim       2 0.876 1.09  Extrap_SC(2n, alpha) 0.053 0.988 1.20
#> 64   Marim       2 0.900 1.09         Extrapolation 0.054 0.987 1.20
#> 65   Marim       2 0.925 1.09         Extrapolation 0.055 0.986 1.20
#> 66   Marim       2 0.950 1.09         Extrapolation 0.056 0.984 1.20
#> 67   Marim       2 0.969 1.09  Extrap_SC(2n, gamma) 0.057 0.982 1.21
#> 68   Marim       2 0.975 1.09         Extrapolation 0.057 0.982 1.21
#> 69   Marim       2 1.000 1.09         Extrapolation 0.061 0.969 1.21

Run the following code to plot the two types of rarefaction and extrapolation (R/E) sampling curves:


## Coverage-based R/E curves for taxonomic gamma, alpha and beta diversity 
ggiNEXTbeta3D(output_TDc_abun, type = 'B')


## Coverage-based R/E curves for four taxonomic dissimilarity indices
ggiNEXTbeta3D(output_TDc_abun, type = 'D')

The following code generates size-based rarefaction and extrapolation (R/E) sampling curves for taxonomic gamma and alpha diversity:


## Size-based R/E curves with taxonomic gamma and alpha diversity
data(Brazil_rainforests)
output_TDs_abun = iNEXTbeta3D(data = Brazil_rainforests, diversity = 'TD', 
                              datatype = 'abundance', base = "size", nboot = 10)

ggiNEXTbeta3D(output_TDs_abun)

EXAMPLE 2: Incidence raw data with default sample sizes or coverage values (all assemblages as a whole)

We use the Second_growth_forests incidence raw data to compute coverage-based standardized taxonomic gamma, alpha, and beta diversity, as well as four dissimilarity indices, by setting base = "coverage". Size-based standardized taxonomic gamma and alpha diversity can also be computed by setting base = "size". Run the following code to perform coverage-based analysis of the incidence raw data. The resulting output is similar in structure to that for abundance data and is therefore omitted.


## Coverage-based R/E Analysis with taxonomic diversity for incidence raw data (all assemblages as a whole)
data(Second_growth_forests)

output_TDc_inci = iNEXTbeta3D(data = Second_growth_forests, diversity = 'TD', 
                              datatype = 'incidence_raw', base = "coverage", nboot = 10)
output_TDc_inci

The following code generates coverage-based R/E sampling curves for taxonomic gamma, alpha, and beta diversity among the three years (2005, 2011, and 2017) in the two second-growth forests (CR and JE):


## Coverage-based R/E curves with taxonomic gamma, alpha and beta diversity 

ggiNEXTbeta3D(output_TDc_inci, type = 'B')

The following code generates size-based R/E sampling curves for taxonomic gamma and alpha diversity:


## Size-based R/E curves with taxonomic gamma and alpha diversity (all assemblages as a whole)
data(Second_growth_forests)

output_TDs_inci = iNEXTbeta3D(data = Second_growth_forests, diversity = 'TD', 
                              datatype = 'incidence_raw', base = "size", nboot = 10)

ggiNEXTbeta3D(output_TDs_inci)

EXAMPLE 3: Incidence raw data with default sample sizes or coverage values for every pair of assemblages

When a dataset contains N assemblages (N > 2), as in the Second_growth_forests incidence dataset (N = 3), the default output from iNEXTbeta3D() provides gamma, alpha, and beta diversity, as well as dissimilarity, for all N assemblages considered as a whole. By setting by_pair = TRUE, users can obtain the corresponding output for every pair of assemblages. Here, we use the Second_growth_forests data to assess pairwise temporal beta diversity among three years (2005 vs. 2011, 2005 vs. 2017, and 2011 vs. 2017) within each forest.


## Coverage-based R/E Analysis for every pair of assemblages with taxonomic diversity for incidence raw data
data(Second_growth_forests)

output_TDc_inci_pair = iNEXTbeta3D(data = Second_growth_forests, diversity = 'TD',
                                   datatype = 'incidence_raw', base = "coverage",
                                   nboot = 10, by_pair = TRUE)
output_TDc_inci_pair

The output provides gamma, alpha, and beta diversity, as well as four dissimilarity indices, for each pair of assemblages. In the output, the Pair column identifies the two assemblages being compared. The corresponding coverage-based R/E sampling curves for taxonomic gamma, alpha, and beta diversity can be plotted using ggiNEXTbeta3D(). When by_pair = TRUE, the function returns a named list of plots, with one plot for each dataset. Within each plot, different pairs of assemblages are distinguished by color.

For example, the following code displays the coverage-based R/E sampling curves for taxonomic gamma, alpha, and beta diversity for every pair of years in the CR forest:


## Coverage-based R/E curves for every pair of assemblages in CR
plot_TDc_inci_pair = ggiNEXTbeta3D(output_TDc_inci_pair, type = 'B')
plot_TDc_inci_pair[["CR_2005_2011_2017"]]

Similarly, the corresponding coverage-based R/E sampling curves for every pair of years in the JE forest can be displayed using the following code:


## Coverage-based R/E curves for every pair of assemblages in JE
plot_TDc_inci_pair[["JE_2005_2011_2017"]]

The following code generates size-based R/E sampling curves for taxonomic gamma and alpha diversity for every pair of assemblages in the CR forest:


## Size-based R/E Analysis for every pair of assemblages with taxonomic diversity
data(Second_growth_forests)

output_TDs_inci_pair = iNEXTbeta3D(data = Second_growth_forests, diversity = 'TD',
                                   datatype = 'incidence_raw', base = "size",
                                   nboot = 10, by_pair = TRUE)

plot_TDs_inci_pair = ggiNEXTbeta3D(output_TDs_inci_pair)
plot_TDs_inci_pair[["CR_2005_2011_2017"]]

PHYLOGENETIC DIVERSITY (PD): RAREFACTION/EXTRAPOLATION VIA EXAMPLES

EXAMPLE 4: Abundance data with default sample sizes or coverage values

As with taxonomic diversity, iNEXTbeta3D() computes coverage-based standardized phylogenetic gamma, alpha, and beta diversity, as well as four dissimilarity indices. It also computes size-based standardized phylogenetic gamma and alpha diversity. To perform PD analysis, a phylogenetic tree spanning all species observed across the input datasets is required. The species names (or identification codes) in the phylogenetic tree must exactly match those in the corresponding abundance or incidence data. Two types of phylogenetic rarefaction and extrapolation (R/E) sampling curves, coverage-based and size-based, are also provided.

The required argument for PD analysis is PDtree. For example, the phylogenetic tree spanning all species observed in the Marim and Rebio2 fragments is available in the R dataset Brazil_tree. We therefore specify PDtree = Brazil_tree. Three optional arguments are PDtype, PDreftime, and PDcut_number. Two phylogenetic diversity measures are available through PDtype: "meanPD" (the effective number of equally divergent lineages at a specified reference time) and "AUC" (an overall measure obtained by integrating and averaging meanPD across reference times from zero to the root age of the phylogenetic tree). The default is PDtype = "AUC".

For PDtype = "meanPD", PDreftime specifies a fixed reference time for phylogenetic diversity. If PDreftime = NULL (default), the root age of the phylogenetic tree spanning the pooled observed species is used. At the specified reference time, meanPD represents the effective number of lineages, whereas PD represents the effective total branch length. The two measures are related by PD = meanPD × Reftime.

When the root age is large, meanPD may yield relatively low values because the effective number of lineages with long lineage lengths is necessarily low. The AUC approach provides an overall measure that integrates information across all reference times and is therefore less dependent on the choice of a single reference time.

For PDtype = "AUC", PDcut_number specifies the number of equally spaced integration points used to numerically compute the AUC over the interval from zero to the root age. The default is 30, and the minimum is 3. A larger number of integration points may improve the accuracy of the numerical approximation.

Run the following code to perform coverage-based PD analysis using the default AUC measure:


## Coverage-based R/E Analysis with phylogenetic diversity for abundance data
data(Brazil_rainforests)
data(Brazil_tree)

output_PDc_abun = iNEXTbeta3D(data = Brazil_rainforests, diversity = 'PD', 
                              datatype = 'abundance', base = "coverage", nboot = 10, 
                              PDtree = Brazil_tree, PDtype = 'AUC', PDcut_number = 30)
output_PDc_abun

The following code generates coverage-based R/E sampling curves for phylogenetic AUC-based gamma, alpha, and beta diversity:


## Coverage-based R/E sampling curves for phylogenetic gamma, alpha and beta diversity
ggiNEXTbeta3D(output_PDc_abun, type = 'B')

The following code generates size-based R/E sampling curves for phylogenetic AUC-based gamma and alpha diversity:


## Size-based R/E curves for phylogenetic gamma and alpha diversity
data(Brazil_rainforests)
data(Brazil_tree)

output_PDs_abun = iNEXTbeta3D(data = Brazil_rainforests, diversity = 'PD', 
                              datatype = 'abundance', base = "size", nboot = 10, 
                              PDtree = Brazil_tree, PDtype = 'AUC', PDcut_number = 30)
ggiNEXTbeta3D(output_PDs_abun)

FUNCTIONAL DIVERSITY (FD): RAREFACTION/EXTRAPOLATION VIA EXAMPLES

EXAMPLE 5: Abundance data with default sample sizes or coverage values

As with taxonomic and phylogenetic diversity, iNEXTbeta3D() computes coverage-based standardized functional gamma, alpha, and beta diversity, as well as four dissimilarity indices. It also computes size-based standardized functional gamma and alpha diversity. To perform FD analysis, a species pairwise distance matrix covering all species observed across the input datasets is required. The species names (or identification codes) in the distance matrix must exactly match those in the corresponding abundance or incidence data. Two types of functional rarefaction and extrapolation (R/E) sampling curves, coverage-based and size-based, are also provided.

The required argument for FD analysis is FDdistM. For example, the species pairwise distance matrix for all species observed in the Marim and Rebio2 fragments is available in the R dataset Brazil_distM. We therefore specify FDdistM = Brazil_distM. Three optional arguments are:

FDtype: Two functional diversity measures are available. FDtype = "tau_value" computes FD at a specified threshold of functional distinctiveness (tau), whereas FDtype = "AUC" provides an overall FD measure obtained by integrating the tau-profile over threshold values from 0 to 1. The default is FDtype = "AUC".
FDtau: A numeric value specifying the threshold of functional distinctiveness used to compute FD when FDtype = "tau_value". The tau value determines the level at which species are considered functionally distinct. If FDtau = NULL (default), the threshold is set to the mean pairwise functional distance between two individuals randomly selected from the pooled data across all datasets (i.e., Rao’s quadratic entropy).
FDcut_number: An integer specifying the number of equally spaced integration points used to numerically compute the AUC over the interval from 0 to 1. The default is 30, and the minimum is 3. A larger number of integration points may improve the accuracy of the numerical approximation.

Run the following code to perform coverage-based FD analysis using the default AUC measure. The numerical output is omitted because its structure is similar to that of the preceding TD and PD examples; the corresponding R/E sampling curves are shown below.


## Coverage-based R/E Analysis with functional diversity for abundance data - FDtype = 'AUC' (area 
## under curve) by considering all threshold values between zero and one
data(Brazil_rainforests)
data(Brazil_distM)

output_FDc_abun = iNEXTbeta3D(data = Brazil_rainforests, diversity = 'FD', 
                              datatype = 'abundance', base = "coverage", nboot = 10, 
                              FDdistM = Brazil_distM, FDtype = 'AUC', FDcut_number = 30)
output_FDc_abun

The following code generates coverage-based R/E sampling curves for functional AUC-based gamma, alpha, and beta diversity:


## Coverage-based R/E sampling curves for functional gamma, alpha and beta diversity
ggiNEXTbeta3D(output_FDc_abun, type = 'B')

The following code generates size-based R/E sampling curves for functional AUC-based gamma and alpha diversity:


## Size-based R/E curves for functional gamma and alpha diversity
data(Brazil_rainforests)
data(Brazil_distM)

output_FDs_abun = iNEXTbeta3D(data = Brazil_rainforests, diversity = 'FD', 
                              datatype = 'abundance', base = "size", nboot = 10, 
                              FDdistM = Brazil_distM, FDtype = 'AUC', FDcut_number = 30)
ggiNEXTbeta3D(output_FDs_abun)

DATA INFORMATION: FUNCTION DataInfobeta3D()

The function DataInfobeta3D() provides basic information about (1) the reference sample from each individual assemblage, (2) the gamma reference sample from the pooled assemblage, and (3) the alpha reference sample from the joint assemblage. The function and its default arguments are shown below:

DataInfobeta3D(data, diversity = "TD", datatype = "abundance",
               PDtree = NULL, PDreftime = NULL, PDtype = "meanPD",
               FDdistM = NULL, FDtype = "AUC", FDtau = NULL, by_pair = FALSE)

Arguments shared by DataInfobeta3D() and iNEXTbeta3D() have the same definitions as described above. The function returns basic information for each of the three diversity dimensions, including the reference sample size, observed species richness, estimated sample coverage at the reference and doubled sample sizes (SC(n) and SC(2n) for abundance data, or SC(T) and SC(2T) for incidence data), and other relevant statistics. We use the Brazil_rainforests dataset to demonstrate the function for each diversity dimension.

## Data information for taxonomic diversity
data(Brazil_rainforests)
DataInfobeta3D(data = Brazil_rainforests, diversity = 'TD', datatype = 'abundance')

#>   Dataset        Assemblage   n S.obs SC(n) SC(2n) f1 f2 f3 f4 f5
#> 1   Marim              Edge 158    84 0.691  0.852 49 18  8  4  1
#> 2   Marim          Interior 144    80 0.704  0.899 43 23  7  5  0
#> 3   Marim Pooled assemblage 302   119 0.855  0.969 44 34 17  9  7
#> 4   Marim  Joint assemblage 302   164 0.696  0.876 92 41 15  9  1
#> 5  Rebio2              Edge 162    70 0.754  0.895 40 17  4  2  0
#> 6  Rebio2          Interior 168    74 0.763  0.877 40 13  8  4  4
#> 7  Rebio2 Pooled assemblage 330   118 0.819  0.901 60 18 15  5  3
#> 8  Rebio2  Joint assemblage 330   144 0.758  0.886 80 30 12  6  4

Output description:

Dataset = the input datasets.
Pair = the pair of assemblages being compared; this column is included only when by_pair = TRUE.
Assemblage = individual assemblages, "Pooled assemblage" (for gamma diversity), or "Joint assemblage" (for alpha diversity).
n = number of observed individuals in the reference sample (sample size).
S.obs = number of observed species in the reference sample.
SC(n) = estimated sample coverage of the reference sample.
SC(2n) = estimated sample coverage of twice the reference sample size.
f1–f5 = the first five species abundance frequency counts in the reference sample.

Below, we use the demo dataset Second_growth_forests to illustrate the output of DataInfobeta3D() for incidence raw data:

## Data information for taxonomic diversity with incidence data
data(Second_growth_forests)

DataInfobeta3D(data = Second_growth_forests, diversity = 'TD', datatype = 'incidence_raw',
               by_pair = FALSE)

#>              Dataset        Assemblage   T    U S.obs SC(T) SC(2T)  Q1 Q2 Q3 Q4 Q5
#> 1  CR_2005_2011_2017         Year_2005 100  787   135 0.919  0.953  64 17 16  6  4
#> 2  CR_2005_2011_2017         Year_2011 100  768   135 0.916  0.952  65 18 12  7  5
#> 3  CR_2005_2011_2017         Year_2017 100  768   134 0.917  0.956  64 20 11  8  3
#> 4  CR_2005_2011_2017 Pooled assemblage 100  929   152 0.926  0.961  69 22 15  6  6
#> 5  CR_2005_2011_2017  Joint assemblage 100 2323   404 0.917  0.953 193 55 39 21 12
#> 6  JE_2005_2011_2017         Year_2005 100  503    71 0.955  0.979  23  9  8  4  0
#> 7  JE_2005_2011_2017         Year_2011 100  631    88 0.942  0.962  37  8  4  6  7
#> 8  JE_2005_2011_2017         Year_2017 100  659    91 0.953  0.979  31 12  8  3  5
#> 9  JE_2005_2011_2017 Pooled assemblage 100  908   108 0.964  0.985  33 14  9  6  8
#> 10 JE_2005_2011_2017  Joint assemblage 100 1793   250 0.950  0.973  91 29 20 13 12

Information description:

Dataset = the input datasets.
Pair = the pair of assemblages being compared; this column is included only when by_pair = TRUE.
Assemblage = individual assemblages, "Pooled assemblage" (for gamma diversity), or "Joint assemblage" (for alpha diversity).
T = number of sampling units in the reference sample (sample size for incidence data).
U = total number of incidences in the reference sample.
S.obs = number of observed species in the reference sample.
SC(T) = estimated sample coverage of the reference sample.
SC(2T) = estimated sample coverage of twice the reference sample size.
Q1–Q5 = the first five species incidence frequency counts in the reference sample.

To obtain the corresponding data information for every pair of assemblages within each dataset, set by_pair = TRUE and run the following code:

## Data information for taxonomic diversity for every pair of assemblages (incidence data)
data(Second_growth_forests)


DataInfobeta3D(data = Second_growth_forests, diversity = 'TD', datatype = 'incidence_raw', 
               by_pair = TRUE)

#>              Dataset                    Pair        Assemblage   T    U S.obs SC(T) SC(2T)  Q1 Q2 Q3 Q4
#> 1  CR_2005_2011_2017                                 Year_2005 100  787   135 0.919  0.953  64 17 16  6
#> 2  CR_2005_2011_2017                                 Year_2011 100  768   135 0.916  0.952  65 18 12  7
#> 3  CR_2005_2011_2017                                 Year_2017 100  768   134 0.917  0.956  64 20 11  8
#> 4  CR_2005_2011_2017 Year_2005 vs. Year_2011 Pooled assemblage 100  860   145 0.920  0.954  69 19 13  5
#> 5  CR_2005_2011_2017 Year_2005 vs. Year_2011  Joint assemblage 100 1555   270 0.917  0.952 129 35 28 13
#> 6  CR_2005_2011_2017 Year_2005 vs. Year_2017 Pooled assemblage 100  923   151 0.925  0.959  70 21 14  6
#> 7  CR_2005_2011_2017 Year_2005 vs. Year_2017  Joint assemblage 100 1555   269 0.918  0.954 128 37 27 14
#> 8  CR_2005_2011_2017 Year_2011 vs. Year_2017 Pooled assemblage 100  837   142 0.923  0.958  65 20 15  7
#> 9  CR_2005_2011_2017 Year_2011 vs. Year_2017  Joint assemblage 100 1536   269 0.917  0.954 129 38 23 15
#> 10 JE_2005_2011_2017                                 Year_2005 100  503    71 0.955  0.979  23  9  8  4
#> 11 JE_2005_2011_2017                                 Year_2011 100  631    88 0.942  0.962  37  8  4  6
#> 12 JE_2005_2011_2017                                 Year_2017 100  659    91 0.953  0.979  31 12  8  3
#> 13 JE_2005_2011_2017 Year_2005 vs. Year_2011 Pooled assemblage 100  757    96 0.951  0.969  37  8  6  7
#> 14 JE_2005_2011_2017 Year_2005 vs. Year_2011  Joint assemblage 100 1134   159 0.947  0.970  60 17 12 10
#> 15 JE_2005_2011_2017 Year_2005 vs. Year_2017 Pooled assemblage 100  864   107 0.963  0.987  32 17  9  4
#> 16 JE_2005_2011_2017 Year_2005 vs. Year_2017  Joint assemblage 100 1162   162 0.954  0.979  54 21 16  7
#> 17 JE_2005_2011_2017 Year_2011 vs. Year_2017 Pooled assemblage 100  788   100 0.958  0.981  33 13  8  6
#> 18 JE_2005_2011_2017 Year_2011 vs. Year_2017  Joint assemblage 100 1290   179 0.948  0.971  68 20 12  9

The output columns have the same definitions as those described above and are therefore not repeated here. When by_pair = TRUE, the Pair column identifies the two assemblages being compared.


## Data information for mean phylogenetic diversity
data(Brazil_rainforests)
data(Brazil_tree)
DataInfobeta3D(data = Brazil_rainforests, diversity = 'PD', datatype = 'abundance', 
               PDtype = 'meanPD', PDtree = Brazil_tree, PDreftime = NULL)

#>   Dataset        Assemblage   n S.obs SC(n) SC(2n) PD.obs meanPD.obs f1* f2*   g1   g2 Reftime
#> 1   Marim              Edge 158    84 0.691  0.852   8805       22.0  49  26 3278 2188     400
#> 2   Marim          Interior 144    80 0.704  0.899   8436       21.1  43  28 2974 1935     400
#> 3   Marim Pooled assemblage 302   119 0.855  0.969  11842       29.6  44  39 3172 2995     400
#> 4   Marim  Joint assemblage 302   164 0.696  0.876  17241       43.1  92  54 6252 4123     400
#> 5  Rebio2              Edge 162    70 0.754  0.895   7874       19.7  40  23 3648 1717     400
#> 6  Rebio2          Interior 168    74 0.763  0.877   8360       20.9  40  17 3365 1954     400
#> 7  Rebio2 Pooled assemblage 330   118 0.819  0.901  11979       29.9  60  23 5063 1637     400
#> 8  Rebio2  Joint assemblage 330   144 0.758  0.886  16234       40.6  80  40 7013 3671     400

For PDtype = "meanPD", the output information is described as follows:

Dataset, Pair, Assemblage, n, S.obs, SC(n), and SC(2n) = definitions are the same as those in the TD output.
PD.obs = the observed total branch length in the phylogenetic tree spanned by all observed species.
meanPD.obs = the observed mean phylogenetic diversity, calculated as PD.obs / Reftime.
f1*, f2* = the numbers of singletons and doubletons, respectively, in the node/branch abundance set.
g1, g2 = the total branch lengths of those singletons and doubletons, respectively, in the node/branch abundance set.
Reftime = the user-specified reference time used for phylogenetic diversity calculations. If PDreftime = NULL (default), the root age of the phylogenetic tree spanning the pooled observed species is used.


## Data information for phylogenetic diversity when all time reference points are considered
data(Brazil_rainforests)
data(Brazil_tree)

DataInfobeta3D(data = Brazil_rainforests, diversity = 'PD',
               datatype = 'abundance', PDtype = 'AUC',
               PDtree = Brazil_tree)

#>   Dataset        Assemblage   n S.obs SC(n) SC(2n) Tmax
#> 1   Marim              Edge 158    84 0.691  0.852  400
#> 2   Marim          Interior 144    80 0.704  0.899  400
#> 3   Marim Pooled assemblage 302   119 0.855  0.969  400
#> 4   Marim  Joint assemblage 302   164 0.696  0.876  400
#> 5  Rebio2              Edge 162    70 0.754  0.895  400
#> 6  Rebio2          Interior 168    74 0.763  0.877  400
#> 7  Rebio2 Pooled assemblage 330   118 0.819  0.901  400
#> 8  Rebio2  Joint assemblage 330   144 0.758  0.886  400

For PDtype = "AUC", the output information is described as follows:

Dataset, Pair, Assemblage, n, S.obs, SC(n), and SC(2n) = definitions are the same as those in the TD output.
Tmax = the integration end time for AUC-based phylogenetic diversity, corresponding to the root age of the phylogenetic tree spanning the pooled observed species. The integration starts at zero.


## Data information for functional diversity (FDtype = 'tau_value', FDtau = NULL)
data(Brazil_rainforests)
data(Brazil_distM)
DataInfobeta3D(data = Brazil_rainforests, diversity = 'FD', datatype = 'abundance', 
               FDdistM = Brazil_distM, FDtype = 'tau_value', FDtau = NULL)

#>   Dataset        Assemblage   n S.obs SC(n) SC(2n) a1* a2* h1 h2   Tau
#> 1   Marim              Edge 158    84 0.691  0.852   0   0  0  0 0.343
#> 2   Marim          Interior 144    80 0.704  0.899   0   0  0  0 0.343
#> 3   Marim Pooled assemblage 302   119 0.855  0.969   0   0  0  0 0.343
#> 4   Marim  Joint assemblage 302   164 0.696  0.876   0   0  0  0 0.343
#> 5  Rebio2              Edge 162    70 0.754  0.895   0   0  0  0 0.343
#> 6  Rebio2          Interior 168    74 0.763  0.877   0   0  0  0 0.343
#> 7  Rebio2 Pooled assemblage 330   118 0.819  0.901   0   0  0  0 0.343
#> 8  Rebio2  Joint assemblage 330   144 0.758  0.886   0   0  0  0 0.343

For FDtype = "tau_value", the output information is described as follows:

Dataset, Pair, Assemblage, n, S.obs, SC(n), and SC(2n) = definitions are the same as those in the TD output.
a1*, a2* = the numbers of singletons and doubletons, respectively, in the functional abundances at the specified threshold of functional distinctiveness (Tau).
h1, h2 = the total contributions of functional singletons and doubletons, respectively, at the specified threshold.
Tau = the user-specified threshold level for functional distinctiveness. If FDtau = NULL (default), the threshold is set to dmean, the mean distance between two individuals randomly selected from the pooled data across all datasets.

For comparison, we also consider a user-specified threshold value of FDtau = 0.2 to illustrate how the data information changes with the selected threshold.

## Data information for functional diversity (tau = 0.2)
data(Brazil_rainforests)
data(Brazil_distM)
DataInfobeta3D(data = Brazil_rainforests, diversity = 'FD',
               datatype = 'abundance', FDdistM = Brazil_distM,
               FDtype = 'tau_value', FDtau = 0.2)

#>   Dataset        Assemblage   n S.obs SC(n) SC(2n) a1* a2*   h1   h2 Tau
#> 1   Marim              Edge 158    84 0.691  0.852   1   2 1.00 1.50 0.2
#> 2   Marim          Interior 144    80 0.704  0.899   0   1 0.00 0.50 0.2
#> 3   Marim Pooled assemblage 302   119 0.855  0.969   0   1 0.00 0.50 0.2
#> 4   Marim  Joint assemblage 302   164 0.696  0.876   5   9 2.12 2.31 0.2
#> 5  Rebio2              Edge 162    70 0.754  0.895   0   1 0.00 0.50 0.2
#> 6  Rebio2          Interior 168    74 0.763  0.877   0   1 0.00 0.50 0.2
#> 7  Rebio2 Pooled assemblage 330   118 0.819  0.901   0   1 0.00 0.50 0.2
#> 8  Rebio2  Joint assemblage 330   144 0.758  0.886   6   7 2.29 1.95 0.2

Compared with the default threshold (Tau = 0.343), at which a1*, a2*, h1, and h2 are all zero, the user-specified threshold (Tau = 0.2) yields nonzero values for some of these quantities. This example illustrates that singleton- and doubleton-related statistics depend on the selected threshold of functional distinctiveness.

## Data information for functional diversity (FDtype = 'AUC')
data(Brazil_rainforests)
data(Brazil_distM)
DataInfobeta3D(data = Brazil_rainforests, diversity = 'FD', datatype = 'abundance', 
               FDdistM = Brazil_distM, FDtype = 'AUC')

#>   Dataset        Assemblage   n S.obs SC(n) SC(2n) dmin dmean  dmax
#> 1   Marim              Edge 158    84 0.691  0.852    0 0.329 0.755
#> 2   Marim          Interior 144    80 0.704  0.899    0 0.313 0.663
#> 3   Marim Pooled assemblage 302   119 0.855  0.969    0 0.323 0.755
#> 4   Marim  Joint assemblage 302   164 0.696  0.876    0 0.323 0.755
#> 5  Rebio2              Edge 162    70 0.754  0.895    0 0.376 0.659
#> 6  Rebio2          Interior 168    74 0.763  0.877    0 0.310 0.660
#> 7  Rebio2 Pooled assemblage 330   118 0.819  0.901    0 0.355 0.770
#> 8  Rebio2  Joint assemblage 330   144 0.758  0.886    0 0.355 0.770

For FDtype = "AUC", the output information is described as follows:

Dataset, Pair, Assemblage, n, S.obs, SC(n), and SC(2n) = definitions are the same as those in the TD output.
dmin = the minimum distance among all non-diagonal elements in the species pairwise distance matrix.
dmean = the mean functional distance between two individuals randomly selected from the assemblage under consideration.
dmax = the maximum pairwise distance among species in the distance matrix.

License and feedback

The iNEXT.beta3D package is licensed under the GPLv3. To help refine iNEXT.beta3D, users’ comments or feedback would be welcome (please send them to Anne Chao or report an issue on the iNEXT.beta3D GitHub iNEXT.beta3D_github).

References

Chao, A., Chiu, C.-H., Hu, K.-H., and Zelený, D. (2023a). Revisiting Alwyn H. Gentry’s forest transect data: a statistical sampling-model-based approach. Japanese Journal of Statistics and Data Science, 6, 861-884. (https://doi.org/10.1007/s42081-023-00214-1)
Chao, A., Chiu, C.-H., Villéger, S., Sun, I.-F., Thorn, S., Lin, Y.-C., Chiang, J. M., and Sherwin, W. B. (2019). An attribute-diversity approach to functional diversity, functional beta diversity, and related (dis)similarity measures. Ecological Monographs, 89, e01343.
Chao, A., Henderson, P. A., Chiu, C.-H., Moyes, F., Hu, K.-H., Dornelas, M. and Magurran, A. E. (2021). Measuring temporal change in alpha diversity: a framework integrating taxonomic, phylogenetic and functional diversity and the iNEXT.3D standardization. Methods in Ecology and Evolution, 12, 1926-1940.
Chao, A., Thorn, S., Chiu, C.-H., Moyes, F., Hu, K.-H., Chazdon, R. L., Wu, J., Magnago, L. F. S., Dornelas, M., Zelený, D., Colwell, R. K., and Magurran, A. E. (2023b). Rarefaction and extrapolation with beta diversity under a framework of Hill numbers: the iNEXT.beta3D standardization. Ecological Monographs e1588. (https://doi.org/10.1002/ecm.1588)