---
title: "Supplementary text for: Bayesian and parsimony approaches reconstruct informative trees from simulated morphological datasets"
author: "Martin R. Smith (martin.smith@durham.ac.uk)"
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/dependent/biology-letters.csl
link-citations: yes
github-repo: ms609/Quartet
vignette: >
%\VignetteIndexEntry{Tree comparison using quartets}
%\VignetteEncoding{UTF-8}
%\VignetteEngine{knitr::rmarkdown}
---
---
This document has been generated from an R markdown file, which contains the source
code used to generate figures and conduct analyses, and is provided in the Electronic Supplementary Material that accompanies the main article [@Smith2019].
---
```{r setup, include=FALSE}
require('ape')
require('Ternary')
require('Quartet')
if (!require('CongreveLamsdell2016')) {
devtools::install_github('ms609/CongreveLamsdell2016')
library('CongreveLamsdell2016')
}
Sys.setlocale("LC_ALL", "English_United Kingdom.1252") # allows PCH=183
data('sq_trees')
n_tip <- 11
ref_tree <- sq_trees[[1]]
tip_colours <- Ternary::cbPalette15[-c(4, 7)] # Rm Tritanopia duplicates of 13 and 3
palette4 <- Ternary::cbPalette8[1:4]
PCH_MK <- 1 # circle
PCH_EQ <- 61 #'='
PCH_XX <- 183 #'.'
PCH_IW <- 3
PCH_IC <- 17 #triupfilled
COL_MK <- paste0(cbPalette8[4], '99')
COL_EQ <- paste0(cbPalette8[8], '99')
COL_1 <- paste0(cbPalette8[6], '42')
COL_2 <- paste0(cbPalette8[6], '42')
COL_3 <- paste0(cbPalette8[6], '42')
COL_5 <- paste0(cbPalette8[6], '99')
COL10 <- paste0(cbPalette8[6], '42')
COL_C <- paste0(cbPalette8[2], '99')
PCH <- c(
brem = 2,
bootG = 0,
bootF = 5,
jackG = 3,
jackF = 4,
none = NA,
mk = 1,
markov = 1,
equal = 61, #'='
eq = 61,
dot = 183, #'.'
k1 = 4,
implied1 = 4,
k2 = 183,
implied2 = 183,
k3 = 183,
implied3 = 183,
k5 = 3,
implied5 = 3,
kX = 183,
implied10 = 183,
implied20 = 183,
implied200 = 183,
kC = 17,
impliedC = 17 #triupfilled
)
COL <- c(
black = paste0(cbPalette8[1], '99'),
markov = paste0(cbPalette8[4], '99'),
mk = paste0(cbPalette8[4], '99'),
equal = paste0(cbPalette8[8], '99'),
eq = paste0(cbPalette8[8], '99'),
implied1 = paste0(cbPalette8[6], '42'),
k1 = paste0(cbPalette8[6], '42'),
implied2 = paste0(cbPalette8[6], '42'),
k2 = paste0(cbPalette8[6], '42'),
implied3 = paste0(cbPalette8[6], '42'),
k3 = paste0(cbPalette8[6], '42'),
implied5 = paste0(cbPalette8[6], '42'),
k5 = paste0(cbPalette8[6], '42'),
implied10 = paste0(cbPalette8[6], '99'),
kX = paste0(cbPalette8[6], '99'),
implied20 = paste0(cbPalette8[6], '42'),
implied200 = paste0(cbPalette8[6], '42'),
kC = paste0(cbPalette8[2], '99'),
impliedC = paste0(cbPalette8[2], '99')
)
ANALYSES_EN <- c('Equal weights', paste0('IW, k=', c(1, 2, 3, 5, 10)), 'IW consensus')
GRID_COL <- rgb(0.92, 0.92, 0.92)
BG_COL <- rgb(0.985, 0.985, 0.992)
rf_max <- (n_tip - 3) * 2
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=tip_colours[as.integer(tr$tip.label)], main=title, cex.main=0.8, font=font, y.lim=c(-3.5, n_tip), ...)
}
rfplot <- function (tr, title=NULL, highlight=NULL, ref=ref_tree, ...) {
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, '/', rf_max), 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)
}
polyplot <- function (tr, title, highlight, ...) {
tree_pair <- lapply(list(ref_tree, tr), ape::unroot)
class(tree_pair) <- 'multiPhylo'
partitions <- SplitStatus(tree_pair)[2, ]
colplot(tr, title, highlight, ...)
x_mid <- par('usr')[2] * 0.62
text(x_mid, -0.5, cex=0.8, pos=2, "Quartets contradicted")
text(x_mid, -1.5, cex=0.8, pos=2, "Quartets unresolved")
text(x_mid, -2.5, cex=0.8, pos=2, "Partitions contradicted")
text(x_mid, -3.5, cex=0.8, pos=2, "Partitions unresolved")
text(x_mid, -0.5, cex=0.8, pos=4,
paste0(QuartetStatus(tree_pair)[2, 'd'], '/', choose(11, 4)))
text(x_mid, -1.5, cex=0.8, pos=4,
paste0(sum(QuartetStatus(tree_pair)[2, c('r1', 'r2', 'u')]), '/', choose(11, 4)))
# The below assumes that ref is bifurcated, so if a partition's in cf not ref it must contradict.
text(x_mid, -2.5, cex=0.8, pos=4, paste0(partitions['cf_not_ref'], '/', rf_max / 2))
text(x_mid, -3.5, cex=0.8, pos=4, paste0(partitions['ref_not_cf'] - partitions['cf_not_ref'], '/', rf_max / 2))
}
data('clPartitions', 'clQuartets')
Quartet2Ternary <- function (item) clQuartets[[item]][, c('r2', 'd', 's'), TREE]
Split2Ternary <- function (item) {
itemData <- clPartitions[[item]][, , TREE]
rbind(itemData[, 'ref'] - itemData[, 'cf'],
itemData[, 'cf_not_ref'],
itemData[, 'cf_and_ref'])
}
TernaryBasis <- function (qMax = 19, zoom = 1, padding=0.1) {
xLim <- c(0, 1 / zoom) - 0.05
yLim <- c(0.5 - (1 / zoom), 0.5)
lab <- if (zoom == 1) c('Unresolved', "Different", 'Same') else rep('', 3)
TernaryPlot(NULL, NULL, NULL,
lab[1], lab[2], lab[3], lab.cex=0.7,
point='right',
col=BG_COL,
grid.lty='solid', grid.col=GRID_COL, grid.lines=19,
grid.minor.lines=0,
axis.labels = round(seq(0, qMax, length.out=20), 0),
ticks.length = if (zoom == 1) 0.025 else 0.006,
axis.col=rgb(0.6, 0.6, 0.6),
padding=padding, xlim=xLim, ylim=yLim)
HorizontalGrid(19) # must equal grid.lines
}
AddArrows <- function (quality) {
arrows(sqrt(3/4) * 0.5, 0.5, sqrt(3/4) * 0.8, 0.5, length=0.08)
text (sqrt(3/4) * 0.65, 0.5, pos=3, 'Decreasing resolution', cex=0.8)
arrows(sqrt(3/4) * 0.98, 0.40, sqrt(3/4) * 0.98, 0.20, length=0.08)
text (sqrt(3/4) * 1.01, 0.30, pos=3, quality, cex=0.8, srt=270)
}
AddLegend <- function(pos='center')
legend(pos, cex=0.8, bty='n',
lty=1,
pch=c(PCH_MK, PCH_EQ, PCH_XX, PCH_IW, PCH_XX, PCH_XX, PCH_IW, PCH_IC), pt.cex=1.1,
col=c(COL_MK, COL_EQ, COL10, COL_5, COL_3, COL_2, COL_1, COL_C),
legend=c('Markov', 'Equal weights', paste0('Implied, k=', c(10, 5, 3, 2, 1, '2..10')))
)
AverageSplits <- function (item) {
itemData <- apply(clPartitions[[item]][, , ], 2, rowMeans)
rbind(itemData['cf_and_ref', ],
itemData['cf_not_ref', ],
itemData['ref', ] - itemData['cf', ])
}
AverageQuarts <- function (item) apply(clQuartets[[item]][c('r2', 'd', 's'), , ], 2, rowMeans)
```
# 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).
## SPR metric
The subtree pruning and regrafting (SPR) distance [@Penny1985] 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] is the square root of the sum of squares
of each of these differences.
## Partition metric
The Robinson-Foulds (RF or 'partition') metric [@Robinson1981; @Steel1993]
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=Ternary::cbPalette8[c(1, 4, 7, 5)], font=2)
plot(ape::read.tree(text='((A, C), (B, D));'), tip.color=Ternary::cbPalette8[c(1, 7, 4, 5)], font=2)
plot(ape::read.tree(text='((A, D), (C, B));'), tip.color=Ternary::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))
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(ref_tree, '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(ref_tree, 'Reference tree')
rfplot(sq_trees$move_one_near, 'Short move 1', 1)
rfplot(sq_trees$move_two_near, 'Short move 2', 10:11)
colplot(ref_tree, 'Reference tree')
rfplot(sq_trees$move_one_mid, 'Medium move 1', 1)
rfplot(sq_trees$move_two_mid, 'Medium move 2', 10:11)
colplot(ref_tree, '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=Ternary::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=Ternary::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(n_tip, 4) * 2 / 3`}{`r choose(n_tip, 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(ref_tree, '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 <- read.tree(text="(((a, b), (c, d)), ((e, f), (g, h)));"),
edge.color = cbPalette8[3])
edgelabels(1:5, c(1, 2, 5, 9, 12))
legend('topleft', 'Balanced tree', bty='n')
plot(caterpillarTree <- read.tree(text="(a, (b, (c, (d, (e, (f, (g, h)))))));"),
edge.color = cbPalette8[2])
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 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
# Calculating resolution and accuracy
One way to modify a tree topology is to reduce its resolution by collapsing nodes, without changing any of the relationships presented within the tree. The trees below have been derived from a reference tree by collapsing one and many nodes:
```{R Polytomies-same-topology, echo=FALSE, message=FALSE, cache=TRUE, fig.asp=1.6/3, fig.width=6, out.width='90%', fig.align='center'}
par(mfrow=c(1, 3), mar=c(0.4, 0.4, 0.4, 0.4))
colplot(ref_tree, 'Reference tree')
polyplot(sq_trees$collapse_one, 'One polytomy', 7:9)
polyplot(sq_trees$collapse_some, 'Much polytomy', 1:9)
```
These trees do not contain any quartets or partitions that are not present in the reference tree, though they do contain a number of unresolved quartets and partitions.
We can alternatively choose to change the topology, and then collapse some nodes. The following trees represent the same loss of resolution as the previous trees, but applied to a tree in which one tip (tip 1) has been moved relative to the reference tree:
```{R Polytomies 1, echo=FALSE, message=FALSE, cache=TRUE, fig.asp=1.6/3, fig.width=6, out.width='90%', fig.align='center'}
par(mfrow=c(1, 3), mar=c(0.4, 0.4, 0.4, 0.4))
polyplot(sq_trees$move_one_mid, 'One moved, no polytomy', 1)
polyplot(sq_trees$m1mid_col1, 'One moved, one polytomy', 1)
polyplot(sq_trees$m1mid_colsome, 'One moved, much polytomy', 1)
```
This causes results in trees that contradict a number of partitions or quartets that occurred in the original reference tree.
We could introduce a larger change to the tree topology by moving a 'cherry', i.e. two adjacent taxa (10 & 11):
```{R Polytomies-2, echo=FALSE, message=FALSE, cache=TRUE, fig.asp=1.6/3, fig.width=6, out.width='90%', fig.align='center'}
par(mfrow=c(1, 3), mar=c(0.4, 0.4, 0.4, 0.4))
polyplot(sq_trees$move_two_mid, 'Two moved, no polytomy' , 10:11)
polyplot(sq_trees$m2mid_col1, 'Two moved, one polytomy' , 10:11)
polyplot(sq_trees$m2mid_colsome, 'Two moved, much polytomy', 10:11)
```
## Visualizing these data
The number of quartets or partitions that are unresolved, different, or identical to the reference tree can be visualized using ternary diagrams:
```{R Plot-some-trees, echo=FALSE, message=FALSE, cache=TRUE, fig.asp=1/2, fig.width=8, out.width='100%', fig.align='center'}
par(mfrow=c(1, 2), mar=rep(0.3, 4))
collapse_perfect <- c('ref_tree', 'collapse_one', 'collapse_some')
collapse_move1 <- c('move_one_mid', 'm1mid_col1', 'm1mid_colsome')
collapse_move2 <- c('move_two_mid', 'm2mid_col1', 'm2mid_colsome')
TernaryBasis(choose(22, 4))
title(main="\nQuartet distances", cex.main=0.8)
AddToTernary(lines, list(c(0, 2/3, 1/3), c(1, 0, 0)), lty='dotted', col=cbPalette8[8], lwd=2)
pchs <- c(1, 4, 3)
quartet_distances <- QuartetPoints(sq_trees)
TernaryPoints(quartet_distances[collapse_perfect, ], col=cbPalette8[2], pch=pchs, cex=1.5)
TernaryLines (quartet_distances[collapse_perfect, ], col=cbPalette8[2])
TernaryPoints(quartet_distances[collapse_move1, ], col=cbPalette8[3], pch=pchs, cex=1.5)
TernaryLines (quartet_distances[collapse_move1, ], col=cbPalette8[3])
TernaryPoints(quartet_distances[collapse_move2, ], col=cbPalette8[4], pch=pchs, cex=1.5)
TernaryLines (quartet_distances[collapse_move2, ], col=cbPalette8[4])
legend('bottomright',
lty=rep('dotted', 2),
lwd=c(1, 2),
col=c('grey', cbPalette8[8]),
legend=c("Lines of equal RF distance\n/ quartet divergence", 'Mean distance of random tree'),
cex=0.8, bty='n')
arrows(sqrt(3/4) * 0.5, 0.5, sqrt(3/4) * 0.8, 0.5, length=0.1)
text (sqrt(3/4) * 0.65, 0.5, pos=3, 'Decreasing resolution', cex=0.8)
arrows(sqrt(3/4) * 0.98, 0.40, sqrt(3/4) * 0.98, 0.20, length=0.1)
text (sqrt(3/4) * 1.01, 0.30, pos=3, 'Increasing normalized SD', cex=0.8, srt=270)
TernaryBasis(19)
title(main="\nPartition distances", cex.main=0.8)
partition_distances <- SplitPoints(sq_trees)
TernaryPoints(partition_distances[collapse_perfect, ], col=cbPalette8[2], pch=pchs, cex=1.5)
TernaryLines (partition_distances[collapse_perfect, ], col=cbPalette8[2])
TernaryPoints(partition_distances[collapse_move1, ], col=cbPalette8[3], pch=pchs, cex=1.5)
TernaryLines (partition_distances[collapse_move1, ], col=cbPalette8[3])
TernaryPoints(partition_distances[collapse_move2, ], col=cbPalette8[4], pch=pchs, cex=1.5)
TernaryLines (partition_distances[collapse_move2, ], col=cbPalette8[4])
arrows(sqrt(3/4) * 0.98, 0.40, sqrt(3/4) * 0.98, 0.20, length=0.1)
text (sqrt(3/4) * 1.01, 0.30, pos=3, 'Increasing normalized SD', cex=0.8, srt=270)
legend('bottom', cex=0.8, bty='n',
pch=20, pt.cex=1.5, col=cbPalette8[2:4],
title='Move:',
legend=c('No tips', 'One tip', 'Two tips'))
legend('bottomright', cex=0.8, bty='n',
pch=pchs, pt.cex=1.5, col='black',
title='Collapse:',
legend=c('No nodes', 'One node', 'Many nodes')
)
```
In these plots, the vertical direction corresponds to the normalized symmetric distance. Collapsing nodes decreases the resolution (movement in the horizontal direction), but can increase accuracy; the balance between resolution lost and accuracy gained determined whether the collapsing of nodes increases or decreases net divergence.
# Using ternary diagrams to inform tree reconstruction techniques
This means of visualization provides a helpful way to understand how effective
different methods of phylogenetic reconstruction are on particular trees.
Here I have taken a representative dataset simulated from a 22-tip reference tree [@Congreve2016],
and analysed the dataset in TNT v1.5 [@Goloboff2016] under equal weights parsimony
and implied weights (with concavity constants of 1, 2, 3, 5 and 10), and in MrBayes v3.2.2
[@Huelsenbeck2001] using the Markov K model [@Lewis2001].
For each parsimony analysis, I recorded a strict consensus of all optimal trees,
then proceeded to collapse groups with a Bootstrap GC support under -95, -90, -85...
For each Bayesian analysis, I collapsed all groups whose posterior probability was
under $x$, at 20 uniformly spaced values of $x$ from 0.5 to 1.0.
These trees represent a progressive loss of resolution (precision) from the
optimal tree, allowing an exploration of the relationship between resolution and
accuracy. In each plot, resolution decreases from left to right.
Congreve and Lamsdell [-@Congreve2016] argue that equal weighting is the optimal
superior method because it resolves the fewest incorrect bipartitions -- that
is, its most-resolved tree is the closest to the top-right edge of the ternary diagram.
By this measure, all methods are improved by collapsing nodes until none remain.
On the view advocated here, the optimal tree is the one that has the lowest
normalized distance from the generative tree, which corresponds to the greatest
position in the vertical direction. (The normalizing constant is the maximum
possible number of partitions or quartets that could be resolved, not the number
that are actually resolved in a pair of trees.)
By this measure, collapsing the least-supported nodes leads to an increase in
tree quality, as predicted by Goloboff [-@Goloboff1995]: nodes with low support
are likely to be incorrect.
Collapsing better-supported nodes, however, reduces tree quality:
nodes with high support are likely to be correct.
```{R Congreve-Lamsdell-tree, echo=FALSE, message=FALSE, cache=TRUE, fig.asp=1/3, fig.width=9, out.width='100%', fig.align='center'}
par(mfrow=c(1, 3), mar=rep(0.3, 4))
TREE <- 2
clInitializeTernaryQuarts()
clPlotQuartets(clBootGcQuartets, TREE, cex=1.1, pch=PCH)
clPlotQuartets(clMkvQuartets, TREE, cex=1.1, pch=PCH['mk'])
title(main="\nQuartets", cex.main=0.8)
arrows(sqrt(3/4) * 0.5, 0.5, sqrt(3/4) * 0.8, 0.5, length=0.1)
text (sqrt(3/4) * 0.65, 0.5, pos=3, 'Decreasing resolution', cex=0.8)
arrows(sqrt(3/4) * 0.98, 0.40, sqrt(3/4) * 0.98, 0.20, length=0.1)
text (sqrt(3/4) * 1.01, 0.30, pos=3, 'Increasing divergence', cex=0.8, srt=270)
clInitializeTernaryQuarts(zoom=6, padding=0.01)
clPlotQuartets(clBootGcQuartets, TREE, cex=1.1, pch=PCH)
clPlotQuartets(clMkvQuartets, TREE, cex=1.1, pch=PCH['mk'])
clInitializeTernarySplits()
title(main="\nPartitions", cex.main=0.8)
clPlotSplits(clBootGcPartitions, TREE, cex=1.1, pch=PCH)
clPlotSplits(clMkvPartitions, TREE, cex=1.1, pch=PCH['mk'])
arrows(sqrt(3/4) * 0.98, 0.40, sqrt(3/4) * 0.98, 0.20, length=0.1)
text (sqrt(3/4) * 1.01, 0.30, pos=3, 'Increasing RF distance', cex=0.8, srt=270)
AddLegend <- function(pos='bottomright')
legend(pos, cex=0.8, bty='n',
lty=1,
pch=PCH[c('mk', 'eq', 'kX', 'k5', 'k3', 'k2', 'k1', 'kC')], pt.cex=1.1,
col=COL[c('mk', 'eq', 'kX', 'k5', 'k3', 'k2', 'k1', 'kC')],
legend=c('Markov', 'Equal weights', paste0('Implied, k=', c(10, 5, 3, 2, 1, '2..10')))
)
AddLegend()
```
```{r cl-best-perfect, echo=FALSE, results="asis"}
sliceSize <- dim(clBootGcQuartets[[1]])[-2]
treeQuality <- vapply(clBootGcQuartets,
function (item) apply(item, 3, function (slice)
apply(slice, 1, QuartetDivergence, similarity = FALSE)),
matrix(double(prod(sliceSize)), sliceSize[1]))
treeResolved <- vapply(clBootGcQuartets,
function (item) apply(item, 3, function (slice)
slice[, 'r2'] == 0),
matrix(logical(prod(sliceSize)), sliceSize[1]))
bestTrees <- apply(treeQuality, 2, which.min)
flatResolution <- apply(treeResolved, 2, unlist)
perfects <- sum(vapply(seq_len(dim(flatResolution)[2]),
function (i) flatResolution[bestTrees[i], i],
logical(1)))
```
With this particular dataset, the optimal tree is not perfectly resolved in any method. In fact, of all 100 datasets, the best available tree was only perfectly resolved
in a `r perfects`% of cases. The best tree available under a specific method
was perfectly bifurcating in 0--22% of cases:
```{r cl-perfect, echo=FALSE, results="asis"}
perfectlyResolved <- vapply(clBootGcQuartets, function (item)
sum(apply(item, 3, function (slice) slice[which.min(apply(slice, 1, QuartetDivergence, similarity = FALSE)), 'r2']) == 0), double(1))
names(perfectlyResolved) <- ANALYSES_EN
knitr::kable(t(perfectlyResolved))
```
With the 100, 350 and 1000 character datasets of O'Reilly _et al._ [-@OReilly2016],
the best available tree was almost never perfectly resolved.
```{r or-perfect, echo=FALSE, results="asis"}
library('OReillyEtAl2016')
data('orQuartets')
perfects <- vapply(orQuartets, function (dataset) {
sliceSize <- dim(dataset[[1]])[-2]
treeQuality <- vapply(dataset,
function (item) apply(item, 3, function (slice)
apply(slice, 1, QuartetDivergence, similarity = FALSE)),
matrix(double(prod(sliceSize)), sliceSize[1]))
treeResolved <- vapply(dataset,
function (item) apply(item, 3, function (slice)
slice[, 'r2'] == 0),
matrix(logical(prod(sliceSize)), sliceSize[1]))
bestTrees <- apply(treeQuality, 2, which.min)
flatResolution <- apply(treeResolved, 2, unlist)
sum(vapply(seq_len(dim(flatResolution)[2]),
function (i) flatResolution[bestTrees[i], i],
logical(1))) / sliceSize[2]
}, double(1))
perfectlyResolved <- rbind(perfects, vapply(orQuartets, function (dataset)
vapply(dataset,
function (item) sum(apply(item, 3, function (slice)
slice[which.min(apply(slice, 1, QuartetDivergence,
similarity = FALSE)), 'r2']) == 0) * 100 / dim(item)[3],
double(1)),
double(length(orQuartets[[1]]))))
rownames(perfectlyResolved) <- c("Best", "Markov", ANALYSES_EN[c(1, 3)],
"k=3", "k=5", "k=10",
"k=20", "k=200", ANALYSES_EN[7])
knitr::kable(t(round(perfectlyResolved, 2)))
```
## Quartet metric
We can also examine the situation if we average across all 100 Congreve &
Lamsdell datasets:
```{R CL-trees-quartets, echo=FALSE, message=FALSE, cache=TRUE, fig.asp=1/2, fig.width=9, out.width='100%', fig.align='center'}
TREES <- 1:100
# # There is essentially no difference if we exclude datasets with a low CI.
# # The lines just slide up a little.
# data(clCI)
# TREES <- TREES[clCI > 0.26]
par(mfrow=c(1, 2), mar=rep(0.3, 4))
clInitializeTernaryQuarts()
clPlotAverageQuartets(clBootGcQuartets, cex=1.1, pch=PCH)
clPlotTheseAverageQuartets(clMkvQuartets, cex=1.1, pch=PCH['mk'], col=COL['mk'])
clInitializeTernaryQuarts(zoom=3.5, padding=0.01)
clPlotAverageQuartets(clBootGcQuartets, cex=1.1, pch=PCH)
clPlotTheseAverageQuartets(clMkvQuartets, cex=1.1, pch=PCH['mk'], col=COL['mk'])
AddLegend('topright')
```
As the worst-supported nodes are progressively collapsed, the accuracy of implied weights trees begins to increase, decreasing the total distance of trees from the generative
tree (thus meaning that the provide a better representation of the generative tree).
One question we can ask is how much we should reduce the resolution. After a certain point, the increase in accuracy gained by collapsing the least supported nodes no longer offsets the information lost by sacrificing resolution.
Bayesian and equal weights trees already produce incompletely resolved trees, and a further reduction of resolution does not improve their quality.
Under implied weights, averaged over these datasets, the optimum trade-off between accuracy and resolution comes when collapsing nodes with a Bootstrap GC support value below:
```{R Quartet-divergence-for-each-point, echo=FALSE, message=FALSE, fig.width=8, results='asis', fig.align='center'}
divergences <- vapply(clBootGcQuartets, function (item) QuartetDivergence(apply(item[, , TREES], 2, rowMeans), similarity = FALSE), double(dim(clBootGcQuartets[[1]])[1]))
rownames(divergences) <- seq(-100, 100, length.out=dim(clBootGcQuartets[[1]])[1])
# Each row corresponds to a certain reduction of resolution. The best row is:
#```{R Quartet divergence - optimals, echo=FALSE}
which_suboptimal <- apply(divergences, 2, which.min) - 1
best_gc <- rownames(divergences)[which_suboptimal]
names(best_gc) <- ANALYSES_EN
knitr::kable(t(best_gc))
```
The only analyses to produce significantly different ($p = 0.01$) results from implied weights (at $k = 3$) or Bayesian are equal weights and implied weights with $k = 1$.
Both these approaches are significantly worse.
```{R Quartet-divergence-significance, include=FALSE}
treeBests <- vapply(clBootGcQuartets,
function (item) apply(
apply(item[, , TREES], 1, function (x) QuartetDivergence(t(x))),
1, max), double(length(TREES)))
tests <- apply(treeBests, 2,
function (column) t.test(column, treeBests[, 'k3'])$p.value)
ANALYSES_EN[tests < 0.01]
```
## Partition metric
We can run the same analysis counting partitions in place of quartets.
```{R CL-trees-partitions, echo=FALSE, message=FALSE, message=FALSE, cache=TRUE, fig.asp=1, fig.width=6, out.width='75%', fig.align='center'}
clInitializeTernarySplits()
clPlotAverageSplits(clBootGcPartitions, col=COL, pch=PCH, cex=1.1)
clPlotTheseAverageSplits(clMkvPartitions, col=COL['mk'], pch=PCH['mk'], cex=1.1)
arrows(sqrt(3/4) * 0.5, 0.5, sqrt(3/4) * 0.8, 0.5, length=0.1)
text (sqrt(3/4) * 0.65, 0.5, pos=3, 'Decreasing resolution', cex=0.8)
arrows(sqrt(3/4) * 0.98, 0.40, sqrt(3/4) * 0.98, 0.20, length=0.1)
text (sqrt(3/4) * 1.01, 0.30, pos=3, 'Increasing SD/MaxI', cex=0.8, srt=270)
AddLegend()
```
Under the partition metric, the most informative trees were found after collapsing nodes with a Bootstrap GC support value of:
```{R Partition-divergence-for-each-point, echo=FALSE, message=FALSE, fig.width=8, results='asis', fig.align='center'}
RFDistances <- function (item) {
colSums(aperm(item[, c('d1', 'd2', 'r1', 'r2'), TREES], c(2, 1, 3)))
}
divergences <- vapply(clBootGcPartitions,
function (item) rowMeans(RFDistances(item)),
double(dim(clBootGcPartitions[[1]])[1]))
rownames(divergences) <- seq(-100, 100, length.out=dim(clBootGcQuartets[[1]])[1])
# Each row corresponds to a certain reduction of resolution. The best row is:
which_suboptimal <- apply(divergences, 2, which.min) - 1
best_gc <- rownames(divergences)[which_suboptimal]
names(best_gc) <- ANALYSES_EN
knitr::kable(t(best_gc))
```
The partition metric advocates a greater loss of resolution than the quartet metric.
Otherwise, it too finds no statistically significant difference between the effectiveness of the methods, except again that equal weights, and implied weights with _k_ = 1, are significantly worse.
```{R Partition-significance, include=FALSE}
# Normalization would be against the same constant for all trees,
# so is unnecessary
treeBests <- vapply(clBootGcPartitions,
function (item) apply(RFDistances(item), 2, min),
double(length(TREES)))
tests <- apply(treeBests, 2,
function (column) t.test(column, treeBests[, 'k3'])$p.value)
ANALYSES_EN[tests < 0.01]
```
# Why small concavity constants are unsuitable
Of the implied weights concavity constants analysed above, $k = 1$ is strikingly
(and significantly) worse than other values. It is worth recalling the mathematical underpinning for implied weights [@Goloboff1993]:
$\textrm{Character penalty} = \frac{e}{e + k}$, where:
- $e$ is the number of additional steps;
- $k$ is the concavity constant.
The penalty can be normalized such that the first extra step in a character incurs a unit cost:
$\textrm{Normalized penalty} = (1 + k)\frac{e}{e + k}$.
As $k\rightarrow\infty$, the penalty
$\frac{e}{e + k}\rightarrow\frac{e}{k}$,
and the normalized penalty
$(1 + k)\frac{e}{e + k} \rightarrow e$.
This is to say, each subsequent step contributes the same amount to a tree's penalty; as $k\rightarrow \infty$, implied weights converges to equal weights.
At lower values of _k_, the penalty for extra steps decreases such that each subsequent additional step is penalized less than the previous one.
As $k \rightarrow 0$, the penalty for the first step converges to one, and the penalty for subsequent steps converges to zero. At this extreme, the optimal tree is the one that maximises the number of characters that are convex. (A convex character is one that can be plotted onto the tree with no additional steps; its derived states each exhibit a unique origin.) All characters that are not convex are uninformative, as the number of steps beyond the first is irrelevant to their total contribution to tree score. This situation corresponds to clique analysis [@Wilkinson1994], a method that is no longer advocated for use in phylogenetic reconstruction.
The value $k=1$ marks a significant point in the transition from parsimony analysis to clique analysis, because the highest cost that can be associated with a single character is less than twice the cost of a single extra step.
The cost associated with the first extra step is
$\frac{1}{1 + 1} = \frac{1}{2}$. As the number of extra steps increases($e\rightarrow \infty$), the penalty increases towards its maximum value of $\frac{\infty}{\infty + 1} \rightarrow 1$, i.e. just under twice the cost of the first step.
As such, given a pair of characters, a reconstruction that assigns infinitely many changes of one character, but no additional steps to the other, will be preferred to a reconstruction in which both characters undertake a single additional step.
---
```{R k-1-poor-choice, echo=FALSE, warning=FALSE, fig.width=7}
Fit <- function (e, k) (e / (e + k))
NormalizedFit <- function (e, k) (1 + k) * Fit(e, k)
e <- 0:8
plot(-1, -1,
xlab='Extra steps', xlim=c(0, max(e)),
ylab='Normalized penalty', ylim=c(0, 6)
)
kValues <- c(0, 1, 2, 3, 5, 10, 100)
kCols <- TreeSearch::brewer[[length(kValues)]]
abline(h=2, lty=3)
for (i in length(kValues):1) {
k <- kValues[i]
lines(NormalizedFit(e, k)~e, col=kCols[i])
x = 2 / (1 - (1 / k))
if (k > 1) lines(c(x, x), c(-1, 2), col=kCols[i], lty=3)
}
text(pos=2, paste('k =', kValues), x=8, y=NormalizedFit(8, kValues) + 0.25, col=kCols)
text(pos=2, paste('k =', max(kValues)), x=8, y=5.9, col=kCols[length(kCols)])
```
> Two characters with one extra step receive a total normalized penalty of two (dotted line).
At progressively smaller values of _k_, a single character must exhibit increasingly more steps before it receives the same penalty.
Once $k \leq 1$, no amount of steps in a single character will elicit a penalty equal to that which would be encountered if a second character undertook a single extra step.
---
Because non-convex characters are not entirely uninformative, this situation is not strictly equivalent to clique analysis. For example, trees that reconstruct fewer steps in a single non-convex character will still be preferred to those that reconstruct more steps in the same character; and a tree that imposes two extra steps on three characters receives the same penalty ($3 \times \frac{2}{2 + 1} = 2$) as one that imposes one extra step on four ($4 \times \frac{1}{1 + 1} = 2$).
Nevertheless, a value of $k=1$ places significantly more emphasis on maximising the number of convex characters than on minimizing the total number of steps in any given (non-convex) character, behaviour that is more characteristic of clique analysis than parsimony analysis. This supports the long-standing recommendation [@Goloboff1995] that
low values of $k$ should be avoided.
As there does not yet exist an objective method for selecting a single value of _k_ for parsimony analysis, it has been proposed that nodes are recovered by a range of concavity constants are likely to be correct [@Goloboff2008].
Studies that take a consensus of all trees found to be optimal under a range of concavity values (e.g. @Mirande2009) should consider discounting topologies that are only recovered under low values of _k_ (_cf_. [@Smith2015;@Zhang2016]).
# References