---
title: "Tree distance metrics"
author: "Martin R. Smith"
date: "`r Sys.Date()`"
output: 
  bookdown::pdf_document2:
    includes:
      in_header: ../inst/preamble.tex
  rmarkdown::html_document:
bibliography: ../inst/REFERENCES.bib
csl: https://raw.githubusercontent.com/citation-style-language/styles/master/apa-old-doi-prefix.csl

vignette: >
  %\VignetteIndexEntry{Tree distance metrics}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{R initialize, echo=FALSE}
library('Quartet')
cbPalette8 <- Ternary::cbPalette8

ColPlot <- function (tr, title=NULL, bold=NULL, ...) {
  tr$edge.length <- rep(1, dim(tr$edge)[1])
  font <- rep(1, length(tr$tip.label))
  if (!is.null(bold)) font[tr$tip.label %in% bold] <- 4
  plot(tr, tip.col=Ternary::cbPalette15[-c(4, 7)][as.integer(tr$tip.label)], 
       main=title, cex.main=0.8, font=font, y.lim=c(-3.5, 11), ...)
}

RFPlot <- function (tr, title=NULL, highlight=NULL, ref=sq_trees[[1]], ...) {
  tree_dist <- phangorn::treedist(tr, ref)
  tree_pair <- lapply(list(tr, ref), ape::root, outgroup='1', resolve.root=FALSE)
  class(tree_pair) <- 'multiPhylo'
  topo_dist <- as.matrix(ape::dist.topo(tree_pair, method='PH85'))[2]
  
  ColPlot(tr, title, highlight, cex=0.8, ...)
  text_x <- par('usr')[2] * 0.5
  
  text(text_x, -0.5, 'Quartet:', cex=0.8, pos=2)
  text(text_x, -1.5, "RF:", cex=0.8, pos=2)
  text(text_x, -2.5, "Path:", cex=0.8, pos=2)
  text(text_x, -3.5, "SPR:", cex=0.8, pos=2)
  text(text_x, -0.5, paste0(QuartetStatus(list(tr, ref))[2, 'd'],
                            '/', choose(11,4)), cex=0.8, pos=4)
  text(text_x, -1.5, paste0(topo_dist, '/', (11L - 3L) * 2L), cex=0.8, pos=4)
  text(text_x, -2.5, paste0(signif(tree_dist[2], 3)), cex=0.8, pos=4)
  text(text_x, -3.5, paste0(phangorn::sprdist(tr, ref)[1]), cex=0.8, pos=4)
}

Phangorn <- function (linkText) {
  paste0('[`', linkText, '`](https://www.rdocumentation.org/packages/phangorn/topics/treedist])')
}
```

# Tree distance metrics

A number of metrics area available to quantify the similarity between two undirected
topologies (i.e. unrooted trees with no edge lengths), many of which are
implemented in the `phangorn` function `r Phangorn('treedist')`.

## SPR metric
The subtree pruning and regrafting (SPR) distance [@Penny1985; 
`r Phangorn('SPR.dist')`] counts the 
number of SPR rearrangements necessary to transform Tree A into Tree B.

## Path difference metric
The length of a path from one tip to another in a tree is the number of edges
within the tree that must be crossed to navigate from one tip to the other.

Given two trees, is possible to calculate the difference in path length between
each pair of tips.

The path difference metric [@Steel1993; `r Phangorn('path.dist')`] is the 
square root of the sum of squares of each of these differences.

The related nodal distance [@Bluis2003] is not considered further here.

## Partition metric
The Robinson-Foulds (RF or 'partition') metric [@Robinson1981; @Steel1993;
phangorn's `r Phangorn('RF.dist')` or this package's `RobinsonFoulds`]
measures the symmetric difference between two trees by adding
the number of bipartitions that are present in tree A (but not tree B)
to the number of bipartitions present in tree B (but not tree A).  

It is most useful when the trees to be compared are very similar; it has
a low range of integer values, limiting its ability to distinguish between
trees [@Steel1993].

## Quartet metric
Instead of partitions, symmetric differences can be measured by counting the 
number of four-taxon statements (quartets) that differ between two trees
[@Estabrook1985;@Day1986].

For any four tips A, B, C and D, a bipartition on a bifurcating tree will separate
tip A and either B, C or D from the other two tips.  That is to say, removing
all other tips from the tree will leave one of these three trees:

```{R Three-four-taxon-trees, echo=FALSE, cache=TRUE, fig.width=6, fig.asp=1/3, out.width='66%', fig.align='center'}
par(mfrow=c(1, 3), mar=c(0.5, 1, 0.5, 1), cex=1)
plot(ape::read.tree(text='((A, B), (C, D));'), tip.color=cbPalette8[c(1, 4, 7, 5)], font=2)
plot(ape::read.tree(text='((A, C), (B, D));'), tip.color=cbPalette8[c(1, 7, 4, 5)], font=2)
plot(ape::read.tree(text='((A, D), (C, B));'), tip.color=cbPalette8[c(1, 5, 7, 4)], font=2)
```

Thus two of the random trees below share the quartet `(A, B), (C, D)`, whereas 
the third does not; these four tips are divided into `(A, D), (B, C)`.

```{R Plot-a-quartet, echo=FALSE, cache=TRUE, fig.asp=1.3/3, fig.width=6, out.width='80%', fig.align='center'}
par(mfrow=c(1, 3))
suppressWarnings(RNGversion("3.5.0")) # Stopgap until R 3.6.0 is widely available 
set.seed(7)
trees7 <- lapply(logical(3), function (X) {
    tr <- ape::rtree(7, br=NULL)
    tr$edge.length <- rep(1, 12)
    tr$tip.label <- LETTERS[1:7]
    tr
  })
