Fast and accurate estimation of the covariance between pairwise maximum likelihood distances

Preliminaries

Denote $A={\left\{x_i\right\}}_{i=1}^n$ a pairwise alignment consisting of n homologous i.i.d. character-pairs x_i (e.g., nucleotides, amino acids, or codons, but no insertion–deletions). The likelihood of having the two sequences in A separated by an evolutionary distance d is (Felsenstein, 1981) (1)${\displaystyle L\left(A|d\right)=\prod _{i=1}^{n}p\left(x_i,d\right),}$ where p(x_i, d) is the probability of the character-pair x_i at d. The ML estimator of the true distance d_t is (2)${\displaystyle \hat{d}=\underset{d\ge 0}{\mathrm{arg\; max}}L\left(A|d\right).}$ While for simple mechanistic substitution models the maximum can be expressed analytically, it is usually found numerically for empirical and complex mechanistic models using the Newton–Raphson method. Let I_n(d) denote the Fisher information for d, i.e., (3)${\displaystyle I_n\left(d\right)=-nE\left[\frac{{\partial }^2}{\partial d^2}log\hspace{0.17em}p\left(X,d\right)\right],}$ where X is a random variable of which the x_i are realizations. The asymptotic variance of $\hat{d}$ is provided by standard theory (e.g., Pawitan, 2001) (4)${\displaystyle V\left(\hat{d}\right)=I_n{\left(d_t\right)}^{-1}.}$ It can be estimated by evaluating the inverse of the Fisher information at $\hat{d}$ (hereafter ML variance) (5)${\displaystyle \hat{V}\left(\hat{d}\right)=I_n{\left(\hat{d}\right)}^{-1},}$ or, alternatively, from A and $\hat{d}$ by a sample average: (6)${\displaystyle {\hat{V}}_A\left(\hat{d}\right)=-\frac{1}{n}{\left({\left.\frac{1}{n}\sum _{i=1}^n\frac{{\partial }^2}{\partial d^2}log\hspace{0.17em}p\left(x_i,d\right)\right|}_{d=\hat{d}}\right)}^{-1}=-{\left({\left.\sum _{i=1}^n\frac{{\partial }^2}{\partial d^2}log\hspace{0.17em}p\left(x_i,d\right)\right|}_{d=\hat{d}}\right)}^{-1}.}$

Dessimoz & Gil (2008) distinguished three topological relations relevant for covariance estimation between any two pairwise distance estimates. First, the relation dependence, where two distances share some common evolution (e.g., d_ij and d_kl in Fig. 1). Second, the similar relation triplet, where two distances additionally share a sequence (e.g., d_ij and d_jk). Third, the case independence, where the distances are independent (e.g., d_ik and d_jl). Note that the second case can be conceptually reduced to the first (with j = l and δ_j = δ_l = 0). Further, we assume that mutation events on different edges of a tree are independent, thus, the distances in the independence-case have zero covariance. Therefore, our derivations will focus on the dependence-case without loss of generality.

General ML theory provides covariance estimates if all unknown parameters are estimated jointly. For instance, in ML tree reconstruction, the variance/covariance matrix can be estimated from the observed Fisher information matrix. However, the estimation of the pairwise distances considered here is done separately for each distance so that general ML theory can not be applied. Susko (2003) has derived an interesting estimator for the covariance of two distances ${\hat{d}}_{ij}$ and ${\hat{d}}_{kl}$ estimated from pairwise alignments which are induced by an MSA. It is a sample average based on the following expression (hereafter Susko-covariance) (7)${\displaystyle \mathrm{cov}\left({\hat{d}}_{ij},{\hat{d}}_{kl}\right)=n{V}_{ij}{V}_{kl}E\left[{\left.\frac{\partial }{\partial d}log\hspace{0.17em}{p}_{ij}\left(X,d\right)\right|}_{d=d_{t,ij}}\cdot {\left.\frac{\partial }{\partial d}log\hspace{0.17em}{p}_{kl}\left(X,d\right)\right|}_{d=d_{t,kl}}\right].}$ Here, d_t,ij and d_t,kl are the true distances, the random-variable X stands for quartets of homologous characters (as opposed to pairs in Eq. (3)), and p_ij(X, d) denotes the probability for the pair {i, j} in X at the distance d. The Susko-covariance has two limitations. First, it requires an MSA, i.e., it is not applicable to distances derived from OPAs (for a discussion see Dessimoz & Gil, 2008). Second, the O(n) computation time may become prohibitive in large scale studies. Alternatively, a nonparametric bootstrap can be used (Efron & Tibshirani, 1993), but it takes substantially longer computation times and requires an MSA too.

