---
title: "Methodological Notes and Limitations"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Methodological Notes and Limitations}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(taxodist)
```

## Scope

`taxodist` represents each taxon as an ordered lineage from the root of The
Taxonomicon classification to the queried node. The resulting distance is a
property of that classification hierarchy. It is not an estimate of divergence
time, genetic distance, or phylogenetic branch length.

## Definition

Let \(L_A\) and \(L_B\) be two ordered lineages. Let \(h(A,B)\) be the length of
their continuous common prefix, equivalent to the depth of their most recent
common ancestor (MRCA). The distance is

\[
d(A,B) =
\begin{cases}
0, & L_A = L_B, \\
1/h(A,B), & L_A \ne L_B.
\end{cases}
\]

If two lineages have no shared root, the package returns an infinite distance.
Names repeated after the first divergence are not considered evidence of shared
ancestry.

Zero is reserved for the same hierarchy node. Consequently, an ancestor and
one of its descendants have positive distance even though they are connected by
ancestry. Membership and distance answer different questions: `is_member()` and
`taxo_path()` should be used for containment and ancestry queries.

## Ultrametric property

For lineages belonging to one connected hierarchy, the distance is
non-negative, symmetric, and satisfies identity of indiscernibles. It also
satisfies the strong triangle inequality

\[
d(A,C) \leq \max\{d(A,B), d(B,C)\}.
\]

The common prefix of \(A\) and \(C\) must be at least as long as the shorter of
the common prefixes of \((A,B)\) and \((B,C)\). Therefore,

\[
h(A,C) \geq \min\{h(A,B),h(B,C)\}.
\]

Taking reciprocals gives the strong triangle inequality.

The property can also be checked on the offline reference matrix distributed
with the package.

```{r ultrametric-check}
m <- as.matrix(taxobase$statistical_matrix)
tolerance <- sqrt(.Machine$double.eps)

ultrametric_ok <- function(index) {
  i <- index[1]
  j <- index[2]
  k <- index[3]

  d_ij <- m[i, j]
  d_ik <- m[i, k]
  d_jk <- m[j, k]

  d_ij <= max(d_ik, d_jk) + tolerance &&
    d_ik <= max(d_ij, d_jk) + tolerance &&
    d_jk <= max(d_ij, d_ik) + tolerance
}

all(combn(seq_len(nrow(m)), 3, FUN = ultrametric_ok))
```

## Alternative metrics

The reciprocal distance weights shallow divergences heavily: two taxa whose
MRCA is at depth 2 are at distance 0.5, while MRCAs at depths 40 and 41 give
0.025 and 0.024. It is also sensitive to the uneven resolution of the
classification: in The Taxonomicon, vertebrate lineages are typically several
times deeper than bacterial ones, so equal values do not correspond to
equivalent taxonomic levels in different groups.

The `metric` argument offers three alternatives. With \(d_A\) and \(d_B\) the
depths of the two lineages:

| Metric | Definition | Properties |
|---|---|---|
| `"reciprocal"` | \(1/h\) | ultrametric |
| `"path"` | \(d_A + d_B - 2h\) | tree metric (number of edges) |
| `"wu_palmer"` | \(1 - 2h/(d_A + d_B)\) | in \([0, 1)\); normalized by depth |
| `"ranked"` | half the path through standard ranks only | Clarke and Warwick step weights |

The `"ranked"` metric reduces each lineage to its root, the nodes of the
Linnaean ranks kingdom, phylum, class, order, family, genus, and species, and
the queried taxon. Unranked clades, which make up most of the nodes of
richly resolved groups, are ignored. For two species, the distance is the
number of rank levels up to their shared rank: 1 in the same genus, 2 in the
same family, and so on. Rank levels can be changed with the
`taxodist.rank_levels` option. Note that The Taxonomicon follows current
phylogenetic classifications, so that, for example, birds are the subclass
Aves of the class Reptilia.

```{r metrics}
lineages <- list(
  Panthera = c("Biota", "Animalia", "Chordata", "Mammalia", "Carnivora",
               "Felidae", "Panthera"),
  Felis    = c("Biota", "Animalia", "Chordata", "Mammalia", "Carnivora",
               "Felidae", "Felis"),
  Canis    = c("Biota", "Animalia", "Chordata", "Mammalia", "Carnivora",
               "Canidae", "Canis"),
  Homo     = c("Biota", "Animalia", "Chordata", "Mammalia", "Primates",
               "Hominidae", "Homo"),
  Musca    = c("Biota", "Animalia", "Arthropoda", "Insecta", "Diptera",
               "Muscidae", "Musca")
)
standard <- c(NA, "Kingdom", "Phylum", "Class", "Order", "Family", "Genus")
example <- taxo_from_lineages(
  lineages,
  ranks = rep(list(standard), length(lineages)) |> setNames(names(lineages))
)

for (metric in c("reciprocal", "path", "wu_palmer", "ranked")) {
  cat("\n", metric, "\n")
  print(round(distance_matrix(example, metric = metric, progress = FALSE), 3))
}
```

Rankings of pairs are the same under all four metrics within a single,
evenly resolved lineage set such as this one; they differ when lineages of
very different depths are compared.

## Taxonomic distinctness

`taxo_distinctness()` computes the indices of Warwick and Clarke: taxonomic
diversity \(\Delta\), taxonomic distinctness \(\Delta^*\), average taxonomic
distinctness \(\Delta^+\), and its variation \(\Lambda^+\). With the
`"ranked"` metric they reproduce the classical definitions based on Linnaean
ranks, and they agree with `vegan::taxondive()`; with the other metrics they
use the full resolution of the classification.

```{r distinctness}
d <- distance_matrix(example, metric = "ranked", progress = FALSE)
communities <- rbind(
  site_1 = c(Panthera = 4, Felis = 2, Canis = 1, Homo = 0, Musca = 0),
  site_2 = c(Panthera = 1, Felis = 0, Canis = 1, Homo = 2, Musca = 5)
)
taxo_distinctness(communities, d)
```

\(\Delta^+\) does not depend on sample size in expectation, so a sample can
be compared with random samples of the same number of taxa drawn from a
master list:

```{r distinctness-test}
set.seed(1)
taxo_distinctness_test(c("Panthera", "Felis", "Canis"), d, n_iter = 499)
```

## Dependence on the source classification

The numerical value depends on the number and arrangement of nodes returned by
The Taxonomicon. Taxonomic revisions, alternative classifications, and
differences in lineage resolution can change both MRCA depth and distance.

Comparisons should therefore use lineages obtained from the same source and,
where possible, the same package and data-retrieval version.

## Ambiguous and missing taxa

Some names correspond to more than one valid biological entry. A numeric
Taxonomicon identifier can be supplied when a specific taxonomic concept is
required.

Missing lineages produce missing pairwise values. Lineages without a common
root produce infinite values. These cases should be resolved before applying
methods that require a complete finite matrix, including clustering and
ordination.

## Interpretation

Smaller values indicate a deeper shared node in the source classification.
They should not be interpreted as:

- elapsed evolutionary time;
- genetic or morphological divergence;
- phylogenetic branch length;
- probability of common ancestry.

Hierarchical clustering of these distances produces a taxonomic similarity
dendrogram, not an independently inferred phylogenetic tree.

## Recommended reporting

Analyses using `taxodist` should report:

- the package version;
- the lineage retrieval date;
- The Taxonomicon as the source classification;
- the distance definition;
- decisions involving ambiguous or missing taxa;
- any clustering, ordination, or permutation method subsequently applied.

Formatted citations for the package and The Taxonomicon are available through
`citation("taxodist")`.