PlotQuartet(trees7, LETTERS[1:4], cex=1.4, font=2)
```

There are $n\choose4$ groups of four taxa in a tree with $n$ tips;
for each of these groups, one of the three trees above will be consistent with
a given tree.  As such, two identical trees will have a quartet distance of
0, and a random pair of trees will have an expected ${n\choose{4}} / 3$
quartets in common. Because quartets are not independent of one another,
no pair of trees with six or more tips can have all $n\choose4$ quartets in
common [@Steel1993].

# Desired behaviour of tree distance metrics

The advantages of the quartet symmetric difference over other tree distance
metrics [@Penny1985] are best illustrated by examining a set of example trees.

## Moving a single taxon

If trees differ only in the location of a single taxon (see taxon 1 in the trees below),
then the distance between two trees should correspond to the distance that this taxon has been moved.

```{R Moving-a-single-taxon, fig.asp=1.8/5, out.width='100%', fig.width=6, echo=FALSE, message=FALSE, cache=TRUE}
par(mfrow=c(1, 4), mar=rep(0.4, 4))
ColPlot(sq_trees[[1]],     'Reference tree', 1)
RFPlot(sq_trees$move_one_near, 'Short move', 1)
RFPlot(sq_trees$move_one_mid, 'Medium move', 1)
RFPlot(sq_trees$move_one_far,   'Long move', 1)
```

The subtree pruning and regrafting (SPR) distance does not distinguish
between these trees, as they differ only in the placement of a single tip.
The Robinson-Foulds, path difference and quartet metrics, 
in contrast, recognize trees in which this tip has been 
moved further as more distant from the starting tree.

## Moving two taxa

Intuitively, moving a pair of tips on a tree should lead to higher tree distances
than moving a single tip.  In the case of a short move, the RF distance
does not differ whether one or two tips are moved.
For larger moves, however, the RF distance is _less_ when two tips are moved than
when a single tip is moved.  The path and quartet metrics perform as expected.

The trees below differ from a reference tree in the position of a single tip (tip 1), or a pair of tips (tips 10 and 11), which have been moved a short, medium or long distance from their original positions.

```{R Moving-a-cherry, echo=FALSE, message=FALSE, cache=TRUE, fig.asp=1.6*3/3, fig.width=6, out.width='66%', fig.align='center'}
par(mfrow=c(3, 3), mar=c(2.4, 0.4, 0.4, 0.4), cex=1)

ColPlot(sq_trees[[1]], 'Reference tree')
RFPlot(sq_trees$move_one_near, 'Short move 1', 1)
RFPlot(sq_trees$move_two_near, 'Short move 2', 10:11)

ColPlot(sq_trees[[1]], 'Reference tree')
RFPlot(sq_trees$move_one_mid, 'Medium move 1', 1)
RFPlot(sq_trees$move_two_mid, 'Medium move 2', 10:11)