Previously, two estimators have been presented do tackle the limitations. They work with both MSAs and OPAs. The first method is a numerical approximation to the variance of the difference between two distances involving a common sequence. It runs in constant time with respect to the sequence-length and takes the following form (Dessimoz et al., 2006) (8)${\displaystyle {\hat{V}}_{NUM}\left({\hat{d}}_{ij}-{\hat{d}}_{jk}\right)=\frac{{\hat{d}}_{ik}^{{\varphi }_1}}{\hat{V}{\left({\hat{d}}_{ik}\right)}^{{\varphi }_2}}\cdot \frac{{\left[\hat{V}\left({\hat{d}}_{ij}\right)+\hat{V}\left({\hat{d}}_{jk}\right)\right]}^{{\varphi }_3}{\left[\hat{V}\left({\hat{d}}_{ij}\right)\hat{V}\left({\hat{d}}_{jk}\right)\left.\right)\right]}^{{\varphi }_4}}{{\left({\hat{d}}_{ij}+{\hat{d}}_{jk}\right)}^{{\varphi }_5}},}$ where the ϕ_i are parameters optimised a priori for a particular substitution model. This leads (through V(X−Y) = V(X) + V(Y)−2cov(X, Y)) to a fast covariance estimator for the triplet-case (hereafter referred to as triplet-covariance) (9)${\displaystyle \mathrm{cov}\left({\hat{d}}_{ij},{\hat{d}}_{jk}\right)=\frac{1}{2}\left[\hat{V}\left({\hat{d}}_{ij}\right)+\hat{V}\left({\hat{d}}_{jk}\right)-{\hat{V}}_{NUM}\left({\hat{d}}_{ij}-{\hat{d}}_{jk}\right)\right].}$ The second method is based on Susko’s theory and shares the linear time complexity (Dessimoz & Gil, 2008, hereafter anchor-covariance). It was specifically designed to bypass the problem of inconsistent homology inference between OPAs using the concept of anchors—a globally consistent subset of aligned character pairs. In this paper, we propose a fast and general approach we term branch-covariance. It is motivated by an analytic result obtained under the simple r-state symmetric model.

r-state symmetric model

To obtain analytic results we will work with the r-state symmetric model, also know as the N_r model (Neyman, 1971). It is a generalization of the Jukes–Cantor model (Jukes & Cantor, 1969), which has four states (r = 4), to r character-states. The N_r model assumes a uniform distribution of states at the root, and equal rates of transitions between any two distinct character states. The probability to observe a mutation after time t is (10)${\displaystyle p_m\left(t\right)=\beta \left(1-{\mathrm{e}}^{-\frac{\alpha t}{\beta }}\right),\beta =\frac{r-1}{r},}$ where α is the total rate of substitution. Thus, if two sequences are separated by t, the distance between them will be d = αt.

Because of the symmetries in the model, the number of differing sites I in a given pairwise alignment of length n is a sufficient statistics for the pairwise ML distance (11)${\displaystyle \hat{d}=-\beta ln\left(1-\frac{I}{n\beta }\right).}$

An estimator for the variance of $\hat{d}$ can be obtained by applying Eq. (4) derived from the likelihood function. Alternatively, for models estimating distances from proportions of differing sites, the variance can be approximated by the Delta technique. This has been done by Kimura & Ohta (1972) for r = 4 and generalized by Tajima & Nei (1984) to (12)${\displaystyle \hat{V}\left(\hat{d}\right)=\beta \left[\left(1-\beta \right){\mathrm{e}}^{\frac{2d}{\beta }}+\left(2\beta -1\right){\mathrm{e}}^{\frac{d}{\beta }}-\beta \right]/n.}$

We are going to use the Delta technique to derive an estimator for the covariance of two N_r ML distances d_ij = δ_i + δ_m + δ_j and d_kl = δ_k + δ_m + δ_l in the dependence-case (Fig. 1). Nei & Jin (1989) used an informal argument to propose the following expression with β = 3/4: (13)${\displaystyle \mathrm{cov}\left({\hat{d}}_{ij},{\hat{d}}_{kl}\right)=\beta \left[\left(1-\beta \right){\mathrm{e}}^{\frac{2{\delta }_m}{\beta }}+\left(2\beta -1\right){\mathrm{e}}^{\frac{{\delta }_m}{\beta }}-\beta \right]/n.}$ The equation originates from the assumption, that the covariance of two distance estimates with an underlying shared path length δ_m, is formally equivalent to the variance (Eq. (12)) of an ML estimate of a pairwise distance δ_m. Indeed, Bulmer (1989) presented a proof for β = 3/4 for the triplet-case and conjectured that Eq. (13) with j = l and δ_j = δ_l = 0 is true for any β. To the best of our knowledge Eq. (13) has not been proven yet in its most general form, i.e., for the dependence-case and any β.

