---
title: "Statistical Applications of taxodist"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Statistical Applications of taxodist}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center"
)

library(taxodist)

vegan_available <- requireNamespace("vegan", quietly = TRUE)
ape_available <- requireNamespace("ape", quietly = TRUE)

if (is.null(taxobase$statistical_matrix)) {
  stop(
    "The installed taxobase object does not contain statistical_matrix. ",
    "Rebuild data/taxobase.rda before building this vignette."
  )
}
```

## 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.

```{r data}
taxa <- taxobase$statistical_taxa
mat <- taxobase$statistical_matrix

taxobase$metadata
length(taxa)
inherits(mat, "dist")
identical(attr(mat, "Labels"), taxa)
```

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.

```{r taxa}
taxa
```

## Inspecting the distance matrix

```{r matrix-summary}
summary(as.vector(mat))
range(mat)

round(
  as.matrix(mat)[1:6, 1:6],
  digits = 5
)
```

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

```{r clustering, fig.width=7, fig.height=5}
clustering <- taxo_cluster(mat, method = "average")

clustering$hclust
summary(clustering$hclust$height)

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.

```{r ape-tree, eval=ape_available, fig.width=7, fig.height=5}
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:

```{r ape-export, eval=FALSE}
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).

```{r pcoa}
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
100 * ordination$GOF[1]
sum(ordination$eig < -sqrt(.Machine$double.eps))
```

```{r pcoa-plot, fig.width=7, fig.height=5}
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.

```{r community}
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
```

```{r taxondive, eval=vegan_available}
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
```

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.

```{r mantel, eval=vegan_available}
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
```

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.

```{r groups}
groups <- factor(c(
  "xenarthran", "xenarthran",
  "carnivoran",
  "ungulate",
  "marsupial",
  "primate", "primate",
  "carnivoran", "carnivoran", "carnivoran",
  "cetacean", "cetacean",
  "sirenian",
  "ungulate", "ungulate"
))

table(groups)
```

```{r permanova, eval=vegan_available}
set.seed(42)

permanova_result <- vegan::adonis2(
  mat ~ groups,
  permutations = 999
)

permanova_result
```

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.

```{r dispersion, eval=vegan_available}
dispersion <- vegan::betadisper(mat, groups)

anova(dispersion)

set.seed(42)
vegan::permutest(
  dispersion,
  permutations = 999
)
```

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.

```{r live-workflow, eval=FALSE}
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:

```{r citations, eval=FALSE}
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.