ColPlot(sq_trees[[1]], 'Reference tree')
RFPlot(sq_trees$move_one_far,   'Long move 1', 1)
RFPlot(sq_trees$move_two_far,   'Long move 2', 10:11)

```

## Maximum distance
A distance metric should distinguish slightly-perturbed trees from
random trees and those that are more different from the starting tree than 
expected by chance.

The Robinson-Foulds metric can reach its maximum value when a single taxon is
relocated from the most basal to the most derived point of a pectinate tree, 
representing a maximal value despite retaining relationship information about 
all other taxa.

```{R Pectinate-to-break-RF, echo=FALSE, message=FALSE, cache=TRUE, fig.asp=1.5/3, fig.width=6, fig.align='center'}
par(mfrow=c(1, 3), mar=rep(0.4, 4))

pectinate_tree  <- ape::read.tree(text='(1, (2, (3, (4, (5, (6, (7, (8, (9, (10, 11))))))))));')
pectinate_unrooted <- ape::unroot(pectinate_tree)
pectinate_move1 <- ape::read.tree(text='(2, (3, (4, (5, (6, (7, (8, (9, (10, (11, 1))))))))));')
ColPlot(pectinate_tree, 'Pectinate tree')
RFPlot(pectinate_move1, 'Move one taxon', 1, ref=pectinate_unrooted)
RFPlot(sq_trees$random_tree, 'Random tree', 1, ref=pectinate_unrooted)
```

A notable proportion of random trees receive a lower RF distance from the original tree, even though they do not show any structural similarity.
This is not the case with the quartet symmetric difference.

```{R Random-trees, echo=FALSE, cache=TRUE, warn=FALSE, message=FALSE, fig.height=4.6, fig.width=5.6, out.width='80%', fig.align='center'}
pectinate_tree  <- ape::read.tree(text='(1, (2, (3, (4, (5, (6, (7, (8, (9, (10, 11))))))))));')
pectinate_unrooted <- ape::unroot(pectinate_tree)
pectinate_move1 <- ape::read.tree(text='(2, (3, (4, (5, (6, (7, (8, (9, (10, (11, 1))))))))));')

n_random <- if (Sys.getenv('TRAVIS') == "") 10000 else 10
random_trees <- lapply(seq_len(n_random), function (x) ape::rtree(11, br=NULL, tip.label=1:11))
random_quartets <- QuartetStatus(random_trees, cf=pectinate_tree)[, 'd'] / choose(11, 4)
random_distances <- vapply(random_trees, phangorn::treedist,
                           tree2=pectinate_tree, double(2))
random_spr <- vapply(random_trees, phangorn::sprdist, tree2=pectinate_tree, double(4))

par(mfrow=c(2, 2), mar=rep(1.9, 4), cex=0.8, oma=c(0, 0, 2, 0))
manyBreaks <- seq(0, 1, length.out=40)

PlotHist <- function (distances, breaks, pectDist, title) {
  hist(distances, xlab=NULL, main=NULL, breaks=breaks, cex=0.8)
  legend('topleft', bty='n', c('', '', title), cex=0.8)
  abline(v=pectDist, col=cbPalette8[3], lwd=2)
}

PlotHist(random_distances[1, ] / 16, manyBreaks, 1, 'Robinson-Foulds distance')
mtext('Distance between pectinate tree (above) and random trees', outer=TRUE, cex=0.8)
legend('topleft', bty='n', col=cbPalette8[3], lwd=2, '"Move one taxon" tree', cex=0.8)

PlotHist(random_distances[2, ], breaks=manyBreaks*max(random_distances[2, ]),
         17.7, "Path distance")

PlotHist(random_quartets, manyBreaks, 120 / 330, 'Quartet dissimilarity')