We will first compute (14)${\displaystyle \mathrm{cov}\left(I_{ij},I_{kl}\right)=n\mathrm{cov}\left(S_{ij},S_{kl}\right)=nE\left[S_{ij}{S}_{kl}\right]-nE\left[S_{ij}\right]E\left[S_{kl}\right],}$ where S_ij is a random variable indicating whether sequences i and j are identical (S_ij = 0) or different (S_ij = 1) at a particular site. Subsequently, we will apply the Delta method to obtain Eq. (13) from cov(I_ij, I_kl). We start by noting that (15)${\displaystyle E\left[S_{ij}\right]=Pr\left(S_{ij}=1\right)=p_m\left(d_{ij}\right).}$ Therefore, the problem reduces to computing (16)${\displaystyle E\left[S_{ij}{S}_{kl}\right]=Pr\left(S_{ij}=1\wedge S_{kl}=1\right).}$ In the following we will represent the quartet from Fig. 1 by the symbol , where the terminal node of the upper left branch corresponds to i. Furthermore, we are going to mark a branch with ○ if a particular site in the evolving sequence is in a different state at the endpoints of the branch. As an example, we look at the pattern . Here, the site in question changed its state on the branches leading to {j, k, l} but did not change its state on the branch leading to i and on the middle branch. A particular labeling of the nodes with characters from the alphabet of N_r for the given pattern has probability

For this pattern, there are r(r−1)³ possible labellings, of which only r(r−1)²(r−2) satisfy S_ij = 1∧S_kl = 1. Symmetrical mutation patterns, like for example have the same number of labelings leading to S_ij = 1∧S_kl = 1 but different mutation probabilities. We consider now all the patterns (grouped by symmetry) and corresponding labelings for the desired condition:

Using Maple (script in Supplemental Materials) we find that the expression simplifies to (19)${\displaystyle E\left[S_{ij}{S}_{kl}\right]=\beta \left(1-\beta {\mathrm{e}}^{-\frac{d_{ij}}{\beta }}-\beta {\mathrm{e}}^{-\frac{d_{kl}}{\beta }}+\left(1-\beta \right){\mathrm{e}}^{-\frac{d_{ij}+d_{kl}-2{\delta }_m}{\beta }}+\left(2\beta -1\right){\mathrm{e}}^{-\frac{d_{ij}+d_{kl}-{\delta }_m}{\beta }}\right).}$ Plugging Eqs. (19) and (15) in (14) we obtain (20)${\displaystyle \mathrm{cov}\left(I_{ij},I_{kl}\right)=n\beta \left(\left(1-\beta \right){\mathrm{e}}^{-\frac{d_{ij}+d_{kl}-2{\delta }_m}{\beta }}+\left(2\beta -1\right){\mathrm{e}}^{-\frac{d_{ij}+d_{kl}-{\delta }_m}{\beta }}-\beta {\mathrm{e}}^{-\frac{d_{ij}+d_{kl}}{\beta }}\right).}$ We turn to the Delta method. The function ${\hat{d}}_{ij}\left(I_{ij}\right)$ can be approximated by a first-order Taylor series around E[I_ij] = np_m(d_ij) (21)${\displaystyle {\hat{d}}_{ij}^{\star }\left(I_{ij}\right)={\hat{d}}_{ij}\left(E\left[I_{ij}\right]\right)+{\hat{d}}_{ij}^{\prime }\left(E\left[I_{ij}\right]\right)\left(I_{ij}-E\left[I_{ij}\right]\right),}$ where (22)${\displaystyle {\hat{d}}_{ij}^{\prime }\left(E\left[I_{ij}\right]\right)={\left.\frac{\partial {\hat{d}}_{ij}\left(I_{ij}\right)}{\partial I_{ij}}\right|}_{I_{ij}=E\left[I_{ij}\right]}={\left.{\left(n-\frac{I_{ij}}{\beta }\right)}^{-1}\right|}_{I_{ij}=E\left[I_{ij}\right]}=\frac{{\mathrm{e}}^{\frac{d_{ij}}{\beta }}}{n}.}$ The covariance of ${\hat{d}}_{ij}$ and ${\hat{d}}_{kl}$ is asymptotically equal to the covariance of ${\hat{d}}_{ij}^{\star }$ and ${\hat{d}}_{kl}^{\star }$ (23)${\displaystyle \mathrm{cov}\left({\hat{d}}_{ij}\left(I_{ij}\right),{\hat{d}}_{kl}\left(I_{kl}\right)\right)\sim \mathrm{cov}\left({\hat{d}}_{ij}^{\star }\left(I_{ij}\right),{\hat{d}}_{kl}^{\star }\left(I_{kl}\right)\left.\right)\right)}$ (24)${\displaystyle ={\hat{d}}_{ij}^{\prime }\left(E\left[I_{ij}\right]\right){\hat{d}}_{kl}^{\prime }\left(E\left[I_{kl}\right]\right)\mathrm{cov}\left(I_{ij},I_{kl}\right)}$ (25)${\displaystyle =\beta \left[\left(1-\beta \right){\mathrm{e}}^{\frac{2{\delta }_m}{\beta }}+\left(2\beta -1\right){\mathrm{e}}^{\frac{{\delta }_m}{\beta }}-\beta \right]/n.\square }$ We provide a Maple program implementing all steps of the proof in Supplemental Materials.

