| Title: | Fast Spatial Species Accumulation Curves |
|---|---|
| Description: | High-performance spatial species accumulation curves using nearest-neighbor algorithms. Implements 'kNN' and 'kNCN' sampling methods with a 'C++' backend for speed. Supports Hill numbers (q=0,1,2), beta diversity partitioning (turnover/nestedness), coverage-based rarefaction and extrapolation, phylogenetic diversity (Faith's PD, mean pairwise distance, mean nearest taxon distance), functional diversity accumulation, diversity-area relationships (DAR), endemism-area curves, sampling-effort correction and fragmentation analysis, and species-area relationship (SAR) models based on extreme value theory (EVT). Multiple starting points (seeds) provide uncertainty quantification. Methods are described in 'Chao' et al. (2014) <doi:10.1890/13-0133.1>, 'Baselga' (2010) <doi:10.1111/j.1466-8238.2009.00490.x>, 'Chao' and 'Jost' (2012) <doi:10.1890/11-1952.1>, 'Faith' (1992) <doi:10.1016/0006-3207(92)91201-3>, 'Ma' (2018) <doi:10.1002/ece3.4526>, 'Borda-de-Agua' et al. (2025) <doi:10.1038/s41467-025-59239-7>, 'Hanski' et al. (2013) <doi:10.1073/pnas.1311190110>, and 'Jost' (2007) <doi:10.1890/06-1736.1>. |
| Authors: | Gilles Colling [aut, cre, cph] (ORCID: <https://orcid.org/0000-0003-3070-6066>) |
| Maintainer: | Gilles Colling <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.8.3 |
| Built: | 2026-06-20 17:32:49 UTC |
| Source: | https://github.com/cran/spacc |
Abundance-based Coverage Estimator (ACE) of total species richness (Chao & Lee 1992). Separates species into rare and abundant groups based on an abundance threshold.
ace(x, threshold = 10L)ace(x, threshold = 10L)
x |
A site-by-species matrix (abundance data). |
threshold |
Integer. Abundance threshold separating rare from abundant species. Default 10. |
ACE partitions species into rare (abundance threshold) and
abundant (abundance > threshold) groups. The estimate is:
where is the sample coverage of rare species and
is the estimated coefficient of variation.
An object of class spacc_estimate.
Chao, A. & Lee, S.M. (1992). Estimating the number of classes via sample coverage. Journal of the American Statistical Association, 87, 210-217.
species <- matrix(rpois(50 * 30, 2), nrow = 50) ace(species)species <- matrix(rpois(50 * 30, 2), nrow = 50) ace(species)
Compute Hill numbers for each site individually.
alphaDiversity(x, q = c(0, 1, 2), coords = NULL)alphaDiversity(x, q = c(0, 1, 2), coords = NULL)
x |
A site-by-species matrix (abundance data). |
q |
Numeric vector. Orders of diversity. Default |
coords |
Optional data.frame with columns |
Alpha diversity represents local (within-site) diversity. For Hill numbers:
q=0: Species richness
q=1: Exponential of Shannon entropy
q=2: Inverse Simpson concentration
If coords is NULL, a data.frame with columns for each q value.
If coords is provided, a spacc_alpha object.
Jost, L. (2007). Partitioning diversity into independent alpha and beta components. Ecology, 88, 2427-2439.
gammaDiversity() for regional diversity, diversityPartition()
for full alpha-beta-gamma decomposition
species <- matrix(rpois(50 * 30, 2), nrow = 50) alpha <- alphaDiversity(species, q = c(0, 1, 2)) head(alpha) # Mean alpha diversity colMeans(alpha)species <- matrix(rpois(50 * 30, 2), nrow = 50) alpha <- alphaDiversity(species, q = c(0, 1, 2)) head(alpha) # Mean alpha diversity colMeans(alpha)
These methods compute expected species accumulation without simulation. They are faster but don't provide spatial information.
Convert metrics to an sf object for spatial analysis and integration with the areaOfEffect package.
as_sf(x, crs = NULL)as_sf(x, crs = NULL)
x |
A |
crs |
Coordinate reference system. Default |
An sf object with POINT geometries and metric columns.
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) metrics <- spaccMetrics(species, coords) if (requireNamespace("sf", quietly = TRUE)) { metrics_sf <- as_sf(metrics) }coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) metrics <- spaccMetrics(species, coords) if (requireNamespace("sf", quietly = TRUE)) { metrics_sf <- as_sf(metrics) }
ggplot2::autoplot() methods for spacc objects. These are thin wrappers
around the corresponding plot() methods, provided for ggplot2 integration.
object |
A spacc object. |
... |
Additional arguments passed to the corresponding |
A ggplot2 object.
if (requireNamespace("ggplot2", quietly = TRUE)) { coords <- data.frame(x = runif(30), y = runif(30)) species <- matrix(rbinom(30 * 15, 1, 0.3), nrow = 30) sac <- spacc(species, coords, n_seeds = 5, progress = FALSE) ggplot2::autoplot(sac) }if (requireNamespace("ggplot2", quietly = TRUE)) { coords <- data.frame(x = runif(30), y = runif(30)) species <- matrix(rbinom(30 * 15, 1, 0.3), nrow = 30) sac <- spacc(species, coords, n_seeds = 5, progress = FALSE) ggplot2::autoplot(sac) }
Compute the distance-decay of community similarity by fitting models to pairwise beta dissimilarity as a function of geographic distance.
betaDecay( x, coords, index = c("sorensen", "jaccard"), model = c("all", "exponential", "power", "linear"), distance = c("euclidean", "haversine"), n_bins = 50L, progress = TRUE )betaDecay( x, coords, index = c("sorensen", "jaccard"), model = c("all", "exponential", "power", "linear"), distance = c("euclidean", "haversine"), n_bins = 50L, progress = TRUE )
x |
A site-by-species matrix (presence/absence or abundance), or a
|
coords |
A data.frame with columns |
index |
Character. Dissimilarity index: |
model |
Character. Decay model to fit: |
distance |
Character. Distance method: |
n_bins |
Integer. Number of distance bins for binned means in plots. Default 50. |
progress |
Logical. Show progress messages? Default |
Three decay models are available:
Exponential:
Power:
Linear:
The half-life (exponential model) is the distance at which similarity
decays to half its initial value: .
An object of class spacc_beta_decay containing:
pairs |
Data.frame with columns |
fits |
Named list of fitted model objects |
best_model |
Name of best model by AIC |
half_life |
Distance at which similarity halves (exponential model only) |
coefficients |
Data.frame of model coefficients and AIC |
index |
Dissimilarity index used |
n_sites |
Number of sites |
n_pairs |
Number of pairwise comparisons |
Nekola, J.C. & White, P.S. (1999). The distance decay of similarity in biogeography and ecology. Journal of Biogeography, 26, 867-878.
Morlon, H., Chuyong, G., Condit, R., et al. (2008). A general framework for the distance-decay of similarity in ecological communities. Ecology Letters, 11, 904-917.
spaccBeta() for spatial beta diversity accumulation,
distanceDecay() for species similarity decay
coords <- data.frame(x = runif(30), y = runif(30)) species <- matrix(rbinom(30 * 20, 1, 0.3), nrow = 30) bd <- betaDecay(species, coords) print(bd) plot(bd)coords <- data.frame(x = runif(30), y = runif(30)) species <- matrix(rbinom(30 * 20, 1, 0.3), nrow = 30) bd <- betaDecay(species, coords) print(bd) plot(bd)
Bootstrap estimator of total species richness (Smith & van Belle 1984). Uses species detection probabilities to estimate undetected species.
bootstrap_richness(x, n_boot = 200L)bootstrap_richness(x, n_boot = 200L)
x |
A site-by-species matrix (presence/absence or abundance). Automatically binarized. |
n_boot |
Integer. Number of bootstrap replicates for SE. Default 200. |
The bootstrap estimator is:
where is the proportion of sites where species occurs.
An object of class spacc_estimate.
Smith, E.P. & van Belle, G. (1984). Nonparametric estimation of species richness. Biometrics, 40, 119-129.
species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) bootstrap_richness(species, n_boot = 100)species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) bootstrap_richness(species, n_boot = 100)
Combine spacc Objects
## S3 method for class 'spacc' c(...)## S3 method for class 'spacc' c(...)
... |
Named |
A grouped spacc object with per-group curves.
Combine multiple spacc_beta objects by stacking their curve matrices.
All objects must have the same number of sites and index.
## S3 method for class 'spacc_beta' c(...)## S3 method for class 'spacc_beta' c(...)
... |
|
A combined spacc_beta object with more seeds.
Combine multiple spacc_coverage objects by stacking their curve matrices.
All objects must have the same number of sites.
## S3 method for class 'spacc_coverage' c(...)## S3 method for class 'spacc_coverage' c(...)
... |
|
A combined spacc_coverage object with more seeds.
Combine multiple spacc_hill objects by stacking their curve matrices.
All objects must have the same number of sites.
## S3 method for class 'spacc_hill' c(...)## S3 method for class 'spacc_hill' c(...)
... |
|
A combined spacc_hill object with more seeds.
Estimate total species richness from abundance data using the Chao1 estimator (Chao 1984). Uses the number of singletons and doubletons to estimate undetected species.
chao1(x)chao1(x)
x |
A site-by-species matrix (abundance data). Columns are pooled across sites. |
The Chao1 estimator is:
where is the number of singletons (species with total abundance 1)
and is the number of doubletons (abundance 2).
When , the bias-corrected form is used:
An object of class spacc_estimate with components:
Name of the estimator ("chao1")
Estimated total richness
Standard error of the estimate
Lower 95 percent confidence bound
Upper 95 percent confidence bound
Observed species richness
List with f1 (singletons) and f2 (doubletons)
Chao, A. (1984). Nonparametric estimation of the number of classes in a population. Scandinavian Journal of Statistics, 11, 265-270.
Chao, A. (1987). Estimating the population size for capture-recapture data with unequal catchability. Biometrics, 43, 783-791.
chao2() for incidence-based estimation, ace() for
abundance-based coverage estimator
species <- matrix(rpois(50 * 30, 2), nrow = 50) chao1(species)species <- matrix(rpois(50 * 30, 2), nrow = 50) chao1(species)
Estimate total species richness from incidence (presence/absence) data using the Chao2 estimator (Chao 1987).
chao2(x)chao2(x)
x |
A site-by-species matrix (presence/absence or abundance). Automatically binarized. |
The Chao2 estimator is the incidence-based analogue of Chao1:
where is the number of uniques (species found at exactly 1 site)
and is the number of duplicates (species found at exactly 2 sites).
An object of class spacc_estimate.
Chao, A. (1987). Estimating the population size for capture-recapture data with unequal catchability. Biometrics, 43, 783-791.
chao1() for abundance-based estimation
species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) chao2(species)species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) chao2(species)
Compute the expected species accumulation curve using the Coleman method (Coleman et al. 1982). This is an analytical formula, no simulation needed.
coleman(x)coleman(x)
x |
A site-by-species matrix (presence/absence or abundance). |
A data.frame with columns: sites, expected, sd
Coleman, B.D., Mares, M.A., Willig, M.R. & Hsieh, Y.H. (1982). Randomness, area, and species richness. Ecology, 63, 1121-1133.
Compute the species accumulation curve in the order sites appear in the data (no randomization). Useful for understanding how data was collected.
collector(x)collector(x)
x |
A site-by-species matrix. |
A data.frame with columns: sites, species
Test whether two species accumulation curves differ significantly.
compare( x, y, method = c("permutation", "bootstrap", "auc"), normalize = FALSE, n_perm = 999L, ... )compare( x, y, method = c("permutation", "bootstrap", "auc"), normalize = FALSE, n_perm = 999L, ... )
x |
A |
y |
A |
method |
Character. Comparison method: |
normalize |
Logical. If |
n_perm |
Integer. Number of permutations/bootstrap replicates. Default 999. |
... |
Additional arguments passed to comparison methods. |
An object of class spacc_comp containing:
x_name, y_name
|
Names of compared objects |
auc_diff |
Difference in area under curve |
p_value |
P-value from permutation test |
saturation_diff |
Difference in saturation points |
method |
Comparison method used |
normalized |
Whether curves were normalized before comparison |
Colwell, R.K., Mao, C.X. & Chang, J. (2004). Interpolating, extrapolating, and comparing incidence-based species accumulation curves. Ecology, 85, 2717-2727.
Gotelli, N.J. & Colwell, R.K. (2011). Estimating species richness. In: Biological Diversity: Frontiers in Measurement and Assessment, pp. 39-54. Oxford University Press.
coords <- data.frame(x = runif(50), y = runif(50)) sp_a <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sp_b <- matrix(rbinom(50 * 30, 1, 0.5), nrow = 50) sac_a <- spacc(sp_a, coords, n_seeds = 10) sac_b <- spacc(sp_b, coords, n_seeds = 10) comp <- compare(sac_a, sac_b) print(comp)coords <- data.frame(x = runif(50), y = runif(50)) sp_a <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sp_b <- matrix(rbinom(50 * 30, 1, 0.5), nrow = 50) sac_a <- spacc(sp_a, coords, n_seeds = 10) sac_b <- spacc(sp_b, coords, n_seeds = 10) comp <- compare(sac_a, sac_b) print(comp)
Fit all (or a subset of) asymptotic species-area models and compare them using AIC, BIC, delta-AIC, and Akaike weights.
compareModels( object, models = c("michaelis-menten", "lomolino", "asymptotic", "weibull", "logistic", "evt"), ... )compareModels( object, models = c("michaelis-menten", "lomolino", "asymptotic", "weibull", "logistic", "evt"), ... )
object |
A |
models |
Character vector of models to fit. Defaults to all six:
|
... |
Additional arguments passed to |
An object of class spacc_model_compare containing:
table |
Data frame with model comparison statistics |
fits |
Named list of |
best_model |
Name of the best model by AIC |
data |
Mean-curve data frame used for fitting |
spacc |
Original spacc object |
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, coords) cm <- compareModels(sac) print(cm)coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, coords) cm <- compareModels(sac) print(cm)
Compute the ratio of observed to estimated diversity across diversity orders, measuring how complete a sample is at each level of the Hill number spectrum.
completenessProfile(x, q = seq(0, 2, by = 0.2), coords = NULL)completenessProfile(x, q = seq(0, 2, by = 0.2), coords = NULL)
x |
A site-by-species matrix (abundance data). |
q |
Numeric vector. Orders of diversity to evaluate. Default
|
coords |
Optional data.frame with columns |
Sample completeness at order q is:
where is the observed Hill number and is
the estimated asymptotic Hill number.
Completeness near 1 means the sample captures most of the true diversity at that order. Completeness typically increases with q because dominant species are detected early.
Asymptotic estimators used:
q = 0: Chao1 estimator
q = 1: Chao & Jost (2015) entropy estimator
q = 2: Inverse Simpson estimator with bias correction
Other q: Interpolated between adjacent integer estimates
When coords is provided, per-site completeness is computed by treating
each site's abundance vector as an independent sample.
An object of class spacc_completeness containing:
completeness |
Named numeric vector of completeness ratios per q |
observed |
Named numeric vector of observed Hill numbers per q |
estimated |
Named numeric vector of estimated asymptotic Hill numbers per q |
per_site |
Matrix of per-site completeness (sites x q values), or |
q |
Vector of diversity orders |
coords |
Coordinates if provided |
n_sites |
Number of sites |
n_species |
Number of species |
Chao, A. & Jost, L. (2012). Coverage-based rarefaction and extrapolation: standardizing samples by completeness rather than size. Ecology, 93, 2533-2547.
Chao, A. & Jost, L. (2015). Estimating diversity and entropy profiles via discovery rates of new species. Methods in Ecology and Evolution, 6, 873-882.
diversityProfile() for observed diversity profiles,
chao1() for richness estimation
species <- matrix(rpois(50 * 30, 2), nrow = 50) comp <- completenessProfile(species) print(comp) plot(comp)species <- matrix(rpois(50 * 30, 2), nrow = 50) comp <- completenessProfile(species) print(comp) plot(comp)
Extend the classic species-area relationship (SAR) to a diversity-area relationship using Hill numbers of any order q. Instead of plotting species richness vs. sites, this plots effective diversity vs. cumulative area.
dar( x, coords, q = c(0, 1, 2), n_seeds = 50L, method = "knn", area_method = c("convex_hull", "voronoi", "count"), distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )dar( x, coords, q = c(0, 1, 2), n_seeds = 50L, method = "knn", area_method = c("convex_hull", "voronoi", "count"), distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )
x |
A site-by-species matrix (abundance data recommended). |
coords |
A data.frame with columns |
q |
Numeric vector. Diversity orders. Default |
n_seeds |
Integer. Number of random starting points. Default 50. |
method |
Character. Accumulation method. Default |
area_method |
Character. How to estimate cumulative area:
|
distance |
Character. Distance method. Default |
parallel |
Logical. Use parallel processing? Default |
n_cores |
Integer. Number of cores. |
progress |
Logical. Show progress? Default |
seed |
Integer. Random seed. |
The DAR (Ma, 2018) generalizes the SAR by replacing species richness (q=0) with Hill numbers of any order. This reveals how different facets of diversity scale with area:
q=0 (richness) recovers the classic SAR
q=1 (Shannon) shows how common species diversity scales
q=2 (Simpson) shows how dominant species diversity scales
An object of class spacc_dar containing:
hill |
A |
area |
Matrix of cumulative areas (n_seeds x n_sites) |
q |
Diversity orders used |
area_method |
Method used for area estimation |
Ma, Z.S. (2018). DAR (diversity-area relationship): extending classic SAR for biodiversity and biogeography analyses. Ecology and Evolution, 8, 10023-10038.
Arrhenius, O. (1921). Species and area. Journal of Ecology, 9, 95-99.
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rpois(50 * 30, 2), nrow = 50) dar <- dar(species, coords, q = c(0, 1, 2)) plot(dar)coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rpois(50 * 30, 2), nrow = 50) dar <- dar(species, coords, q = c(0, 1, 2)) plot(dar)
Analyze how species richness changes with distance from focal points.
distanceDecay( x, coords, n_seeds = 50L, breaks = NULL, distance = c("euclidean", "haversine"), progress = TRUE, seed = NULL )distanceDecay( x, coords, n_seeds = 50L, breaks = NULL, distance = c("euclidean", "haversine"), progress = TRUE, seed = NULL )
x |
A site-by-species matrix. |
coords |
A data.frame with x and y columns. |
n_seeds |
Integer. Number of focal points. |
breaks |
Numeric vector. Distance thresholds. Default auto-calculated. |
distance |
Character. Distance method. |
progress |
Logical. Show progress? |
seed |
Integer. Random seed. |
An object of class spacc_decay with distance-species relationships.
Nekola, J.C. & White, P.S. (1999). The distance decay of similarity in biogeography and ecology. Journal of Biogeography, 26, 867-878.
Soininen, J., McDonald, R. & Hillebrand, H. (2007). The distance decay of similarity in ecological communities. Ecography, 30, 3-12.
Pre-compute pairwise distances between sites for reuse across multiple
spacc() calls. Supports sf objects with accurate geodesic distances
for global-scale studies.
distances(x, method = NULL, fun = NULL, which = NULL)distances(x, method = NULL, fun = NULL, which = NULL)
x |
Site locations. Can be:
|
method |
Character. Distance method:
|
fun |
Optional custom distance function. Must take two coordinate
vectors (x, y) and return a distance matrix. Overrides |
which |
For sf objects, column name containing the geometry. Default uses active geometry. |
For continental and global-scale studies, use sf objects with geographic CRS (e.g., EPSG:4326). The function will automatically use accurate geodesic distances via the S2 spherical geometry library.
For smaller study areas with projected coordinates (UTM, etc.), Euclidean distance is appropriate and faster.
An object of class spacc_dist containing the distance matrix
with coordinates stored as an attribute.
coords <- data.frame(x = runif(50), y = runif(50)) d <- distances(coords) # Reuse for multiple analyses species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, d)coords <- data.frame(x = runif(50), y = runif(50)) d <- distances(coords) # Reuse for multiple analyses species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, d)
Decompose regional (gamma) diversity into local (alpha) and turnover (beta) components using multiplicative partitioning of Hill numbers.
diversityPartition(x, q = c(0, 1, 2), weights = "equal", coords = NULL)diversityPartition(x, q = c(0, 1, 2), weights = "equal", coords = NULL)
x |
A site-by-species matrix (abundance data). |
q |
Numeric vector. Orders of diversity. Default |
weights |
Character or numeric. Site weights for alpha calculation:
|
coords |
Optional data.frame with columns |
This function implements Jost (2007) multiplicative partitioning:
Where:
Alpha: Mean effective number of species per site
Beta: Effective number of distinct communities (1 = all identical, n_sites = all completely different)
Gamma: Total effective number of species in the region
Beta diversity is interpreted as the effective number of communities:
Beta = 1: All sites have identical composition
Beta = n_sites: Sites share no species
An object of class spacc_partition containing:
alpha |
Mean alpha diversity (effective number of species per site) |
beta |
Beta diversity (effective number of communities) |
gamma |
Gamma diversity (regional species pool) |
q |
Orders of diversity |
n_sites |
Number of sites |
n_species |
Total species count |
Jost, L. (2007). Partitioning diversity into independent alpha and beta components. Ecology, 88, 2427-2439.
Chao, A., Chiu, C.H. & Jost, L. (2014). Unifying species diversity, phylogenetic diversity, functional diversity, and related similarity and differentiation measures through Hill numbers. Annual Review of Ecology, Evolution, and Systematics, 45, 297-324.
alphaDiversity(), gammaDiversity(), spaccBeta() for
spatial beta diversity accumulation
# Simulated community data species <- matrix(rpois(50 * 30, 2), nrow = 50) # Partition diversity part <- diversityPartition(species, q = c(0, 1, 2)) print(part) # Beta near 1 = homogeneous region # Beta near n_sites = heterogeneous region# Simulated community data species <- matrix(rpois(50 * 30, 2), nrow = 50) # Partition diversity part <- diversityPartition(species, q = c(0, 1, 2)) print(part) # Beta near 1 = homogeneous region # Beta near n_sites = heterogeneous region
Compute Hill numbers across a continuous range of diversity orders (q), producing a diversity profile for each site and the regional pool.
diversityProfile( x, q = seq(0, 3, by = 0.1), type = c("both", "per_site", "regional"), coords = NULL )diversityProfile( x, q = seq(0, 3, by = 0.1), type = c("both", "per_site", "regional"), coords = NULL )
x |
A site-by-species matrix (abundance data). |
q |
Numeric vector. Orders of diversity to evaluate. Default
|
type |
Character. What to compute: |
coords |
Optional data.frame with columns |
A diversity profile plots effective number of species as a function of the
order q. The key property is that Hill numbers are non-increasing in q:
for .
q = 0: Species richness (insensitive to abundance)
q = 1: Exponential of Shannon entropy (all species weighted equally by their proportional abundance)
q = 2: Inverse Simpson concentration (emphasizes dominant species)
q > 2: Increasingly dominated by common species
An object of class spacc_profile containing:
per_site |
Matrix of per-site diversity (sites x q values), or |
regional |
Named numeric vector of gamma diversity per q, or |
q |
Vector of diversity orders used |
coords |
Coordinates if provided |
n_sites |
Number of sites |
n_species |
Number of species |
Leinster, T. & Cobbold, C.A. (2012). Measuring diversity: the importance of species similarity. Ecology, 93, 477-489.
Chao, A., Chiu, C.H. & Jost, L. (2014). Unifying species diversity, phylogenetic diversity, functional diversity, and related similarity and differentiation measures through Hill numbers. Annual Review of Ecology, Evolution, and Systematics, 45, 297-324.
alphaDiversity() for per-site values at specific q,
gammaDiversity() for regional diversity, evenness() for evenness
profiles
species <- matrix(rpois(50 * 30, 2), nrow = 50) prof <- diversityProfile(species) print(prof) plot(prof)species <- matrix(rpois(50 * 30, 2), nrow = 50) prof <- diversityProfile(species) print(prof) plot(prof)
Compute functional Hill numbers (Leinster & Cobbold 2012) across a continuous range of diversity orders (q), producing a functional diversity profile based on trait similarity.
diversityProfileFunc( x, traits, q = seq(0, 3, by = 0.1), type = c("both", "per_site", "regional"), dist_method = c("euclidean", "gower"), normalize = TRUE, coords = NULL )diversityProfileFunc( x, traits, q = seq(0, 3, by = 0.1), type = c("both", "per_site", "regional"), dist_method = c("euclidean", "gower"), normalize = TRUE, coords = NULL )
x |
A site-by-species matrix (abundance data). Column names must match
row names in |
traits |
A data.frame of species traits. Row names must match column
names in |
q |
Numeric vector. Orders of diversity. Default |
type |
Character. What to compute: |
dist_method |
Character. Distance method for trait matrix:
|
normalize |
Logical. Normalize distances to [0, 1]? Default |
coords |
Optional data.frame with |
Functional Hill numbers (Leinster & Cobbold 2012) incorporate trait similarity via a similarity matrix Z = 1 - D. When all species are maximally dissimilar (Z = identity), this reduces to standard Hill numbers.
An object of class spacc_profile with $profile_type = "functional".
Leinster, T. & Cobbold, C.A. (2012). Measuring diversity: the importance of species similarity. Ecology, 93, 477-489.
diversityProfile() for taxonomic profiles,
diversityProfilePhylo() for phylogenetic profiles
species <- matrix(rpois(20 * 10, 2), nrow = 20, dimnames = list(NULL, paste0("sp", 1:10))) traits <- data.frame( body_size = rnorm(10), diet = rnorm(10), row.names = paste0("sp", 1:10) ) prof <- diversityProfileFunc(species, traits) print(prof)species <- matrix(rpois(20 * 10, 2), nrow = 20, dimnames = list(NULL, paste0("sp", 1:10))) traits <- data.frame( body_size = rnorm(10), diet = rnorm(10), row.names = paste0("sp", 1:10) ) prof <- diversityProfileFunc(species, traits) print(prof)
Compute phylogenetic Hill numbers (Chao et al. 2010) across a continuous range of diversity orders (q), producing a phylogenetic diversity profile.
diversityProfilePhylo( x, tree, q = seq(0, 3, by = 0.1), type = c("both", "per_site", "regional"), coords = NULL )diversityProfilePhylo( x, tree, q = seq(0, 3, by = 0.1), type = c("both", "per_site", "regional"), coords = NULL )
x |
A site-by-species matrix (abundance data). Column names must match tip labels in the phylogeny. |
tree |
An |
q |
Numeric vector. Orders of diversity. Default |
type |
Character. What to compute: |
coords |
Optional data.frame with |
Phylogenetic Hill numbers (Chao et al. 2010) weight branches by their evolutionary distance. At q=0 this approximates normalized Faith's PD. Higher q values increasingly emphasize common lineages.
An object of class spacc_profile with $profile_type = "phylogenetic".
Chao, A., Chiu, C.H. & Jost, L. (2010). Phylogenetic diversity measures based on Hill numbers. Philosophical Transactions of the Royal Society B, 365, 3599-3609.
diversityProfile() for taxonomic profiles,
diversityProfileFunc() for functional profiles
if (requireNamespace("ape", quietly = TRUE)) { species <- matrix(rpois(20 * 10, 2), nrow = 20, dimnames = list(NULL, paste0("sp", 1:10))) tree <- ape::rcoal(10, tip.label = paste0("sp", 1:10)) prof <- diversityProfilePhylo(species, tree) print(prof) }if (requireNamespace("ape", quietly = TRUE)) { species <- matrix(rpois(20 * 10, 2), nrow = 20, dimnames = list(NULL, paste0("sp", 1:10))) tree <- ape::rcoal(10, tip.label = paste0("sp", 1:10)) prof <- diversityProfilePhylo(species, tree) print(prof) }
Compute species evenness across sites using Hill-based, Pielou, or Simpson evenness measures.
evenness( x, q = seq(0.1, 3, by = 0.1), type = c("hill", "pielou", "simpson"), coords = NULL )evenness( x, q = seq(0.1, 3, by = 0.1), type = c("hill", "pielou", "simpson"), coords = NULL )
x |
A site-by-species matrix (abundance data). |
q |
Numeric vector. Orders of diversity for Hill evenness.
Default |
type |
Character. Evenness type: |
coords |
Optional data.frame with columns |
All evenness measures are bounded in :
0 = maximally uneven (one dominant species)
1 = perfectly even (all species equally abundant)
Hill evenness (Jost 2010):
This divides the effective number of species at order q by species richness.
Pielou's J (Pielou 1966):
Simpson evenness:
An object of class spacc_evenness containing:
per_site |
Matrix or vector of per-site evenness values |
regional |
Regional (pooled) evenness |
q |
Orders used (for Hill type) |
type |
Evenness type |
coords |
Coordinates if provided |
n_sites |
Number of sites |
n_species |
Number of species |
Jost, L. (2010). The relation between evenness and diversity. Diversity, 2, 207-232.
Pielou, E.C. (1966). The measurement of diversity in different types of biological collections. Journal of Theoretical Biology, 13, 131-144.
diversityProfile() for Hill number profiles,
alphaDiversity() for raw diversity values
species <- matrix(rpois(50 * 30, 2), nrow = 50) # Hill evenness profile ev <- evenness(species) print(ev) # Pielou's J ev_j <- evenness(species, type = "pielou") print(ev_j)species <- matrix(rpois(50 * 30, 2), nrow = 50) # Hill evenness profile ev <- evenness(species) print(ev) # Pielou's J ev_j <- evenness(species, type = "pielou") print(ev_j)
Fit an asymptotic model to estimate total species richness beyond the observed sampling effort.
extrapolate( object, model = c("michaelis-menten", "lomolino", "asymptotic", "weibull", "logistic", "evt"), ... )extrapolate( object, model = c("michaelis-menten", "lomolino", "asymptotic", "weibull", "logistic", "evt"), ... )
object |
A |
model |
Character. Model to fit: |
... |
Additional arguments passed to |
An object of class spacc_fit containing:
asymptote |
Estimated total species richness |
asymptote_ci |
Confidence interval for asymptote |
model |
Model name |
fit |
The nls fit object |
aic |
AIC of the model |
Lomolino, M.V. (2000). Ecology's most general, yet protean pattern: the species-area relationship. Journal of Biogeography, 27, 17-26.
Flather, C.H. (1996). Fitting species-accumulation functions and assessing regional land use impacts on avian diversity. Journal of Biogeography, 23, 155-168.
Borda-de-Agua, L., Whittaker, R.J., Cardoso, P., et al. (2025). Extreme value theory explains the topography and scaling of the species-area relationship. Nature Communications, 16, 5346.
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, coords) fit <- extrapolate(sac) print(fit)coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, coords) fit <- extrapolate(sac) print(fit)
Predict species richness at coverage levels beyond the empirical maximum, following the Chao et al. (2014) framework. Provides seamless interpolation and extrapolation as a function of sample coverage.
extrapolateCoverage(x, target_coverage = c(0.9, 0.95, 0.99), q = 0)extrapolateCoverage(x, target_coverage = c(0.9, 0.95, 0.99), q = 0)
x |
A |
target_coverage |
Numeric vector of target coverage levels (0 to 1).
Can exceed observed coverage for extrapolation. Default |
q |
Numeric. Diversity order for extrapolation: 0 (richness, default), 1 (Shannon), or 2 (Simpson). |
For targets within observed coverage, linear interpolation is used. For targets beyond observed coverage, asymptotic estimators are applied:
q = 0: Chao1 estimator: S_est = S_obs + f1^2 / (2 * f2), where f1/f2 are singleton/doubleton counts. Extrapolation via coverage deficit.
q = 1: Shannon extrapolation based on the Good-Turing frequency formula.
q = 2: Simpson extrapolation using the unbiased estimator.
An object of class spacc_coverage_ext containing:
richness |
Matrix of interpolated/extrapolated richness (n_seeds x n_targets) |
target_coverage |
Target coverage levels |
q |
Diversity order used |
observed_coverage |
Mean observed final coverage |
observed_richness |
Mean observed final richness |
Chao, A. & Jost, L. (2012). Coverage-based rarefaction and extrapolation: standardizing samples by completeness rather than size. Ecology, 93, 2533-2547.
Chao, A., Gotelli, N.J., Hsieh, T.C., et al. (2014). Rarefaction and extrapolation with Hill numbers: a framework for sampling and estimation in species diversity studies. Ecological Monographs, 84, 45-67.
spaccCoverage(), interpolateCoverage()
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rpois(50 * 30, 2), nrow = 50) cov <- spaccCoverage(species, coords) ext <- extrapolateCoverage(cov, target_coverage = c(0.95, 0.99)) print(ext)coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rpois(50 * 30, 2), nrow = 50) cov <- spaccCoverage(species, coords) ext <- extrapolateCoverage(cov, target_coverage = c(0.95, 0.99)) print(ext)
Compute Hill numbers for the pooled community across all sites.
gammaDiversity(x, q = c(0, 1, 2))gammaDiversity(x, q = c(0, 1, 2))
x |
A site-by-species matrix (abundance data). |
q |
Numeric vector. Orders of diversity. Default |
Gamma diversity represents regional (total) diversity across all sites. It is computed by pooling abundances across all sites and calculating Hill numbers on the combined community.
A named numeric vector with gamma diversity for each q.
Jost, L. (2007). Partitioning diversity into independent alpha and beta components. Ecology, 88, 2427-2439.
alphaDiversity() for local diversity, diversityPartition()
for full alpha-beta-gamma decomposition
species <- matrix(rpois(50 * 30, 2), nrow = 50) gammaDiversity(species, q = c(0, 1, 2))species <- matrix(rpois(50 * 30, 2), nrow = 50) gammaDiversity(species, q = c(0, 1, 2))
Estimate total species richness using the improved Chao1 estimator (Chiu et al. 2014). Uses singletons through quadrupletons (f1–f4) to reduce bias for small samples.
iChao1(x)iChao1(x)
x |
A site-by-species matrix (abundance data). Columns are pooled across sites. |
The improved Chao1 estimator adds a correction term using f3 and f4:
When , the estimator collapses to Chao1.
An object of class spacc_estimate with components:
Name of the estimator ("iChao1")
Estimated total richness
Standard error of the estimate
Lower 95 percent confidence bound
Upper 95 percent confidence bound
Observed species richness
List with f1, f2, f3, f4
Chiu, C.H., Wang, Y.T., Walther, B.A. & Chao, A. (2014). An improved nonparametric lower bound of species richness via a modified Good-Turing frequency formula. Biometrics, 70, 671-682.
chao1() for the standard estimator, iChao2() for
incidence-based version
species <- matrix(rpois(50 * 30, 2), nrow = 50) iChao1(species)species <- matrix(rpois(50 * 30, 2), nrow = 50) iChao1(species)
Estimate total species richness from incidence data using the improved Chao2 estimator (Chiu et al. 2014). Uses uniques through quadruplicates (Q1–Q4) to reduce bias for small samples.
iChao2(x)iChao2(x)
x |
A site-by-species matrix (presence/absence or abundance). Automatically binarized. |
The improved Chao2 estimator is the incidence-based analogue of iChao1:
When , the estimator collapses to Chao2.
An object of class spacc_estimate with components:
Name of the estimator ("iChao2")
Estimated total richness
Standard error of the estimate
Lower 95 percent confidence bound
Upper 95 percent confidence bound
Observed species richness
List with Q1, Q2, Q3, Q4, n_sites
Chiu, C.H., Wang, Y.T., Walther, B.A. & Chao, A. (2014). An improved nonparametric lower bound of species richness via a modified Good-Turing frequency formula. Biometrics, 70, 671-682.
chao2() for the standard estimator, iChao1() for
abundance-based version
species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) iChao2(species)species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) iChao2(species)
Estimate species richness at specified coverage levels by interpolation.
interpolateCoverage(x, target = c(0.9, 0.95, 0.99))interpolateCoverage(x, target = c(0.9, 0.95, 0.99))
x |
A |
target |
Numeric vector of target coverage levels (0 to 1).
Default |
A data.frame with columns for each target coverage level, showing interpolated richness for each seed.
Chao, A. & Jost, L. (2012). Coverage-based rarefaction and extrapolation: standardizing samples by completeness rather than size. Ecology, 93, 2533-2547.
First- or second-order jackknife estimator of total species richness (Burnham & Overton 1978, 1979).
jackknife(x, order = 1L)jackknife(x, order = 1L)
x |
A site-by-species matrix (presence/absence or abundance). Automatically binarized. |
order |
Integer. Jackknife order: 1 (default) or 2. |
First-order jackknife:
Second-order jackknife:
An object of class spacc_estimate.
Burnham, K.P. & Overton, W.S. (1978). Estimation of the size of a closed population when capture probabilities vary among animals. Biometrika, 65, 625-633.
Burnham, K.P. & Overton, W.S. (1979). Robust estimation of population size when capture probabilities vary among animals. Ecology, 60, 927-936.
species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) jackknife(species, order = 1) jackknife(species, order = 2)species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) jackknife(species, order = 1) jackknife(species, order = 2)
Compute the expected species accumulation curve using sample-based rarefaction (Mao Tau estimator). This is analytically identical to the expected curve from random permutations.
mao_tau(x)mao_tau(x)
x |
A site-by-species matrix (presence/absence or abundance). |
A data.frame with columns: sites, expected, sd, lower, upper
Colwell, R.K., Mao, C.X. & Chang, J. (2004). Interpolating, extrapolating, and comparing incidence-based species accumulation curves. Ecology, 85, 2717-2727.
Create a ggplot2 visualization of species accumulation curves.
For grouped spacc objects, curves are overlaid with different colors.
## S3 method for class 'spacc' plot( x, ci = TRUE, ci_level = 0.95, ci_alpha = 0.3, show_seeds = FALSE, saturation = FALSE, saturation_level = 0.9, facet = FALSE, ... )## S3 method for class 'spacc' plot( x, ci = TRUE, ci_level = 0.95, ci_alpha = 0.3, show_seeds = FALSE, saturation = FALSE, saturation_level = 0.9, facet = FALSE, ... )
x |
A |
ci |
Logical. Show confidence interval ribbon? Default |
ci_level |
Numeric. Confidence level for interval. Default 0.95. |
ci_alpha |
Numeric. Transparency of CI ribbon. Default 0.3. |
show_seeds |
Logical. Show individual seed curves? Default |
saturation |
Logical. Mark saturation point? Default |
saturation_level |
Numeric. Proportion for saturation. Default 0.9. |
facet |
Logical. Use faceted panels for grouped objects? Default |
... |
Additional arguments (ignored). |
A ggplot2 object.
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, coords) plot(sac)coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, coords) plot(sac)
Create visualizations of per-site accumulation metrics.
## S3 method for class 'spacc_metrics' plot( x, metric = NULL, type = c("heatmap", "points", "histogram"), point_size = 3, palette = "viridis", ... )## S3 method for class 'spacc_metrics' plot( x, metric = NULL, type = c("heatmap", "points", "histogram"), point_size = 3, palette = "viridis", ... )
x |
A |
metric |
Character. Which metric to plot. Default is first metric. |
type |
Character. Plot type:
|
point_size |
Numeric. Size of points in heatmap. Default 3. |
palette |
Character. Color palette for heatmap. One of |
... |
Additional arguments (unused). |
A ggplot2 object.
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) metrics <- spaccMetrics(species, coords, metrics = c("slope_10", "auc")) plot(metrics, metric = "slope_10", type = "heatmap")coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) metrics <- spaccMetrics(species, coords, metrics = c("slope_10", "auc")) plot(metrics, metric = "slope_10", type = "heatmap")
Interpolate the mean empirical accumulation curve at arbitrary site counts
using linear interpolation. Unlike the predict method for spacc_fit
objects, this does not use a fitted model; it interpolates the observed
curve directly.
## S3 method for class 'spacc' predict(object, n = NULL, ci = TRUE, ci_level = 0.95, warn = TRUE, ...)## S3 method for class 'spacc' predict(object, n = NULL, ci = TRUE, ci_level = 0.95, warn = TRUE, ...)
object |
A |
n |
Numeric vector of site counts at which to interpolate. Defaults to 25%, 50%, and 100% of total sites. |
ci |
Logical. If |
ci_level |
Confidence level for the interval (default 0.95). |
warn |
Logical. Warn when |
... |
Ignored. |
A data frame (if ci = TRUE) or named numeric vector (if ci = FALSE).
Out-of-range values return NA.
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, coords, n_seeds = 10, progress = FALSE) predict(sac, n = c(10, 25, 50)) predict(sac, n = c(10, 25), ci = FALSE)coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, coords, n_seeds = 10, progress = FALSE) predict(sac, n = c(10, 25, 50)) predict(sac, n = c(10, 25), ci = FALSE)
Compute individual-based rarefaction curves for Hill numbers at any order q.
This complements the sample-based accumulation in spacc().
rarefy(x, n_individuals = NULL, q = 0, n_boot = 100L)rarefy(x, n_individuals = NULL, q = 0, n_boot = 100L)
x |
A site-by-species matrix with abundance data (not presence/absence). |
n_individuals |
Integer vector. Sample sizes to compute expected
diversity for. Default |
q |
Numeric. Order of Hill number; any value |
n_boot |
Integer. Number of bootstrap replicates for CI. Default 100. |
For q=0 (species richness): uses the Hurlbert (1971) formula.
For q=1 (Shannon diversity): rarefied Shannon entropy is computed and converted to effective number of species via exponentiation.
For q=2 (Simpson diversity): rarefied Simpson concentration is computed and converted to effective number of species via inversion.
An object of class spacc_rare containing:
n |
Sample sizes |
expected |
Expected diversity (Hill number of order q) |
sd |
Standard deviation |
lower, upper
|
95 percent confidence bounds |
q |
Order of diversity used |
Hurlbert, S.H. (1971). The nonconcept of species diversity: a critique and alternative parameters. Ecology, 52, 577-586.
Chao, A., Gotelli, N.J., Hsieh, T.C., et al. (2014). Rarefaction and extrapolation with Hill numbers: a framework for sampling and estimation in species diversity studies. Ecological Monographs, 84, 45-67.
abundance_matrix <- matrix(rpois(50 * 30, 2), nrow = 50) rare <- rarefy(abundance_matrix) plot(rare) # Shannon rarefaction rare_q1 <- rarefy(abundance_matrix, q = 1) plot(rare_q1)abundance_matrix <- matrix(rpois(50 * 30, 2), nrow = 50) rare <- rarefy(abundance_matrix) plot(rare) # Shannon rarefaction rare_q1 <- rarefy(abundance_matrix, q = 1) plot(rare_q1)
Compute standardized effect sizes by comparing observed diversity metrics against a null distribution from randomized communities. Supports multiple null model algorithms and works with most spacc output classes.
ses( x, species, coords = NULL, metric = NULL, null_model = c("frequency", "richness", "both", "curveball", "torus", "spatial_swap"), n_perm = 999L, parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )ses( x, species, coords = NULL, metric = NULL, null_model = c("frequency", "richness", "both", "curveball", "torus", "spatial_swap"), n_perm = 999L, parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )
x |
A spacc output object. Supported classes: |
species |
A site-by-species matrix (required). The species matrix used
to produce |
coords |
Optional data.frame with columns |
metric |
Character or |
null_model |
Character. Null model algorithm:
|
n_perm |
Integer. Number of permutations. Default 999. |
parallel |
Logical. Use parallel processing for the underlying analysis?
Default |
n_cores |
Integer. Number of cores. Default |
progress |
Logical. Show progress? Default |
seed |
Integer. Random seed for reproducibility. |
SES is computed as:
A two-tailed p-value is calculated as the proportion of null values at least as extreme as the observed value:
where is the rank of the observed value among null values.
Null model algorithms:
"frequency": Tests whether species composition matters given observed
species frequencies
"richness": Tests whether species identity matters given observed site
richness
"both": Maintains both marginal totals; tests non-random species
co-occurrence patterns
"curveball": Efficient alternative to "both" with proven uniform
sampling properties
An object of class spacc_ses containing:
observed |
Numeric vector of observed metric values |
null_mean |
Mean of null distribution |
null_sd |
Standard deviation of null distribution |
ses |
Standardized effect size: (observed - null_mean) / null_sd |
p_value |
Two-tailed p-value |
n_perm |
Number of permutations |
null_model |
Null model algorithm used |
metric |
Metric name |
input_class |
Class of input object |
Gotelli, N.J. (2000). Null model analysis of species co-occurrence patterns. Ecology, 81, 2606-2621.
Strona, G., Nappo, D., Boccacci, F., Fattorini, S. & San-Miguel-Ayanz, J. (2014). A fast and unbiased procedure to randomize ecological binary matrices with fixed row and column totals. Nature Communications, 5, 4114.
spaccHill(), spaccBeta(), spaccMetrics()
coords <- data.frame(x = runif(20), y = runif(20)) species <- matrix(rbinom(20 * 15, 1, 0.3), nrow = 20) sac <- spacc(species, coords, n_seeds = 10) result <- ses(sac, species, n_perm = 19) print(result)coords <- data.frame(x = runif(20), y = runif(20)) species <- matrix(rbinom(20 * 15, 1, 0.3), nrow = 20) sac <- spacc(species, coords, n_seeds = 10) result <- ses(sac, species, n_perm = 19) print(result)
Model the joint effect of sampling effort and area on species richness. Corrects for unequal survey intensity across sites, common in atlas data and citizen science datasets.
sesars(object, effort, model = c("power", "additive"), ...)sesars(object, effort, model = c("power", "additive"), ...)
object |
A |
effort |
Numeric vector. Sampling effort per site (e.g., hours, visits, trap-nights). Must have length equal to number of sites. |
model |
Character. SESARS model:
|
... |
Additional arguments passed to |
Standard SARs assume complete sampling within each area unit. SESARS incorporates sampling effort (E) alongside area (A) to provide unbiased richness estimates across regions with unequal survey intensity.
An object of class spacc_sesars containing:
model |
Model type |
fit |
Fitted model object |
coef |
Model coefficients |
data |
Data frame used for fitting |
Dennstadt, F., Horak, J. & Martin, M.D. (2019). Predictive sampling effort and species-area relationship models for estimating richness in fragmented landscapes. Diversity and Distributions, 26, 1112-1123.
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, coords) effort <- rpois(50, 10) ses <- sesars(sac, effort, model = "power") print(ses)coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, coords) effort <- rpois(50, 10) ses <- sesars(sac, effort, model = "power") print(ses)
Separate the effects of habitat loss (area reduction) from fragmentation (splitting into patches) on species richness. Extends the classic power-law SAR with an explicit fragmentation term.
sfar(object, patches, model = c("power", "log"), ...)sfar(object, patches, model = c("power", "log"), ...)
object |
A |
patches |
Factor or integer vector assigning each site to a habitat fragment (patch). Must have length equal to the number of sites. |
model |
Character. SFAR model:
|
... |
Additional arguments. |
The SFAR (Hanski et al. 2013) extends the power-law SAR to quantify the additional effect of habitat fragmentation on species richness. The model S = c * A^z * n^(-f) adds a penalty term for fragmentation (n = number of fragments), where f > 0 indicates that fragmentation reduces richness beyond what area loss alone would predict.
An object of class spacc_sfar containing:
fit |
Fitted model object |
coef |
Coefficients: c (intercept), z (area exponent), f (fragmentation exponent) |
n_patches |
Number of habitat fragments |
Hanski, I., Zurita, G.A., Bellocq, M.I. & Rybicki, J. (2013). Species-fragmented area relationship. Proceedings of the National Academy of Sciences, 110, 12715-12720.
Rybicki, J. & Hanski, I. (2013). Species-area relationships and extinctions caused by habitat loss and fragmentation. Ecology Letters, 16, 27-38.
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, coords) patches <- kmeans(coords, centers = 5)$cluster sfar_result <- sfar(sac, patches) print(sfar_result)coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) sac <- spacc(species, coords) patches <- kmeans(coords, centers = 5)$cluster sfar_result <- sfar(sac, patches) print(sfar_result)
Compute species accumulation curves using various spatial sampling methods with C++ backend for performance.
spacc( x, coords, n_seeds = 50L, method = c("knn", "kncn", "random", "radius", "gaussian", "cone", "collector"), distance = c("euclidean", "haversine"), backend = c("auto", "exact", "kdtree"), support = NULL, include_halo = TRUE, sigma = NULL, cone_width = pi/4, parallel = TRUE, n_cores = NULL, progress = TRUE, groups = NULL, time = NULL, w_space = 1, w_time = 1, seed = NULL, order = NULL )spacc( x, coords, n_seeds = 50L, method = c("knn", "kncn", "random", "radius", "gaussian", "cone", "collector"), distance = c("euclidean", "haversine"), backend = c("auto", "exact", "kdtree"), support = NULL, include_halo = TRUE, sigma = NULL, cone_width = pi/4, parallel = TRUE, n_cores = NULL, progress = TRUE, groups = NULL, time = NULL, w_space = 1, w_time = 1, seed = NULL, order = NULL )
x |
A site-by-species matrix (rows = sites, cols = species) with presence/absence (0/1) or abundance data. Can also be a data.frame. |
coords |
Site coordinates. Can be:
|
n_seeds |
Integer. Number of random starting points for uncertainty quantification. Default 50. |
method |
Character. Accumulation method:
|
distance |
Character. Distance method: |
backend |
Character. Nearest-neighbor backend for
|
support |
Optional. Spatial support for core/halo classification via
|
include_halo |
Logical. When |
sigma |
Numeric. Bandwidth for Gaussian method. Default auto-calculated. |
cone_width |
Numeric. Half-width in radians for cone method. Default pi/4. |
parallel |
Logical. Use parallel processing? Default |
n_cores |
Integer. Number of cores. Default |
progress |
Logical. Show progress bar? Default |
groups |
Optional. A factor, character, or integer vector of length
|
time |
Optional. Numeric vector of length |
w_space |
Numeric. Weight for spatial distance when |
w_time |
Numeric. Weight for temporal distance when |
seed |
Integer. Random seed for reproducibility. Default |
order |
Optional user-defined accumulation order(s). When supplied,
|
When groups = NULL, an object of class spacc containing:
curves |
Matrix of cumulative species counts (n_seeds x n_sites) |
coords |
Original coordinates |
n_seeds |
Number of seeds used |
method |
Method used |
n_species |
Total species in dataset |
Arrhenius, O. (1921). Species and area. Journal of Ecology, 9, 95-99.
Scheiner, S.M. (2003). Six types of species-area curves. Global Ecology and Biogeography, 12, 441-447.
Chiarucci, A., Bacaro, G., Scheiner, S.M. (2011). Old and new challenges in using species diversity for assessing biodiversity. Philosophical Transactions of the Royal Society B, 366, 2426-2437.
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) # Basic usage sac <- spacc(species, coords) plot(sac) # Different methods sac_knn <- spacc(species, coords, method = "knn") sac_rand <- spacc(species, coords, method = "random") comp <- compare(sac_knn, sac_rand)coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) # Basic usage sac <- spacc(species, coords) plot(sac) # Different methods sac_knn <- spacc(species, coords, method = "knn") sac_rand <- spacc(species, coords, method = "random") comp <- compare(sac_knn, sac_rand)
Analyze how beta diversity changes as sites are accumulated spatially. Partitions beta diversity into turnover (species replacement) and nestedness (species loss) components following Baselga (2010).
spaccBeta( x, coords, n_seeds = 50L, method = "knn", index = c("sorensen", "jaccard"), distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL, map = FALSE )spaccBeta( x, coords, n_seeds = 50L, method = "knn", index = c("sorensen", "jaccard"), distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL, map = FALSE )
x |
A site-by-species matrix (presence/absence or abundance). |
coords |
A data.frame with columns |
n_seeds |
Integer. Number of random starting points. Default 50. |
method |
Character. Accumulation method. Default |
index |
Character. Dissimilarity index: |
distance |
Character. Distance method: |
parallel |
Logical. Use parallel processing? Default |
n_cores |
Integer. Number of cores. |
progress |
Logical. Show progress? Default |
seed |
Integer. Random seed. |
map |
Logical. If |
At each step of spatial accumulation, beta diversity is calculated between the accumulated species pool and the newly added site. This reveals how species composition changes as you expand spatially.
Interpretation:
High turnover: New sites contribute different species (replacement)
High nestedness: New sites contribute subsets of existing species (loss)
The sum of turnover and nestedness equals total beta diversity.
An object of class spacc_beta containing:
beta_total |
Matrix of total beta diversity (n_seeds x n_sites-1) |
beta_turnover |
Matrix of turnover component |
beta_nestedness |
Matrix of nestedness component |
distance |
Matrix of cumulative distances |
n_seeds, n_sites, method, index
|
Parameters used |
Baselga, A. (2010). Partitioning the turnover and nestedness components of beta diversity. Global Ecology and Biogeography, 19, 134-143.
betapart::beta.pair() for pairwise beta diversity
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) beta <- spaccBeta(species, coords, n_seeds = 30) plot(beta) # Compare Sorensen vs Jaccard beta_jac <- spaccBeta(species, coords, index = "jaccard")coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) beta <- spaccBeta(species, coords, n_seeds = 30) plot(beta) # Compare Sorensen vs Jaccard beta_jac <- spaccBeta(species, coords, index = "jaccard")
Compute spatial accumulation of functional beta diversity, partitioned into turnover and nestedness components. Measures how functional trait space composition changes as sites are accumulated spatially.
spaccBetaFunc( x, coords, traits, n_seeds = 50L, method = "knn", index = c("sorensen", "jaccard"), distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )spaccBetaFunc( x, coords, traits, n_seeds = 50L, method = "knn", index = c("sorensen", "jaccard"), distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )
x |
A site-by-species matrix (presence/absence or abundance). |
coords |
A data.frame with columns |
traits |
A species-by-traits matrix. Row names should match species. |
n_seeds |
Integer. Number of random starting points. Default 50. |
method |
Character. Accumulation method. Default |
index |
Character. Dissimilarity index: |
distance |
Character. Distance method. Default |
parallel |
Logical. Use parallel processing? Default |
n_cores |
Integer. Number of cores. |
progress |
Logical. Show progress? Default |
seed |
Integer. Random seed. |
Functional beta diversity quantifies the turnover of functional traits across space. At each accumulation step, beta is computed based on the overlap of trait ranges (functional space) between the accumulated pool and the newly added site.
An object of class spacc_beta with additional attribute
beta_type = "functional".
Baselga, A. (2012). The relationship between species replacement, dissimilarity derived from nestedness, and nestedness. Global Ecology and Biogeography, 21, 1223-1232.
Cardoso, P., Rigal, F. & Carvalho, J.C. (2015). BAT – Biodiversity Assessment Tools. Methods in Ecology and Evolution, 6, 232-236.
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 20, 1, 0.3), nrow = 50) traits <- matrix(rnorm(20 * 3), nrow = 20) rownames(traits) <- colnames(species) <- paste0("sp", 1:20) beta_func <- spaccBetaFunc(species, coords, traits) plot(beta_func)coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 20, 1, 0.3), nrow = 50) traits <- matrix(rnorm(20 * 3), nrow = 20) rownames(traits) <- colnames(species) <- paste0("sp", 1:20) beta_func <- spaccBetaFunc(species, coords, traits) plot(beta_func)
Compute spatial accumulation of phylogenetic beta diversity, partitioned into turnover and nestedness components. Measures how evolutionary composition changes as sites are accumulated spatially.
spaccBetaPhylo( x, coords, tree, n_seeds = 50L, method = "knn", index = c("sorensen", "jaccard"), distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )spaccBetaPhylo( x, coords, tree, n_seeds = 50L, method = "knn", index = c("sorensen", "jaccard"), distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )
x |
A site-by-species matrix (presence/absence or abundance). |
coords |
A data.frame with columns |
tree |
A phylogenetic tree of class |
n_seeds |
Integer. Number of random starting points. Default 50. |
method |
Character. Accumulation method. Default |
index |
Character. Dissimilarity index: |
distance |
Character. Distance method. Default |
parallel |
Logical. Use parallel processing? Default |
n_cores |
Integer. Number of cores. |
progress |
Logical. Show progress? Default |
seed |
Integer. Random seed. |
Phylogenetic beta diversity quantifies evolutionary turnover across space. The PhyloSor index (phylogenetic Sorensen) is used: the fraction of branch length shared between two communities relative to total branch length. Partitioned into replacement (turnover) and loss (nestedness) components.
An object of class spacc_beta with additional attribute
beta_type = "phylogenetic".
Baselga, A. (2010). Partitioning the turnover and nestedness components of beta diversity. Global Ecology and Biogeography, 19, 134-143.
Chao, A., Chiu, C.H., Villeger, S., et al. (2023). Rarefaction and extrapolation with beta diversity under a framework of Hill numbers: the iNEXT.beta3D standardization. Ecological Monographs, 93, e1588.
if (requireNamespace("ape", quietly = TRUE)) { tree <- ape::rtree(20) coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 20, 1, 0.3), nrow = 50) colnames(species) <- tree$tip.label beta_phylo <- spaccBetaPhylo(species, coords, tree) plot(beta_phylo) }if (requireNamespace("ape", quietly = TRUE)) { tree <- ape::rtree(20) coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 20, 1, 0.3), nrow = 50) colnames(species) <- tree$tip.label beta_phylo <- spaccBetaPhylo(species, coords, tree) plot(beta_phylo) }
Compute spatial accumulation curves with sample coverage tracking. Allows standardization by completeness (coverage) rather than sample size, following Chao & Jost (2012) or the sample-based estimator of Chiu (2023).
spaccCoverage( x, coords, n_seeds = 50L, method = "knn", distance = c("euclidean", "haversine"), coverage = c("chao", "chiu"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL, map = FALSE )spaccCoverage( x, coords, n_seeds = 50L, method = "knn", distance = c("euclidean", "haversine"), coverage = c("chao", "chiu"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL, map = FALSE )
x |
A site-by-species matrix with abundance data. |
coords |
A data.frame with columns |
n_seeds |
Integer. Number of random starting points. Default 50. |
method |
Character. Accumulation method. Default |
distance |
Character. Distance method: |
coverage |
Character. Coverage estimator to use: |
parallel |
Logical. Use parallel processing? Default |
n_cores |
Integer. Number of cores. |
progress |
Logical. Show progress? Default |
seed |
Integer. Random seed. |
map |
Logical. If |
Sample coverage estimates the proportion of the total community abundance represented by observed species. It provides a measure of sampling completeness that is independent of sample size.
The Chao-Jost (2012) estimator counts singletons (f1) and doubletons (f2) in the cumulative abundance vector. It assumes individuals are sampled independently, which may not hold for plot-based spatial data.
The Chiu (2023) estimator uses incidence frequency counts instead: Q1 (species in exactly 1 site), Q2 (species in exactly 2 sites), and G1 (total abundance of Q1 species). It gives near-unbiased coverage estimates when organisms are spatially aggregated across sampling units.
Coverage-based rarefaction allows fair comparison of diversity across communities with different abundances by standardizing to equal completeness rather than equal sample size.
An object of class spacc_coverage containing:
richness |
Matrix of species richness (n_seeds x n_sites) |
individuals |
Matrix of individual counts |
coverage |
Matrix of coverage estimates |
coverage_type |
Coverage estimator used ( |
coords, n_seeds, n_sites, method
|
Parameters used |
Chao, A. & Jost, L. (2012). Coverage-based rarefaction and extrapolation: standardizing samples by completeness rather than size. Ecology, 93, 2533-2547.
Chiu, C.-H. (2023). A sample-based estimator for sample coverage. Ecology, 104, e4099.
iNEXT::iNEXT() for coverage-based rarefaction without spatial structure
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rpois(50 * 30, 2), nrow = 50) cov <- spaccCoverage(species, coords) plot(cov) # Sample-based coverage (recommended for spatial data) cov_chiu <- spaccCoverage(species, coords, coverage = "chiu") # Interpolate richness at 90% and 95% coverage interp <- interpolateCoverage(cov, target = c(0.90, 0.95))coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rpois(50 * 30, 2), nrow = 50) cov <- spaccCoverage(species, coords) plot(cov) # Sample-based coverage (recommended for spatial data) cov_chiu <- spaccCoverage(species, coords, coverage = "chiu") # Interpolate richness at 90% and 95% coverage interp <- interpolateCoverage(cov, target = c(0.90, 0.95))
Accumulate any user-supplied diversity index along a spatial ordering of
sites. At each accumulation step the cumulative community is passed to
fun, which returns a single number. This is the general escape hatch
behind the built-in metric functions: use it for indices that spacc does
not implement directly.
spaccDiversity( x, coords, fun, ..., method = c("knn", "kncn", "random", "radius", "collector"), incidence = FALSE, n_seeds = 50L, distance = c("euclidean", "haversine"), progress = TRUE, seed = NULL )spaccDiversity( x, coords, fun, ..., method = c("knn", "kncn", "random", "radius", "collector"), incidence = FALSE, n_seeds = 50L, distance = c("euclidean", "haversine"), progress = TRUE, seed = NULL )
x |
A site-by-species matrix (rows = sites, cols = species), abundance or presence/absence. |
coords |
A data.frame with columns |
fun |
A function applied to the cumulative community at each step. It
receives a named numeric vector of length |
... |
Additional arguments passed to |
method |
Character. Spatial ordering of sites: |
incidence |
Logical. If |
n_seeds |
Integer. Number of random starting points / orderings.
Ignored for |
distance |
Character. |
progress |
Logical. Show progress? Default |
seed |
Integer. Random seed for reproducibility. |
The site ordering reuses the same spatial traversals as the built-in
methods (nearest-neighbour, nearest-centroid, random, distance-rank, or
data order), then evaluates fun on the accumulating community. Because the
index is an arbitrary R function, this trades the speed of the compiled
metrics for full flexibility.
An object of class spacc_diversity that inherits from spacc, so
the standard summary(), plot(), as.data.frame() and predict()
methods apply. curves is an n_seeds x n_sites matrix of the metric
along the accumulation.
spaccHill(), spaccPhylo(), spaccFunc() for built-in metrics.
coords <- data.frame(x = runif(40), y = runif(40)) species <- matrix(rpois(40 * 20, 2), nrow = 40) # Shannon entropy along the accumulation shannon <- function(comm) { p <- comm[comm > 0] / sum(comm) -sum(p * log(p)) } div <- spaccDiversity(species, coords, shannon, n_seeds = 20) plot(div)coords <- data.frame(x = runif(40), y = runif(40)) species <- matrix(rpois(40 * 20, 2), nrow = 40) # Shannon entropy along the accumulation shannon <- function(comm) { p <- comm[comm > 0] / sum(comm) -sum(p * log(p)) } div <- spaccDiversity(species, coords, shannon, n_seeds = 20) plot(div)
Compute the number of endemic species (species found only within the accumulated area) as sites are added spatially. Complements the standard SAC by tracking species unique to each spatial extent.
spaccEndemism( x, coords, n_seeds = 50L, method = "knn", distance = c("euclidean", "haversine"), map = FALSE, parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )spaccEndemism( x, coords, n_seeds = 50L, method = "knn", distance = c("euclidean", "haversine"), map = FALSE, parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )
x |
A site-by-species matrix (presence/absence or abundance). |
coords |
A data.frame with columns |
n_seeds |
Integer. Number of random starting points. Default 50. |
method |
Character. Accumulation method. Default |
distance |
Character. Distance method: |
map |
Logical. If |
parallel |
Logical. Use parallel processing? Default |
n_cores |
Integer. Number of cores. |
progress |
Logical. Show progress? Default |
seed |
Integer. Random seed. |
At each accumulation step k, an endemic species is one that is present in the accumulated sites (1..k) but absent from all remaining unvisited sites (k+1..n). This tracks how many species are unique to the area sampled so far.
The endemism curve typically starts low (few endemics at small areas), increases as the region grows, and eventually equals total richness when all sites are included.
An object of class spacc_endemism containing:
richness |
Matrix of cumulative richness (n_seeds x n_sites) |
endemism |
Matrix of endemic species count (n_seeds x n_sites) |
site_values |
Per-site endemism data.frame (if |
coords, n_seeds, n_sites, method
|
Parameters used |
Kier, G., Kreft, H., Lee, T.M., et al. (2009). A global assessment of endemism and species richness across island and mainland regions. Proceedings of the National Academy of Sciences, 106, 9322-9327.
May, F., Gerstner, K., McGlinn, D.J., et al. (2018). mobsim: an R package for the simulation and measurement of biodiversity across spatial scales. Methods in Ecology and Evolution, 9, 1401-1408.
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) end <- spaccEndemism(species, coords, n_seeds = 30) plot(end)coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) end <- spaccEndemism(species, coords, n_seeds = 30) plot(end)
Compute spatial accumulation of functional diversity metrics based on traits.
spaccFunc( x, coords, traits, metric = c("fdis", "fric"), n_seeds = 50L, method = "knn", distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL, map = FALSE )spaccFunc( x, coords, traits, metric = c("fdis", "fric"), n_seeds = 50L, method = "knn", distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL, map = FALSE )
x |
A site-by-species matrix (abundance data recommended). |
coords |
A data.frame with columns |
traits |
A species-by-traits matrix. Row names should match species (columns of x). |
metric |
Character vector. Metrics to compute:
|
n_seeds |
Integer. Number of random starting points. Default 50. |
method |
Character. Accumulation method. Default |
distance |
Character. Site distance method: |
parallel |
Logical. Use parallel processing? Default |
n_cores |
Integer. Number of cores. |
progress |
Logical. Show progress? Default |
seed |
Integer. Random seed. |
map |
Logical. If |
Functional diversity metrics quantify trait space occupation:
FDis (Functional Dispersion): Abundance-weighted mean distance from the community centroid in trait space. Captures functional divergence.
FRic (Functional Richness): Volume of trait space occupied (convex hull). Requires more species than traits to compute.
An object of class spacc_func containing:
curves |
Named list of matrices, one per metric (n_seeds x n_sites) |
metric |
Metrics computed |
coords, n_seeds, n_sites, method
|
Parameters used |
Laliberté, E. & Legendre, P. (2010). A distance-based framework for measuring functional diversity from multiple traits. Ecology, 91, 299-305.
FD::dbFD() for comprehensive functional diversity analysis
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rpois(50 * 20, 2), nrow = 50) # Trait matrix (species x traits) traits <- matrix(rnorm(20 * 3), nrow = 20) rownames(traits) <- paste0("sp", 1:20) colnames(species) <- rownames(traits) func <- spaccFunc(species, coords, traits, metric = c("fdis", "fric")) plot(func)coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rpois(50 * 20, 2), nrow = 50) # Trait matrix (species x traits) traits <- matrix(rnorm(20 * 3), nrow = 20) rownames(traits) <- paste0("sp", 1:20) colnames(species) <- rownames(traits) func <- spaccFunc(species, coords, traits, metric = c("fdis", "fric")) plot(func)
Compute spatial species accumulation curves using Hill numbers (effective number of species) instead of raw richness. Hill numbers unify diversity measures: q=0 is richness, q=1 is exponential Shannon, q=2 is inverse Simpson.
spaccHill( x, coords, q = c(0, 1, 2), n_seeds = 50L, method = "knn", distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL, map = FALSE )spaccHill( x, coords, q = c(0, 1, 2), n_seeds = 50L, method = "knn", distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL, map = FALSE )
x |
A site-by-species matrix (rows = sites, cols = species) with presence/absence (0/1) or abundance data. |
coords |
A data.frame with columns |
q |
Numeric vector. Orders of diversity to compute. Default
|
n_seeds |
Integer. Number of random starting points. Default 50. |
method |
Character. Accumulation method: |
distance |
Character. Distance method: |
parallel |
Logical. Use parallel processing? Default |
n_cores |
Integer. Number of cores. Default |
progress |
Logical. Show progress bar? Default |
seed |
Integer. Random seed for reproducibility. |
map |
Logical. If |
Hill numbers (Chao et al. 2014) provide a unified framework for diversity measurement. Unlike raw richness (q=0), higher-order Hill numbers (q=1, q=2) down-weight rare species, providing different perspectives on diversity.
The spatial accumulation of Hill numbers can reveal scale-dependent diversity patterns missed by richness alone.
An object of class spacc_hill containing:
curves |
Named list of matrices, one per q value (n_seeds x n_sites) |
q |
Vector of q values used |
coords |
Original coordinates |
n_seeds |
Number of seeds |
n_sites |
Number of sites |
n_species |
Total species |
method |
Method used |
Chao, A., Gotelli, N.J., Hsieh, T.C., Sander, E.L., Ma, K.H., Colwell, R.K. & Ellison, A.M. (2014). Rarefaction and extrapolation with Hill numbers: a framework for sampling and estimation in species diversity studies. Ecological Monographs, 84, 45-67.
spacc() for richness-only accumulation, iNEXT::iNEXT() for
non-spatial Hill number rarefaction
# Compare diversity at different orders coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rpois(50 * 30, 2), nrow = 50) hill <- spaccHill(species, coords, q = c(0, 1, 2)) plot(hill) # Extract summary at final site summary(hill)# Compare diversity at different orders coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rpois(50 * 30, 2), nrow = 50) hill <- spaccHill(species, coords, q = c(0, 1, 2)) plot(hill) # Extract summary at final site summary(hill)
Compute multisite beta diversity as gamma/alpha decomposition of Hill numbers along the spatial accumulation curve (Jost 2007 framework).
spaccHillBeta( x, coords, q = c(0, 1, 2), n_seeds = 50L, distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )spaccHillBeta( x, coords, q = c(0, 1, 2), n_seeds = 50L, distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )
x |
A site-by-species matrix (rows = sites, cols = species). |
coords |
A data.frame with columns |
q |
Numeric vector. Orders of diversity. Default |
n_seeds |
Integer. Number of random starting points. Default 50. |
distance |
Character. |
parallel |
Logical. Use parallel processing? Default |
n_cores |
Integer. Number of cores. Default |
progress |
Logical. Show progress? Default |
seed |
Integer. Random seed. |
At each accumulation step k, the function computes:
Gamma: Hill number of the pooled community (all k sites combined)
Alpha: Generalized mean of per-site Hill numbers (Jost's power mean)
Beta: gamma / alpha (effective number of distinct communities)
Beta = 1 means all sites are identical; beta = k means all sites are
completely different. This provides the Hill-number analogue of the
Baselga-based spaccBeta().
An object of class spacc_hill_beta containing:
gamma |
Named list of n_seeds x n_sites matrices (one per q) |
alpha |
Named list of n_seeds x n_sites matrices (one per q) |
beta |
Named list of n_seeds x n_sites matrices (one per q) |
q |
Vector of q values |
coords |
Original coordinates |
n_seeds, n_sites, n_species
|
Dimensions |
Jost, L. (2007). Partitioning diversity into independent alpha and beta components. Ecology, 88, 2427-2439.
spaccBeta() for P/A-based Baselga partitioning,
spaccHill() for Hill accumulation without beta decomposition
coords <- data.frame(x = runif(40), y = runif(40)) species <- matrix(rpois(40 * 20, 2), nrow = 40) hb <- spaccHillBeta(species, coords, n_seeds = 10, progress = FALSE) plot(hb)coords <- data.frame(x = runif(40), y = runif(40)) species <- matrix(rpois(40 * 20, 2), nrow = 40) hb <- spaccHillBeta(species, coords, n_seeds = 10, progress = FALSE) plot(hb)
Compute spatial accumulation of Hill numbers alongside sample coverage, enabling standardized comparison at equal completeness levels (iNEXT-style analysis applied spatially).
spaccHillCoverage( x, coords, q = c(0, 1, 2), target_coverage = NULL, n_seeds = 50L, distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )spaccHillCoverage( x, coords, q = c(0, 1, 2), target_coverage = NULL, n_seeds = 50L, distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL )
x |
A site-by-species matrix (rows = sites, cols = species) with abundance data. |
coords |
A data.frame with columns |
q |
Numeric vector. Orders of diversity. Default |
target_coverage |
Numeric vector. Target coverage levels for
standardization. Default |
n_seeds |
Integer. Number of random starting points. Default 50. |
distance |
Character. Distance method: |
parallel |
Logical. Use parallel processing? Default |
n_cores |
Integer. Number of cores. Default |
progress |
Logical. Show progress? Default |
seed |
Integer. Random seed for reproducibility. |
Combines spatial kNN accumulation with simultaneous Hill number and coverage computation. At each accumulation step, both the Chao-Jost sample coverage and Hill numbers for all requested orders are calculated.
When target_coverage is specified, Hill numbers are interpolated at
those coverage levels using the existing interpolate_at_coverage() C++
function, enabling fair comparison between sites with different sampling
completeness.
An object of class spacc_hill_coverage containing:
hills |
Named list of n_seeds x n_sites matrices (one per q) |
coverage |
n_seeds x n_sites matrix of Chao-Jost coverage |
standardized |
List of per-q values interpolated at target coverage,
or |
q |
Vector of q values |
target_coverage |
Target coverage levels used |
coords |
Original coordinates |
n_seeds |
Number of seeds |
n_sites |
Number of sites |
n_species |
Total species |
Chao, A. & Jost, L. (2012). Coverage-based rarefaction and extrapolation: standardizing samples by completeness rather than size. Ecology, 93, 2533-2547.
Chao, A., Gotelli, N.J., Hsieh, T.C., Sander, E.L., Ma, K.H., Colwell, R.K. & Ellison, A.M. (2014). Rarefaction and extrapolation with Hill numbers. Ecological Monographs, 84, 45-67.
spaccHill() for Hill accumulation without coverage,
spaccCoverage() for coverage accumulation without Hill numbers
coords <- data.frame(x = runif(40), y = runif(40)) species <- matrix(rpois(40 * 20, 2), nrow = 40) hc <- spaccHillCoverage(species, coords, n_seeds = 10, progress = FALSE) plot(hc) # Standardize at 90% coverage hc2 <- spaccHillCoverage(species, coords, target_coverage = 0.9, n_seeds = 10, progress = FALSE) summary(hc2)coords <- data.frame(x = runif(40), y = runif(40)) species <- matrix(rpois(40 * 20, 2), nrow = 40) hc <- spaccHillCoverage(species, coords, n_seeds = 10, progress = FALSE) plot(hc) # Standardize at 90% coverage hc2 <- spaccHillCoverage(species, coords, target_coverage = 0.9, n_seeds = 10, progress = FALSE) summary(hc2)
Compute spatial accumulation metrics for each site as a starting point. Useful for identifying sites with high or low accumulation rates, visualizing spatial patterns in diversity, and understanding edge effects.
spaccMetrics( x, coords, metrics = c("slope_10", "half_richness", "auc"), method = c("knn", "kncn", "random"), distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE )spaccMetrics( x, coords, metrics = c("slope_10", "half_richness", "auc"), method = c("knn", "kncn", "random"), distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE )
x |
A site-by-species matrix (rows = sites, cols = species). |
coords |
A data.frame with columns |
metrics |
Character vector. Metrics to compute. Options include:
|
method |
Character. Accumulation method: |
distance |
Character. Distance method: |
parallel |
Logical. Use parallel processing? Default |
n_cores |
Integer. Number of cores for parallel processing. |
progress |
Logical. Show progress bar? Default |
This function runs a spatial accumulation curve starting from each site individually, then extracts summary metrics from each curve. This allows you to identify:
Sites in species-rich areas (high initial slope)
Core vs edge sites (fast vs slow accumulation)
Spatial patterns in community structure
The metrics can be plotted as a heatmap using plot(result, type = "heatmap"),
which requires the ggplot2 package. For more sophisticated spatial
visualization with study area boundaries, see the areaOfEffect package.
An object of class spacc_metrics containing:
metrics |
Data frame with one row per site and columns for each metric |
coords |
Original coordinates |
metric_names |
Names of computed metrics |
n_sites |
Number of sites |
n_species |
Total species count |
Soberon, J.M. & Llorente, J.B. (1993). The use of species accumulation functions for the prediction of species richness. Conservation Biology, 7, 480-488.
spacc() for standard accumulation curves
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) metrics <- spaccMetrics(species, coords, metrics = c("slope_10", "auc")) metrics$metrics$slope_10coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) metrics <- spaccMetrics(species, coords, metrics = c("slope_10", "auc")) metrics$metrics$slope_10
Compute spatial accumulation of phylogenetic diversity metrics (MPD, MNTD, PD).
spaccPhylo( x, coords, tree, metric = c("mpd", "mntd"), n_seeds = 50L, method = "knn", distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL, map = FALSE )spaccPhylo( x, coords, tree, metric = c("mpd", "mntd"), n_seeds = 50L, method = "knn", distance = c("euclidean", "haversine"), parallel = TRUE, n_cores = NULL, progress = TRUE, seed = NULL, map = FALSE )
x |
A site-by-species matrix. |
coords |
A data.frame with columns |
tree |
A phylogenetic tree of class |
metric |
Character vector. Metrics to compute:
|
n_seeds |
Integer. Number of random starting points. Default 50. |
method |
Character. Accumulation method. Default |
distance |
Character. Site distance method: |
parallel |
Logical. Use parallel processing? Default |
n_cores |
Integer. Number of cores. |
progress |
Logical. Show progress? Default |
seed |
Integer. Random seed. |
map |
Logical. If |
Phylogenetic diversity metrics incorporate evolutionary relationships:
MPD (Mean Pairwise Distance): Average phylogenetic distance between all pairs of species. Sensitive to tree-wide patterns.
MNTD (Mean Nearest Taxon Distance): Average distance to closest relative. Sensitive to terminal clustering.
PD (Faith's Phylogenetic Diversity): Total branch length connecting species. Requires full tree object.
An object of class spacc_phylo containing:
curves |
Named list of matrices, one per metric (n_seeds x n_sites) |
metric |
Metrics computed |
coords, n_seeds, n_sites, method
|
Parameters used |
Faith, D.P. (1992). Conservation evaluation and phylogenetic diversity. Biological Conservation, 61, 1-10.
Webb, C.O. (2000). Exploring the phylogenetic structure of ecological communities: an example for rain forest trees. American Naturalist, 156, 145-155.
picante::mpd(), picante::mntd(), picante::pd()
if (requireNamespace("ape", quietly = TRUE)) { # Create random tree tree <- ape::rtree(30) coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) colnames(species) <- tree$tip.label phylo <- spaccPhylo(species, coords, tree, metric = c("mpd", "mntd")) plot(phylo) }if (requireNamespace("ape", quietly = TRUE)) { # Create random tree tree <- ape::rtree(30) coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) colnames(species) <- tree$tip.label phylo <- spaccPhylo(species, coords, tree, metric = c("mpd", "mntd")) plot(phylo) }
Compute spatial eigenvectors from site coordinates using Principal Coordinates of Neighbour Matrices (PCNM) or distance-based Moran's Eigenvector Maps (dbMEM).
spatialEigenvectors( coords, threshold = NULL, method = c("pcnm", "dbmem"), distance = c("euclidean", "haversine") )spatialEigenvectors( coords, threshold = NULL, method = c("pcnm", "dbmem"), distance = c("euclidean", "haversine") )
coords |
A data.frame with columns |
threshold |
Numeric. Distance threshold for truncation. Default |
method |
Character. |
distance |
Character. |
PCNM (Borcard & Legendre 2002) decomposes spatial structure into orthogonal eigenvectors representing patterns at different spatial scales. Large eigenvalues correspond to broad-scale patterns (positive spatial autocorrelation), while small eigenvalues represent fine-scale patterns.
The algorithm:
Compute pairwise distances
Truncate distances beyond threshold (set to 4 * threshold)
Double-centre the squared truncated distance matrix
Extract eigenvectors with positive eigenvalues
For "dbmem": additionally filter to Moran's I > 0
An object of class spacc_mem containing:
vectors |
Matrix of eigenvectors (sites x n_vectors) |
eigenvalues |
Positive eigenvalues |
moran_i |
Moran's I for each eigenvector |
threshold |
Truncation distance used |
coords |
Original coordinates |
n_sites |
Number of sites |
method |
Method used |
Borcard, D. & Legendre, P. (2002). All-scale spatial analysis of ecological data by means of principal coordinates of neighbour matrices. Ecological Modelling, 153, 51-68.
Dray, S., Legendre, P. & Peres-Neto, P.R. (2006). Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM). Ecological Modelling, 196, 483-493.
spatialPartition() for variance partitioning with MEMs
coords <- data.frame(x = runif(30), y = runif(30)) mem <- spatialEigenvectors(coords) print(mem)coords <- data.frame(x = runif(30), y = runif(30)) mem <- spatialEigenvectors(coords) print(mem)
Partition diversity variation into spatial and non-spatial components using forward selection of Moran's Eigenvector Maps.
spatialPartition(x, mem, metric = NULL, forward = TRUE, alpha = 0.05)spatialPartition(x, mem, metric = NULL, forward = TRUE, alpha = 0.05)
x |
A spacc result object (e.g., |
mem |
A |
metric |
Character. Which metric to extract from |
forward |
Logical. Use forward selection? Default |
alpha |
Numeric. Significance threshold for forward selection. Default 0.05. |
Forward selection of MEMs proceeds by adding the MEM that most improves the model AIC at each step, stopping when no MEM improves AIC by more than 2 units or when p > alpha.
An object of class spacc_mem_partition containing:
r_squared_spatial |
R-squared of the spatial model |
r_squared_total |
Total R-squared (same as spatial here) |
selected_mems |
Names of selected MEM vectors |
n_selected |
Number of selected MEMs |
anova_table |
ANOVA table from the final model |
coefficients |
Model coefficients |
spatialEigenvectors() for computing MEMs
coords <- data.frame(x = runif(30), y = runif(30)) species <- matrix(rpois(30 * 15, 2), nrow = 30) mem <- spatialEigenvectors(coords) alpha <- alphaDiversity(species, q = 0) part <- spatialPartition(alpha$q0, mem) print(part)coords <- data.frame(x = runif(30), y = runif(30)) species <- matrix(rpois(30 * 15, 2), nrow = 30) mem <- spatialEigenvectors(coords) alpha <- alphaDiversity(species, q = 0) part <- spatialPartition(alpha$q0, mem) print(part)
Compute expected species richness accounting for spatial autocorrelation (Chiarucci et al. 2009). Uses distance-weighted sampling probabilities.
spatialRarefaction(x, coords, n_perm = 100, bandwidth = NULL)spatialRarefaction(x, coords, n_perm = 100, bandwidth = NULL)
x |
A site-by-species matrix. |
coords |
A data.frame with x and y columns. |
n_perm |
Number of permutations. Default 100. |
bandwidth |
Distance bandwidth for spatial weighting. |
A data.frame with columns: sites, mean, sd, lower, upper
Chiarucci, A., Bacaro, G., Rocchini, D. & Fattorini, L. (2009). Discovering and rediscovering the sample-based rarefaction formula in the ecological literature. Community Ecology, 10, 195-199.
Reduce spatial autocorrelation by subsampling sites using various methods.
subsample( coords, n = NULL, method = c("grid", "random", "thinning"), cell_size = NULL, min_dist = NULL, seed = NULL )subsample( coords, n = NULL, method = c("grid", "random", "thinning"), cell_size = NULL, min_dist = NULL, seed = NULL )
coords |
A data.frame with columns |
n |
Integer. Target number of sites to retain. If |
method |
Character. Subsampling method: |
cell_size |
Numeric. Grid cell size for |
min_dist |
Numeric. Minimum distance between retained sites for
|
seed |
Integer. Random seed for reproducibility. |
Methods:
"grid": Overlay a grid and retain one random site per cell.
"random": Simple random subsample of n sites.
"thinning": Iteratively remove sites until minimum distance is achieved.
Integer vector of row indices to retain.
Aiello-Lammens, M.E., Boria, R.A., Radosavljevic, A., et al. (2015). spThin: an R package for spatial thinning of species occurrence records for use in ecological niche models. Ecography, 38, 541-545.
Lennon, J.J., Koleff, P., Greenwood, J.J.D. & Gaston, K.J. (2004). Contribution of rarity and commonness to patterns of species richness. Ecology Letters, 7, 81-87.
coords <- data.frame(x = runif(50) * 100, y = runif(50) * 100) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) keep <- subsample(coords, method = "grid", cell_size = 20) sac <- spacc(species[keep, ], coords[keep, ])coords <- data.frame(x = runif(50) * 100, y = runif(50) * 100) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) keep <- subsample(coords, method = "grid", cell_size = 20) sac <- spacc(species[keep, ], coords[keep, ])
Accumulate species within an expanding radius from seed points. Models invasion spread from introduction points.
wavefront( x, coords, n_seeds = 50L, r0 = 0, dr = NULL, n_steps = 50L, distance = c("euclidean", "haversine"), progress = TRUE, seed = NULL )wavefront( x, coords, n_seeds = 50L, r0 = 0, dr = NULL, n_steps = 50L, distance = c("euclidean", "haversine"), progress = TRUE, seed = NULL )
x |
A site-by-species matrix. |
coords |
A data.frame with x and y columns, or a |
n_seeds |
Integer. Number of random starting points. |
r0 |
Numeric. Initial radius. Default 0. |
dr |
Numeric. Radius increment per step. Default auto-calculated. |
n_steps |
Integer. Number of expansion steps. Default 50. |
distance |
Character. Distance method. |
progress |
Logical. Show progress? |
seed |
Integer. Random seed. |
An object of class spacc_wavefront containing:
curves |
Matrix of species counts (n_seeds x n_steps) |
radius |
Vector of radius values |
sites_included |
Matrix of sites included at each step |
Shigesada, N. & Kawasaki, K. (1997). Biological Invasions: Theory and Practice. Oxford University Press.
coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) wf <- wavefront(species, coords, n_seeds = 20, n_steps = 50) plot(wf)coords <- data.frame(x = runif(50), y = runif(50)) species <- matrix(rbinom(50 * 30, 1, 0.3), nrow = 50) wf <- wavefront(species, coords, n_seeds = 20, n_steps = 50) plot(wf)
Compute zeta diversity — the mean number of species shared across k sites — for increasing orders of k. The zeta decline curve reveals community assembly processes: exponential decline suggests stochastic assembly, while power-law decline indicates niche-based assembly.
zetaDiversity( x, coords, orders = 1:10, n_samples = 100L, method = c("knn", "random"), distance = c("euclidean", "haversine"), seed = NULL, progress = TRUE )zetaDiversity( x, coords, orders = 1:10, n_samples = 100L, method = c("knn", "random"), distance = c("euclidean", "haversine"), seed = NULL, progress = TRUE )
x |
A site-by-species matrix (presence/absence or abundance). Automatically binarized. |
coords |
A data.frame with columns |
orders |
Integer vector. Orders of zeta diversity to compute (number of
sites in each combination). Default |
n_samples |
Integer. Number of random combinations to sample per order. Default 100. |
method |
Character. Method for selecting k-site combinations:
|
distance |
Character. Distance method: |
seed |
Integer. Random seed for reproducibility. Default |
progress |
Logical. Show progress? Default |
Zeta diversity of order k () is the mean number of species
shared across k sites. Key properties:
= mean species richness per site
= mean number of species shared by any two sites
decreases monotonically with k
The zeta decline ratio () is diagnostic:
Constant ratio: exponential decline (stochastic assembly)
Increasing ratio: power-law decline (deterministic/niche-based assembly)
The knn method selects spatially nearest k sites from each focal site,
which is ecologically meaningful for testing spatial turnover. The random
method samples random k-site combinations, providing a null expectation.
An object of class spacc_zeta containing:
zeta |
Mean zeta values per order |
zeta_sd |
Standard deviations per order |
orders |
The k values |
n_samples |
Number of samples per order |
ratio |
Zeta ratio: zeta_k / zeta_(k-1) |
decline |
Data.frame with exponential and power-law fit statistics |
method |
Method used |
n_sites |
Number of sites |
n_species |
Total species count |
Hui, C. & McGeoch, M.A. (2014). Zeta diversity as a concept and metric that unifies incidence-based biodiversity patterns. The American Naturalist, 184, 684-694.
Latombe, G., McGeoch, M.A., Nipperess, D.A. & Hui, C. (2018). zetadiv: an R package for computing compositional change across multiple sites, assemblages or cases. bioRxiv, 324897.
spaccBeta() for pairwise beta diversity, distanceDecay() for
distance-decay relationships
coords <- data.frame(x = runif(30), y = runif(30)) species <- matrix(rbinom(30 * 20, 1, 0.3), nrow = 30) zeta <- zetaDiversity(species, coords, orders = 1:5, n_samples = 50) plot(zeta)coords <- data.frame(x = runif(30), y = runif(30)) species <- matrix(rbinom(30 * 20, 1, 0.3), nrow = 30) zeta <- zetaDiversity(species, coords, orders = 1:5, n_samples = 50) plot(zeta)