PlotHist(random_spr[1, ], max(random_spr[1, ]) * manyBreaks, 1, 'SPR distance')
```

An advantage of the quartet symmetric distance is that the normalized metric of a random
tree is $\frac{2}{3}$ [@Day1986; @Steel1993].  As such, trees that are more different
than expected by chance can be readily recognized, as their distance metric
will be greater than $\frac{`r choose(11, 4) * 2 / 3`}{`r choose(11, 4)`}$.
The 'maximum distance' tree depicted below was identified using the R package
`TreeSearch` [@TreeSearch], using the quartet difference from the reference tree 
as an optimality criterion.

```{R Increasing-distances, echo=FALSE, message=FALSE, cache=TRUE, fig.asp=1.4/3, fig.width=6, out.width='90%', fig.align='center'}
par(mfrow=c(1, 3), mar=rep(0.4, 4))
ColPlot(sq_trees[[1]], 'Reference tree')
RFPlot(sq_trees$random_tree,   'Random tree'     )
RFPlot(sq_trees$opposite_tree, 'Maximum distance')
```

## Unit equivalence

A further shortcoming of the RF metric is that not all partitions represent an equivalent amount of information.  A partition distance of 1 could mean that two trees differ in an uninformative partition, or a more informative partition.  All quartets,
in contrast, are equally informative.

Consider a balanced and an unbalanced eight-taxon tree:

```{r two-trees, echo=FALSE, fig.height = 3, fig.width=6}
par(mfrow=c(1, 2), mar=rep(0.8, 4), cex=0.9)
plot(balancedTree <- ape::read.tree(text="(((a, b), (c, d)), ((e, f), (g, h)));"),
     edge.color = cbPalette8[3])
ape::edgelabels(1:5, c(1, 2, 5, 9, 12))
legend('topleft', 'Balanced tree', bty='n')

plot(caterpillarTree <- ape::read.tree(text="(a, (b, (c, (d, (e, (f, (g, h)))))));"),
     edge.color = cbPalette8[2])
ape::edgelabels(1:5, c(4, 6, 8, 10, 12))
legend('topleft', 'Asymmetric tree', bty='n')
```

Each tree divides the eight taxa into five bipartition splits.  
The information content (Shannon entropy) of a split can be calculated based on 
what proportion of eight-tip trees contain the split in question. 
This is a function of the evenness of the split:

```{r ic-of-splits, display='asis', echo=FALSE}
splitSmall <- 2:4
splitLarge <- 8L - splitSmall

rootedTrees <- c(
  '1' = 1,
  '2' = 1 * 1,
  '3' = 3 * 1,
  '4' = 5 * 3 * 1,
  '5' = 7 * 5 * 3 * 1,
  '6' = 9 * 7 * 5 * 3 * 1
)

matchingTrees <- rootedTrees[splitSmall] * rootedTrees[splitLarge]
names(matchingTrees) <- paste0("Partition size: ", splitSmall, ':', splitLarge)

matchingP <- matchingTrees / (11 * 9 * 7 * 5 * 3 * 1 * 1)

ic <- -log(matchingP) / log(2)

knitr::kable(cbind(
      'Matching trees' = matchingTrees,
      'p(Match in random tree)' = signif(matchingP, 3),
      'Information content / bits' = ic))
```

In the first tree, split 1 is even, dividing four taxa from four others (`4:4`); splits 2--5 are maximally uneven (`2:6`).  The total information content of these five splits is `r signif(sum(ic[c(4, 2, 2, 2, 2) - 1]), 4)`, whereas that of the five splits in 
the second tree, of sizes `2:6`, `3:5`, `4:4`, `3:5` and `2:6`, is
`r signif(sum(ic[c(2, 3, 4, 3, 2) - 1]), 4)`. Put another way, a random tree will
on average share more partitions with the balanced tree (whose partitions are 
predominantly uneven and thus likely to be matched) than the asymmetric tree 
(which contains more even partitions that are less likely to occur in a random tree).

Of the 10&nbsp;395 eight-tip trees, many more bear at least one partition in 
common with a balanced tree than with an asymmetric tree, whereas the distribution of quartets is identical:

```{r all-8-tip-trees, echo=FALSE, cache=TRUE}
#calculate <- FALSE
#if (calculate) {
#  all8 <- phytools::allFurcTrees(8, letters[1:8], FALSE)
#  allBif8 <- all8[vapply(all8, function(x) x$Nnode, 1L) == 6]
#  inBalanced <- SplitStatus(allBif8, balancedTree)[, 's']
#  inCaterpillar <- SplitStatus(allBif8, caterpillarTree)[, 's']
#  qInBalanced <- QuartetStatus(allBif8, balancedTree)[, 's']
#  qInCaterpillar <- QuartetStatus(allBif8, balancedTree)[, 's']
#} else {
  inBalanced <- rep(0:5, c(7088, 2708, 512, 76, 10, 1))
  inCaterpillar <- rep(0:5, c(8162, 1808, 350, 64, 10, 1))
  qInBalanced <- qInCaterpillar <- rep(c(10:42, 44:46, 49, 50, 53, 54, 57, 62, 70),
    c(32, 256, 512, 384, 512, 320, 512, 448, 512, 512, 688, 160, 512, 448, 256, 464, 416, 448, 576, 256, 192, 320, 160, 288, 128, 192, 132, 80, 32, 32, 224, 64, 8, 64, 112, 16, 64, 16, 4, 18, 16, 8, 1))