Covariance under general Markov models

We have derived under N_r that the covariance of two distance estimates with an underlying shared path length δ_m is asymptotically equal to the variance of an ML estimate of a true pairwise distance δ_m, i.e., (26)${\displaystyle \mathrm{cov}\left({\hat{d}}_{ij},{\hat{d}}_{kl}\right)\sim V\left(\hat{d}|d_t={\delta }_m\right).}$ We conjecture now that the relationship holds for general substitution models and apply it to derive a covariance estimator. To this end, we first discuss how the covariance can be written in terms of a rate matrix Q, an equilibrium frequency vector π (which together fully specify a substitution model), and δ_m. In the next section, we will show how δ_m can be estimated from the input distances by the method of weighted least squares (WLS).

According to the conjecture and Eq. (4) we express the desired covariance by (27)${\displaystyle \mathrm{cov}\left({\hat{d}}_{ij},{\hat{d}}_{kl}\right)=-\frac{1}{n}E{\left[{\left.\frac{{\partial }^2}{\partial d^2}log\hspace{0.17em}p\left(X,d\right)\right|}_{d={\delta }_m}\right]}^{-1}.}$ The expected value expression can be written as (28)${\displaystyle \sum _{\forall \left(u,v\right)}{\left[p\left(x_{uv},d\right)\frac{{\partial }^2}{\partial d^2}log\hspace{0.17em}p\left(x_{uv},d\right)\right]}_{d={\delta }_m}=\sum _{\forall \left(u,v\right)}{\left[\frac{{\partial }^2}{\partial d^2}p\left(x_{uv},d\right)-\frac{{\left(\frac{\partial }{\partial d}p\left(x_{uv},d\right)\right)}^2}{p\left(x_{uv},d\right)}\right]}_{d={\delta }_m},}$ where the summation goes over all possible character pairs. In terms of a rate matrix Q and an equilibrium frequency vector π the probability of a character pair (k = 0) and the derivatives (k > 0) of the probability with respect to the distance are (29)${\displaystyle \frac{{\partial }^{\left(k\right)}}{\partial d^{\left(k\right)}}p\left(x_{uv},d\right)={\pi }_u{\left[Q^k{\mathrm{e}}^{Qd}\right]}_{uv}.}$ Plugging Eqs. (28) and (29) in (27) we obtain (30)${\displaystyle \mathrm{cov}\left({\hat{d}}_{ij},{\hat{d}}_{kl}\right)=-\frac{1}{n}{\left[\sum _{\forall \left(u,v\right)}{\pi }_u\left({\left[Q^2{\mathrm{e}}^{Q{\delta }_m}\right]}_{uv}-\frac{{\left[Q{\mathrm{e}}^{Q{\delta }_m}\right]}_{uv}^2}{{\left[{\mathrm{e}}^{Q{\delta }_m}\right]}_{uv}}\right)\right]}^{-1}.}$ A plot of this expression as a function of δ_m can be found in Fig. S1. A covariance estimator is obtained by substituting δ_m in Eq. (30) by its WLS estimate, which we derive in the next section. To save computation time, we can discretise the distance space in the relevant range to some desired level of accuracy, precompute the expected value expressions, and store them in a hash table.

