Review: A Simulation Comparison of Phylogeny Algorithms under Equal and Unequal Evolutionary Rates¶
Citation
- Kuhner, M. K., & Felsenstein, J. (1994). A simulation comparison of phylogeny algorithms under equal and unequal evolutionary rates. Molecular Biology and Evolution, 11(3), 459–468.
- DOI
Abstract¶
Using simulated data, we compared five methods of phylogenetic tree estimation: parsimony, compatibility, maximum likelihood, Fitch-Margoliash, and neighbor joining. For each combination of substitution rates and sequence length, 100 data sets were generated for each of 50 trees, for a total of 5,000 replications per condition. Accuracy was measured by two measures of the distance between the true tree and the estimate of the tree, one measure sensitive to accuracy of branch lengths and the other not.
Kuhner and Felsenstein (1994) introduce the primary evaluation framework used in this paper: a simulation study comparing five phylogenetic reconstruction methods across 5,000 replicates per condition. Two accuracy measures are used: the Robinson-Foulds \( dT \) score for topological accuracy, and the branch score \( Bs \) for joint topology and branch length accuracy. A clade-indexed form of the branch score is the tree-distance metric behind Hifuku's two-dimensional chart.
The Branch Score¶
The branch score is defined over the full set of possible bipartitions (splits) of the leaf set. Every edge in an unrooted binary tree induces a bipartition of the leaves into two groups, one on each side. Let \( P_1, P_2, \ldots, P_N \) enumerate all possible such bipartitions, and let \( b_i \) be the branch length assigned to bipartition \( P_i \) in tree \( T \), with \( b_i = 0 \) if that bipartition does not exist in \( T \). The branch score between two trees \( T \) and \( T' \) is:
This is the squared Euclidean distance between the trees' branch-length vectors in the ambient space of all bipartitions. Branches that appear in one tree but not the other are compared against zero, so the score penalizes both branch length error and topological disagreement simultaneously.
Key properties noted by Kuhner and Felsenstein:
- \( Bs = 0 \) for identical trees; it increases as the match worsens.
- \( Bs \) can be 0 for two non-identical topologies if all discordant branches have length 0 (a degenerate case).
- The branch score is sensitive to the absolute size of branch lengths, so scores cannot be directly compared across trees with different substitution rates.
- \( dT \) is a special case of \( Bs \) where all \( b_i \in \{0, 1\} \).
For rooted trees the same construction applies using ordered partitions (clades): each non-root node defines a clade (the leaf set of its subtree), and branches are indexed by clade rather than bipartition.
Hifuku's Tree-Distance Metric¶
Hifuku uses the clade branch score (CBS), the square root of the clade-indexed branch score, as its tree-distance metric, \( s(T, T') = \sqrt{Bs(T, T')} \). The map is not built from a single such distance. The chart is fixed by three anchor trees, and a tree is located by its CBS distances \( s_1, s_2, s_3 \) to the three anchors, which give its barycentric coordinates in the two-dimensional chart. The metric also reports a saturation score \( B_\mathrm{frac} = 2B/(A + C) \in [0, 1] \) used to reject anchor pairs that share too few splits.
Each CBS distance is maintained incrementally. Writing \( s^2 = A - 2B + C \) for the distance to one anchor:
where the sums are over all non-root nodes and clades are the indexing scheme (the rooted variant of equation 1). \( C \) is constant per anchor. \( A \) and \( B \) are running scalars updated after each MCMC move:
| Move | A update | B update |
|---|---|---|
| branch_slide (leaf \( k \), \( \Delta l \)) | \( A \mathrel{+}= (l_{\mathrm{new}}^2 - l_{\mathrm{old}}^2) \) | \( B \mathrel{+}= \Delta l \cdot r_k \) |
| branch_scale (\( \times c \)) | \( A \mathrel{\times}= c^2 \) | \( B \mathrel{\times}= c \) |
| NNI (node \( u \)) | unchanged | \( B \mathrel{+}= l_u(r_{\mathrm{new}} - r_{\mathrm{old}}) \); update \( r_u \) |
Here \( r_k \) is the anchor branch length for node \( k \)'s clade (looked up once at initialization and refreshed after topology-changing moves). All three updates are \( O(1) \) in branch-length operations; the NNI update also requires an \( O(n) \) DFS to reidentify the clade at the rearranged node. The sampler carries one \( (A, B) \) pair per anchor, so a single fused traversal updates all three distances a step needs.
Simulation Results (summary)¶
Maximum likelihood was the most accurate method overall across all rate and sequence length conditions. The distance methods (Fitch-Margoliash, neighbor joining) were slightly inferior under equal rates but similarly accurate under unequal rates per branch; both were unbiased. Parsimony and compatibility showed bias when rates varied across branches, tending to incorrectly group taxa with long branches (long-branch attraction, Felsenstein (1981)). When rates varied across sites, all methods showed inaccuracy and bias.
The branch score results (Tables 1–4, part C) showed that the three branch-length-estimating methods (ML, F-M, NJ) performed similarly in absolute accuracy, consistent with the metric being dominated by the longest branches.
Relationship to Robinson-Foulds and BHV¶
\( Bs \) sits between the purely topological RF distance (Robinson & Foulds (1981)) and the BHV geodesic (Billera et al. (2001)). RF ignores branch lengths entirely; BHV measures the true geodesic through tree space, which may pass through multiple topology orthants. The branch score is the Euclidean distance in the ambient bipartition space, a Euclidean approximation that is cheap to compute and maintain incrementally. That makes it a good choice for the CBS distances that fix the chart coordinates.