#}

par(mfrow=c(1, 2))
plot.new()
plot.window(main='8-tip trees with N common splits', cex.main=1,
            xlim=c(0, 6), ylim=c(0, 2708),
            xlab='Splits in common', ylab='Number of trees')
sch <- hist(inCaterpillar + 0.5, breaks = 0:12 / 2 - 0.0001, plot=FALSE)
sbh <- hist(inBalanced, breaks = 0:12 / 2 - 0.0001, plot=FALSE)
plot(sch, col=paste0(cbPalette8[2], '44'), add=TRUE)
plot(sbh, col=paste0(cbPalette8[3], '44'), add=TRUE)
text(0.0, 100, paste0('Balanced: ', sum(inBalanced == 0)), pos=4, srt=90, cex=0.7)
text(0.5, 100, paste0('Asymmetric: ', sum(inCaterpillar == 0)), pos=4, srt=90, cex=0.7)

legend('topright', pch=22,
       pt.cex = 2, col='black',
       pt.bg = paste0(cbPalette8[2:3], '44'), bty='n',
       c('Asymmetric', 'Balanced'))
axis(1, at=0:5 + 0.5, labels=0:5)
axis(2)

qch <- hist(qInCaterpillar, breaks = 0:36 * 2 - 0.0001, plot=FALSE)
qbh <- hist(qInBalanced,    breaks = 0:36 * 2 - 0.0001, plot=FALSE)

plot(qch, col=paste0(cbPalette8[2], '44'), main='8-tip trees with N common quartets',
     xlab='Quartets in common', ylab='Number of trees', axes=FALSE, cex.main=1)
plot(qbh, col=paste0(cbPalette8[3], '44'), add=TRUE)
axis(1, at=0:36 * 2, labels=(0:36) * 2)
axis(2)

```

## Unresolved trees

Whereas the path distance and SPR metrics are only defined on bifurcating trees,
symmetric difference approaches can be applied to trees that contain polytomies -- i.e. 
not every node is resolved as bifurcating.

## Conclusion

Quartet dissimilarity is the only available metric of tree distance that fulfils 
all of the following desiderata:

 - Allocates trees higher distances if a clade moves greater distances
 - Allocates trees higher distances if a the clade that is moved is larger
 - Distinguishes contradicted from unresolved information in trees that 
   are not fully bifurcating (resolved)
 - Identifies pairs of trees that are more random than expected by chance
 - Does not reach its maximum value after relatively trivial rearrangements
 
# Improvements to partition metrics

The primary issue with the Robinson-Foulds distance is its conservative 
pairing of clades: a grouping is either identical in both trees, or is 
treated as different.  
This makes no account for the possibility that one tree contains a grouping that
is similar (though not quite identical) to a grouping in another.

In order to incorporate some measure of similar-but-not-quite-identical clades,
it is necessary to construct an optimal matching that pairs each clade in each
tree with either a suitable clade in the other, or with no clade at all.

Various approaches to this problem have been proposed; these typically search
for a matching that minimises the total dissimilarity summed for all pairings.
Approaches differ in how they penalize differences between paired groups.
Bogdanowicz and Giaro [-@Bogdanowicz2012;-@Bogdanowicz2013;-@Bogdanowicz2017]
calculate the size
of the largest bipartition that is consistent with both input groups; 
Nye _et al_. [-@Nye2006] normalize this number using the Jaccard index;
Böcker _et al_. [-@Bocker2013] raise this value to an arbitrary exponent, and
add a constraint that prohibits pairings that comprise contradictory groups.

These approaches represent a refinement of the Robinson-Foulds metric,
circumventing many of the issues raised above.  They are, however, more 
computationally expensive, and lack an objective measure of clade similarity.
I am not presently aware of an R implementation of any of these measures,
though a related metric is planned for a future version of this package.

# References
 