Topological relation and path length

Equation (30) expresses the covariance in the dependence case as a function of the shared path length δ_m. To obtain a covariance estimator, we determine first whether the two distances in question are dependent and, in case they are, estimate δ_m. We will do that by WLS using the six pairwise distance estimates $\left\{{\hat{d}}_{uv}\right\}$ between the four sequences i, j, k, l and their variances {v_uv}. The sequences can be related by three topological configurations: ${\displaystyle T_1:\left(\left(i,k\right),\left(j,l\right)\right),T_2:\left(\left(i,l\right),\left(j,k\right)\right),T_3:\left(\left(i,j\right),\left(k,l\right)\right),}$ where T₁, T₂ map to the dependence case and T₃ corresponds to the independence case. An argument set out in Supplemental Materials shows that the weighted sum of squares (S) for each of the topologies can be expressed in a simple form which is fast to compute: (31)${\displaystyle \begin{array}{l}S\left(T_1\right)=\frac{{\hat{d}}_{ij}+{\hat{d}}_{kl}-{\hat{d}}_{il}-{\hat{d}}_{jk}}{v_{ij}+v_{kl}+v_{il}+v_{jk}},S\left(T_2\right)=\frac{{\hat{d}}_{ij}+{\hat{d}}_{kl}-{\hat{d}}_{ik}-{\hat{d}}_{jl}}{v_{ij}+v_{kl}+v_{ik}+v_{jl}},\hfill \\ S\left(T_3\right)=\frac{{\hat{d}}_{ik}+{\hat{d}}_{jl}-{\hat{d}}_{il}-{\hat{d}}_{jk}}{v_{ik}+v_{jl}+v_{il}+v_{jk}}.\hfill \end{array}}$ The best fitting topology is determined by argmin_Ti{S(T_i)}. If this results in T₃ the desired covariance is zero, otherwise we need to estimate δ_m. The WLS estimates are (see Supplemental Materials): (32)${\displaystyle 2{\hat{\delta }}_m\left(T_1\right)=\frac{\left({\hat{d}}_{ij}+{\hat{d}}_{kl}\right)\left(v_{il}+v_{jk}\right)+\left({\hat{d}}_{il}+{\hat{d}}_{jk}\right)\left(v_{ij}+v_{kl}\right)}{v_{ij}+v_{kl}+v_{il}+v_{jk}}-\left({\hat{d}}_{ik}+{\hat{d}}_{jl}\right),}$ (33)${\displaystyle 2{\hat{\delta }}_m\left(T_2\right)=\frac{\left({\hat{d}}_{ij}+{\hat{d}}_{kl}\right)\left(v_{ik}+v_{jl}\right)+\left({\hat{d}}_{ik}+{\hat{d}}_{jl}\right)\left(v_{ij}+v_{kl}\right)}{v_{ij}+v_{kl}+v_{ik}+v_{jl}}-\left({\hat{d}}_{il}+{\hat{d}}_{jk}\right).}$ Since these are the estimators for unconstrained WLS they can result in negative values, in which case we estimate the covariance to be zero.

Simulation settings

The performance of the various covariance estimators was evaluated with the same simulation approach as in a previous study (Dessimoz & Gil, 2008). Specifically, 100 random quartets were sampled from a tree of life on 352 species. The tree was inferred by the LeastSquaresTree function (Gil & Gonnet, 2009) included in the Darwin package (Gonnet et al., 2000) using pairwise distance and variance data from the OMA project (Dessimoz et al., 2005). A uniformly distributed U(0.5, 2) expansion/contraction factor was applied on each quartet to also explore extremer regions of the branch-length space, while preserving the relative branch-length structure of the original tree.

For each dilated model quartet 10,000 times three random amino-acid sequences of length m = {200, 500, 800} were generated and mutated along the quartet assuming the GCB substitution model (Gonnet, Cohen & Benner, 1992). The entire simulation procedure was run twice, once without any insertion–deletions to produce ungapped alignments to test the methods under the true models (i.e., without the effect of alignment errors), and once introducing gaps of Zipfian distributed length (Benner, Cohen & Gonnet, 1993). The gapped sequences were aligned by global pairwise dynamic programming with the Align function from Darwin to obtain OPAs. ML pairwise distances were estimated with the EstimatePam function. The sample variance and covariance over the 10,000 samples (hereafter Monte Carlo-variance and Monte Carlo-covariance) served as a reference values, as they are unbiased estimators of the true values.

