Statistical Applications of taxodist
Overview
A distance matrix produced by taxodist is a base R
dist object. It can therefore be supplied to clustering,
ordination, permutation tests, and taxonomic-diversity functions that
accept pairwise dissimilarities.
These analyses operate on distances derived from classification
depth. They do not convert the source hierarchy into evolutionary time,
genetic distance, or phylogenetic branch length. Results should be
described as analyses of taxonomic hierarchy distance.
This vignette uses the reference matrix stored in
taxobase. The matrix is built in advance from The
Taxonomicon and permits the vignette to run without network access.
taxa <- taxobase$statistical_taxa
mat <- taxobase$statistical_matrix
taxobase$metadata
#> $source
#> [1] "The Taxonomicon"
#>
#> $source_url
#> [1] "http://taxonomicon.taxonomy.nl"
#>
#> $generated_on
#> [1] "2026-08-06"
#>
#> $package_version
#> [1] "0.6.0"
#>
#> $distance_definition
#> [1] "0 for identical nodes; otherwise 1 / depth(MRCA)"
length(taxa)
#> [1] 15
inherits(mat, "dist")
#> [1] TRUE
identical(attr(mat, "Labels"), taxa)
#> [1] TRUE
The example contains 15 mammalian genera representing several broad
taxonomic groups. The selection is illustrative and is not intended to
constitute a formal ecological sample or a conservation-status
dataset.
taxa
#> [1] "Priodontes" "Myrmecophaga" "Chrysocyon" "Tapirus"
#> [5] "Didelphis" "Leontopithecus" "Brachyteles" "Panthera"
#> [9] "Pteronura" "Puma" "Sotalia" "Pontoporia"
#> [13] "Trichechus" "Mazama" "Blastocerus"
Inspecting the distance matrix
summary(as.vector(mat))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.01370 0.01538 0.01613 0.01593 0.01639 0.01695
range(mat)
#> [1] 0.01369863 0.01694915
round(
as.matrix(mat)[1:6, 1:6],
digits = 5
)
#> Priodontes Myrmecophaga Chrysocyon Tapirus Didelphis
#> Priodontes 0.00000 0.01587 0.01639 0.01639 0.01695
#> Myrmecophaga 0.01587 0.00000 0.01639 0.01639 0.01695
#> Chrysocyon 0.01639 0.01639 0.00000 0.01538 0.01695
#> Tapirus 0.01639 0.01639 0.01538 0.00000 0.01695
#> Didelphis 0.01695 0.01695 0.01695 0.01695 0.00000
#> Leontopithecus 0.01639 0.01639 0.01613 0.01613 0.01695
#> Leontopithecus
#> Priodontes 0.01639
#> Myrmecophaga 0.01639
#> Chrysocyon 0.01613
#> Tapirus 0.01613
#> Didelphis 0.01695
#> Leontopithecus 0.00000
For this dataset, pairwise distances range from approximately 0.01370
to 0.01695. Smaller values indicate a deeper MRCA in the source
hierarchy.
The absolute scale depends on classification resolution. A value such
as 0.014 should not be interpreted as a proportion of evolutionary
divergence.
Hierarchical clustering
clustering <- taxo_cluster(mat, method = "average")
clustering$hclust
#>
#> Call:
#> stats::hclust(d = d, method = method)
#>
#> Cluster method : average
#> Number of objects: 15
summary(clustering$hclust$height)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.01370 0.01429 0.01504 0.01509 0.01607 0.01695
plot(
clustering,
main = "Average-linkage clustering of taxonomic hierarchy distances",
xlab = "",
sub = ""
)
Average-linkage clustering summarizes similarities in the supplied
matrix. The topology and node heights also depend on the selected
linkage method.
The resulting dendrogram is not an independently inferred
phylogenetic tree. It contains no substitution model, fossil
calibration, branch-length estimation, or statistical reconstruction of
evolutionary history.
Tree representation with ape
An hclust object can be converted to the
phylo data structure used by ape. This can be
useful for plotting and file-format interoperability.
tree <- ape::as.phylo(clustering$hclust)
plot(
tree,
main = "Tree representation of a taxonomic distance dendrogram",
cex = 0.8
)
Conversion to class phylo changes the representation of
the clustering result, not its scientific interpretation. The object
remains derived from a classification-based distance matrix.
Export can be performed when required:
ape::write.tree(tree, file = "taxonomic-distance-dendrogram.nwk")
Principal coordinates analysis
Principal coordinates analysis, or classical multidimensional
scaling, represents pairwise dissimilarities as coordinates in a
Euclidean space (Gower, 1966).
ordination <- taxo_ordinate(mat, k = 2)
positive_eigenvalues <- ordination$eig[ordination$eig > 0]
variance_percent <- 100 * positive_eigenvalues / sum(positive_eigenvalues)
ordination_summary <- data.frame(
Axis = c("PC1", "PC2"),
Eigenvalue = ordination$eig[1:2],
Variance_percent = variance_percent[1:2],
Cumulative_percent = cumsum(variance_percent)[1:2]
)
ordination_summary
#> Axis Eigenvalue Variance_percent Cumulative_percent
#> 1 PC1 0.0002040109 11.459192 11.45919
#> 2 PC2 0.0001613024 9.060276 20.51947
100 * ordination$GOF[1]
#> [1] 20.51947
sum(ordination$eig < -sqrt(.Machine$double.eps))
#> [1] 0
plot(
ordination,
main = "PCoA of taxonomic hierarchy distances"
)
The first two axes explain approximately 11.46% and 9.06% of the
positive eigenvalue total, respectively. Their cumulative representation
is therefore about 20.52%.
Consequently, the two-dimensional plot is a limited projection of the
complete distance structure. Proximity in this graph should not be
interpreted without considering its relatively low goodness-of-fit.
No materially negative eigenvalues were observed for this matrix.
This does not imply that every matrix produced from missing,
disconnected, or differently processed lineages will have the same
property.
Taxonomic diversity and distinctness
Clarke and Warwick (1998, 2001) developed indices based on averaged
taxonomic distances among species or individuals.
vegan::taxondive() accepts a community matrix together with
a dist object.
The following communities are hypothetical presence-absence
assemblages. They are constructed only to demonstrate the interface.
comm <- matrix(
c(
1,1,0,0,1,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,1,1,1,1,1
),
nrow = 3,
byrow = TRUE,
dimnames = list(
c("community_A", "community_B", "community_C"),
taxa
)
)
comm
#> Priodontes Myrmecophaga Chrysocyon Tapirus Didelphis Leontopithecus
#> community_A 1 1 0 0 1 0
#> community_B 0 0 0 0 0 1
#> community_C 0 0 0 0 0 0
#> Brachyteles Panthera Pteronura Puma Sotalia Pontoporia Trichechus
#> community_A 0 0 0 0 0 0 0
#> community_B 1 1 1 1 0 0 0
#> community_C 0 0 0 0 1 1 1
#> Mazama Blastocerus
#> community_A 0 0
#> community_B 0 0
#> community_C 1 1
taxonomic_diversity <- vegan::taxondive(comm, mat)
taxonomic_diversity_table <- data.frame(
Species = taxonomic_diversity$Species,
Delta = taxonomic_diversity$D,
Delta_star = taxonomic_diversity$Dstar,
Lambda_plus = taxonomic_diversity$Lambda,
Delta_plus = taxonomic_diversity$Dplus,
SD_Delta_plus = taxonomic_diversity$sd.Dplus
)
taxonomic_diversity_table
#> Species Delta Delta_star Lambda_plus Delta_plus SD_Delta_plus
#> community_A 3 0.01659044 0.01659044 2.573489e-07 0.01659044 0.0005404065
#> community_B 5 0.01539343 0.01539343 8.359778e-07 0.01539343 0.0003567428
#> community_C 5 0.01526725 0.01526725 1.046190e-06 0.01526725 0.0003567428
Here:
- \(\Delta\) is taxonomic diversity
and incorporates abundances;
- \(\Delta^*\) is taxonomic
distinctness conditional on individuals belonging to different
taxa;
- \(\Delta^+\) is average taxonomic
distinctness for presence-absence data;
- \(\Lambda^+\) describes variation
in pairwise taxonomic distinctness.
Because these examples use presence-absence data, \(\Delta\), \(\Delta^*\), and \(\Delta^+\) coincide within each
community.
community_A has the largest mean distance, approximately
0.01659. community_B and community_C have
values of approximately 0.01539 and 0.01527. These differences describe
only the deliberately constructed example communities.
Traditional applications often calculate taxonomic distances from
integer or weighted steps through standardized ranks. The distances
supplied here instead use reciprocal MRCA depth.
taxondive() permits other dissimilarities, but the
resulting values are specific to this definition and should not be
compared numerically with conventional Clarke-Warwick indices without
accounting for the different distance scale. Multiplying all distances
by a constant also rescales the Delta indices and rescales \(\Lambda^+\) by the square of that
constant.
Mantel test
The Mantel test evaluates association between two distance matrices
using permutations (Mantel, 1967). Both matrices must refer to the same
observations in the same order.
The coordinates below are simulated and have no biological
interpretation.
set.seed(42)
coordinates <- matrix(
rnorm(2 * length(taxa)),
ncol = 2,
dimnames = list(taxa, c("x", "y"))
)
geographic_distance <- stats::dist(coordinates)
mantel_result <- vegan::mantel(
mat,
geographic_distance,
method = "pearson",
permutations = 999
)
mantel_result
#>
#> Mantel statistic based on Pearson's product-moment correlation
#>
#> Call:
#> vegan::mantel(xdis = mat, ydis = geographic_distance, method = "pearson", permutations = 999)
#>
#> Mantel statistic r: -0.05619
#> Significance: 0.655
#>
#> Upper quantiles of permutations (null model):
#> 90% 95% 97.5% 99%
#> 0.187 0.234 0.268 0.336
#> Permutation: free
#> Number of permutations: 999
The simulated example gives a Mantel correlation of approximately
-0.056 and a permutation p-value of approximately 0.655. It therefore
provides no evidence of association between taxonomic and simulated
geographic distance.
This is the expected interpretation of this example: the coordinates
were generated randomly. It is not evidence for or against a biological
biogeographic process.
In an empirical analysis, the permutation design must reflect the
sampling design. Spatial dependence, repeated measurements, blocks, and
nested sampling can invalidate unrestricted permutations.
PERMANOVA
PERMANOVA partitions variation in a distance matrix and assesses a
pseudo-\(F\) statistic using
permutations (Anderson, 2001).
The groups below are illustrative taxonomic categories. They are
uneven in size, and some contain only one observation. They are
therefore suitable for demonstrating the software interface but not for
drawing substantive inferential conclusions.
groups <- factor(c(
"xenarthran", "xenarthran",
"carnivoran",
"ungulate",
"marsupial",
"primate", "primate",
"carnivoran", "carnivoran", "carnivoran",
"cetacean", "cetacean",
"sirenian",
"ungulate", "ungulate"
))
table(groups)
#> groups
#> carnivoran cetacean marsupial primate sirenian ungulate xenarthran
#> 4 2 1 2 1 3 2
set.seed(42)
permanova_result <- vegan::adonis2(
mat ~ groups,
permutations = 999
)
permanova_result
#> Permutation test for adonis under reduced model
#> Permutation: free
#> Number of permutations: 999
#>
#> vegan::adonis2(formula = mat ~ groups, permutations = 999)
#> Df SumOfSqs R2 F Pr(>F)
#> Model 6 0.00093360 0.5244 1.4701 0.001 ***
#> Residual 8 0.00084673 0.4756
#> Total 14 0.00178033 1.0000
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For these constructed groups, the model accounts for approximately
52.44% of the total sum of squares and gives a permutation p-value of
0.001.
That result must not immediately be interpreted as evidence that
group centroids differ. PERMANOVA can respond to differences in group
location, within-group dispersion, or both, particularly in an
unbalanced design.
Multivariate dispersion
Distance-based analysis of multivariate dispersion examines distances
from observations to their group centroid or spatial median (Anderson,
2006). It provides an important companion analysis for PERMANOVA.
dispersion <- vegan::betadisper(mat, groups)
anova(dispersion)
#> Analysis of Variance Table
#>
#> Response: Distances
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Groups 6 1.1906e-04 1.9844e-05 84.335 8.584e-07 ***
#> Residuals 8 1.8820e-06 2.3530e-07
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
set.seed(42)
vegan::permutest(
dispersion,
permutations = 999
)
#>
#> Permutation test for homogeneity of multivariate dispersions
#> Permutation: free
#> Number of permutations: 999
#>
#> Response: Distances
#> Df Sum Sq Mean Sq F N.Perm Pr(>F)
#> Groups 6 1.1906e-04 1.9844e-05 84.335 999 0.001 ***
#> Residuals 8 1.8820e-06 2.3530e-07
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The ordinary analysis of distances to group centers is highly
significant in this example. Thus the groups differ strongly in
dispersion.
The PERMANOVA result is consequently confounded with dispersion
differences and cannot be presented as an unambiguous difference among
group centroids. Moreover, singleton and very small groups make this
example unsuitable for formal biological inference.
A real study should use independently justified groups, adequate
replication, and a permutation structure consistent with the sampling
design.
Complete live workflow
The following workflow computes a new matrix from taxon names. It is
not evaluated while the vignette is built because it requires access to
The Taxonomicon.
library(taxodist)
taxa <- c(
"Tyrannosaurus",
"Velociraptor",
"Spinosaurus",
"Allosaurus"
)
coverage <- check_coverage(taxa)
stopifnot(all(coverage))
mat <- distance_matrix(taxa)
clustering <- taxo_cluster(mat, method = "average")
ordination <- taxo_ordinate(mat, k = 2)
plot(clustering)
plot(ordination)
Reporting recommendations
Analyses based on taxodist should report:
- the package version;
- the lineage retrieval or dataset-generation date;
- The Taxonomicon as the classification source;
- the distance definition;
- decisions involving ambiguous or missing taxa;
- the clustering linkage or ordination method;
- the number and structure of permutations;
- checks of multivariate dispersion when PERMANOVA is used;
- the distinction between hierarchy distance and evolutionary branch
length.
The exact package and data-source citations are available with:
citation("taxodist")
citation("vegan")
citation("ape")
References
Anderson, M.J. (2001). A new method for non-parametric multivariate
analysis of variance. Austral Ecology, 26, 32–46. doi:10.1046/j.1442-9993.2001.01070.x.
Anderson, M.J. (2006). Distance-based tests for homogeneity of
multivariate dispersions. Biometrics, 62, 245–253. doi:10.1111/j.1541-0420.2005.00440.x.
Brands, S.J. (1989 onwards). Systema Naturae 2000.
Amsterdam, The Netherlands. The Taxonomicon.
Clarke, K.R. and Warwick, R.M. (1998). A taxonomic distinctness index
and its statistical properties. Journal of Applied Ecology, 35,
523–531. doi:10.1046/j.1365-2664.1998.3540523.x.
Clarke, K.R. and Warwick, R.M. (2001). A further biodiversity index
applicable to species lists: variation in taxonomic distinctness.
Marine Ecology Progress Series, 216, 265–278. doi:10.3354/meps216265.
Gower, J.C. (1966). Some distance properties of latent root and
vector methods used in multivariate analysis. Biometrika, 53,
325–338. doi:10.1093/biomet/53.3-4.325.
Mantel, N. (1967). The detection of disease clustering and a
generalized regression approach. Cancer Research, 27,
209–220.
Paradis, E. and Schliep, K. (2019). ape 5.0: an environment for
modern phylogenetics and evolutionary analyses in R.
Bioinformatics, 35, 526–528. doi:10.1093/bioinformatics/bty633.