To test the asymptotic conjecture (Eq. (26)) the ungapped simulation was repeated with very long sequences (10,000 amino-acids). Then, the length δ_m of the middle branch of each model quartet was extracted. Subsequently, 10,000 random sequences were generated, and mutated each with a distance δ_m according to the GCB model. Finally, the ML distance and ML variance between each resulting pair of sequences, and the Monte Carlo-variance were computed.

Evaluation of basic components of branch-covariance

We have tested the validity of the conjecture in Eq. (26), which was derived under the N_r model, and the accuracy of the ML variance (Eq. (5)) in a simulation under the GCB model with long sequences (10,000 amino-acids).

A plot of the Monte Carlo-variance of the ML estimate of a pairwise distance δ_m versus the Monte Carlo-covariance between the ML estimates of two pairwise distances in the dependence case with a shared path of length δ_m corroborates the conjecture and suggests that the result is valid in general (Fig 2A). The branch-covariance relies on the ML variance to approximate the true variance. A comparison with the Monte Carlo-variance shows that it underestimates the variance for large samples and with correct alignments (Fig 2B). Therefore, we expected the branch-covariance to inherit the negative bias.

Evaluation of estimators

To evaluate and compare the various estimators under ideal conditions we have carried out two types of simulations: one without introducing indel events (ungapped simulation), and one with indel events (gapped simulation).

The gapped simulation did not require explicit alignment; the MSA is trivial without gaps. This situation corresponds to the model assumed by the Susko-, branch-, and triplet-covariance. Note that in practice with real data, there are usually gaps in an MSA estimate. In this case, the distance and covariance estimators are applied assuming the MSA to be correct. Gaps are either excluded from the analysis, or treated as unknown characters.

The ungapped simulation required alignment which was carried out by OPA. Here, the Susko-covariance can not be applied; it requires an MSA. The setting corresponds to the model assumed by the anchor-covariance, which is an adaptation of the Susko-covariance, specifically designed to bypass the problem of inconsistent homology inference between OPAs. The triplet- and branch-covariance are applicable to OPAs, as these estimators operate directly on the distances and not on the alignments (but the assumption of a gap-free and correct underlying alignment is violated).

Gapped simulation

The branch-, anchor-, and triplet-covariance have been tested on distances derived from OPAs.

The branch-covariance has a lower bias than the anchor-covariance for all three topological relations (Fig. 4). In the dependence case the differences are minor, though the branch-covariance has a consistently higher correlation with the Monte Carlo covariance. A big difference is visible for the two other topological relations (independence and triplet), where the anchor-covariance’s bias is up to twice the branch-covariance’s. The superiority of the branch-covariance is also reflected by the MSEs (Table 2). The triplet-covariance has a greater bias than the branch-covariance for sequences of length 200; for length 500 the two estimators have quantitatively a similar bias, but in opposite directions; and for length 800 the triplet-covariance has clearly a smaller systematic error.

Table 2:

Gapped simulation.

Average MSE, median MSE (italic), and average of ratios between MSE of branch-covariance and MSE of anchor-covariance (bold). A ratio of 0.22, for instance, indicates that the MSE of the branch-covariance was—on average—22% of the MSE of the anchor-covariance. Dependence (D), triplet (T) and independence (I) case for sequence lengths 200, 500 and 800 amino-acids.

		200			500			800
		D	T	I	D	T	I	D	T	I
Branch	Avg	36.19	170.04	3.63	3.99	27.23	0.21	1.38	9.68	0.06
	Median	19.13	110.60	1.23	1.95	6.26	0.03	0.48	1.88	0.00
Anchor	Avg	159.07	242.56	95.79	11.60	95.62	7.68	3.19	42.87	2.09
	Median	102.13	170.27	59.02	6.66	12.30	3.51	1.67	4.55	0.92
	Avg ratio	0.22	0.67	0.06	0.22	0.36	0.06	0.20	0.34	0.04

DOI: 10.7717/peerj.583/table-2

Note that the anchor- and triplet-covariance have been evaluated under the same simulation conditions in a previous work (Dessimoz & Gil, 2008). They have been reproduced to evaluate the branch-covariance. The results on the triplet- and anchor-covariance reported here are in agreement with our previous results.