iGLASS: An Improvement to the GLASS Method for Estimating Species Trees from Gene Trees

3.1. Derivation of \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$f_{V_{AB}} (v \mid {\bf k}_A , {\bf k}_B)$$\end{document} for one locus.

Consider a single locus and let k_A and k_B denote the numbers of lineages from taxa A and B remaining at the divergence time τ_AB. The quantity \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$f_{V_{AB}} (v \mid k_A , k_B)$$\end{document} is then the distribution of the time to the first interspecific coalescence at the locus, measuring from time τ_AB.

To derive \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$f_{V_{AB}} (v \mid k_A , k_B)$$\end{document} , recall that the time T_i until i lineages coalesce to i − 1 lineages is exponentially distributed with mean \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$1 / {i \choose 2}$$\end{document} Thus, if k = k_A + k_B lineages remain at the divergence time, and if the first interspecific coalescence occurs on the Mth coalescence past the divergence time, then the waiting time V _AB until this coalescence can be expressed as the summation \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$V_{AB} = \sum_{i = 1}^M T_{k - (i - 1)}$$\end{document} .

The location M in the sequence of coalescences of the first interspecific coalescence is itself a random variable and hence, V _AB has a Coxian distribution (Ross, 2007) with probability density function given by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*}f_{V_{AB}} (v \mid k_A , k_B) = \sum_{m = 1}^{k - 1} Pr (M = m) \sum_{i = 1}^m c_{i , m} \gamma_{i} e^{- \gamma_i v}. \tag{6}\end{align*}\end{document}

In Equation (6), \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$c_{i , m} = \prod_{j = 1 , j \neq i}^{m} \gamma_j / (\gamma_j - \gamma_i)$$\end{document} , where \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\gamma_i = {k - (i - 1) \choose 2}$$\end{document} is the parameter of the ith waiting time T_k−(i−1). For m = 1, we define c_i,m to be unity.

The distribution in Equation (6) was derived by Takahata (1989). In Takahata's result, the probability Pr(M = m) is obtained recursively; however, it is possible to derive a closed-form solution. Rosenberg (2003) derived a closed-form expression that is equivalent to the cumulative distribution function of M. In Appendix B, we derive the closed-form of the probability mass function of M from Equation A8 of Rosenberg (2003). We obtain \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} Pr (M = m) = \frac {I_{k - m , 1}} {2^{m - 1} k_A k_B I_{k , 1}} \sum_{\eta = \max \{0 , m - k_B \}}^{\min \{m - 1 , k_A - 1 \}} {m - 1 \choose \eta} {k_A}^2_{[\eta + 1]} {k_B}^2_{[m - \eta]} \tag {7}\end{align*}\end{document}

whenever m ≤ k_A + k_B − 1, where \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$k_{[i]} = \frac {k!} {(k - i) !}$$\end{document} , and where, as in Rosenberg (2003), \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$I_{k , m} = \frac {k! (k - 1) !} {2^{k - m} m! (m - 1) !}$$\end{document} is the number of ways in which k lineages can coalesce down to m lineages. Plugging expression (7) into (6) gives the formula for the distribution of V _AB in the case of one locus.

3.2. Derivation of \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$f_{V_{AB}} (v \mid {\bf k}_A , {\bf k}_B)$$\end{document} for L loci.

We now extend formula (6) to multiple loci. Let \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$V_{AB}^{(\ell)}$$\end{document} be the random variable describing the time to the first interspecific coalescence at locus \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\ell \ ({\ell} = 1 , \ldots , L)$$\end{document} . We assume that all loci are independent, conditional on the species tree and its parameter values. Therefore, the cumulative distribution function \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$F_{V_{AB}} (v \mid {\bf k}_A , {\bf k}_B)$$\end{document} of the minimum interspecific coalescence time \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$V_{AB} = \min_{\ell} V_{AB}^{(\ell)}$$\end{document} is given by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*}F_{V_{AB}} (v \mid {\bf k}_A , {\bf k}_B) & = 1 - Pr \bigg (V_{AB}^{(\ell)} \geq v , \forall \ell = 1 , \ldots , L \mid {\bf k}_A , {\bf k}_B \bigg) \\ & = 1 - \prod_{{\ell} = 1}^L Pr \left(V_{AB}^{(\ell)} \geq v \mid k_{A_{\ell}} , k_{B_{\ell}} \right) . & (8) \end{align*}\end{document}

Here, \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$Pr (V_{AB}^{(\ell)} \geq v \mid k_{A_{\ell}} , k_{B_{\ell}})$$\end{document} is given by integrating Equation (6): \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*}Pr (V_{AB}^{(\ell)} \geq v \mid k_{A_{\ell}} , k_{B_{\ell}}) & = \int_{t = v}^{\infty} f_{V_{AB}^{(\ell)}} (t \mid k_{A_{\ell}} , k_{B_{\ell}}) dt \\ & = \sum_{m_{\ell} = 1}^{k_{\ell} - 1} Pr (M_{\ell} = m_{\ell}) \sum_{i_{\ell} = 1}^{m_{\ell}} c_{i_{\ell} , m_{\ell}} e^{- \gamma_{i_{\ell}} v} , & (9) \end{align*}\end{document}

where k_ℓ = k_Aℓ + k_Bℓ is the total number of lineages remaining at the divergence time at locus ℓ. Plugging (9) into (8) and differentiating gives the density function \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$f_{V_{AB}} (v \mid {\bf k}_A , {\bf k}_B)$$\end{document} : \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} f_{V_{AB}} (v \mid {\bf k}_A , {\bf k}_B) & = \frac {d} {dv} \Bigg(1 - \prod_{\ell = 1}^L Pr (V_{AB}^{(\ell)} \geq v \mid k_{A_{\ell}} , k_{B_{\ell}}) \Bigg) \\ & = - \sum_{{\ell} = 1}^L \left[\prod_{j = 1 \atop j \neq \ell}^L Pr (V_{AB}^{(j)} \geq v \mid k_{A_{j}} , k_{B_{j}}) \right] \frac {d} {dv} Pr \left(V_{AB}^{(\ell)} \geq v \mid k_{A_{\ell}} , k_{B_{\ell}} \right) \\ & = \sum_{{\ell} = 1}^L \left[\sum_{m_1 = 1}^{k_1 - 1} \sum_{i_1 = 1}^{m_1} \cdots \sum_{m_L = 1}^{k_L - 1} \sum_{i_L = 1}^{m_L} Pr (M_1 = m_1) \cdots Pr (M_L = m_L)\right. \\ & \left. \quad \qquad \qquad \qquad \qquad \times \ c_{i_1 , m_1} \cdots c_{i_L , m_L} \gamma_{i_{\ell}} e^{- (\gamma_{i_1} + \cdots + \gamma_{i_L}) v} \right] \\ & = \sum_{m_1 = 1}^{k_1 - 1} \sum_{i_1 = 1}^{m_1} \cdots \sum_{m_L = 1}^{k_L - 1} \sum_{i_L = 1}^{m_L} Pr (M_1 = m_1) \cdots Pr (M_L = m_L) \\ & \qquad \qquad \quad \qquad \times c_{i_1 , m_1} \cdots c_{i_L , m_L} (\gamma_{i_1} + \cdots + \gamma_{i_L}) e^{- (\gamma_{i_1} + \cdots + \gamma_{i_L}) v} . & (10)\end{align*}\end{document}

In the last equality, we have brought the outer summation inside.

3.3. Closed-form expressions for \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$f_{V_{AB}} (v)$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB}]$$\end{document} .

Closed-form expressions for \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$f_{V_{AB}} (v)$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB}]$$\end{document} can now be computed using Equation (10). The unconditional density \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$f_{V_{AB}} (v)$$\end{document} is given by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*}f_{V_{AB}} (v) = \sum_{{\bf k}_A , {\bf k}_B} f_{V_{AB}} (v \mid {\bf k}_A , {\bf k}_B) \prod_{\ell = 1}^L h_{n_{A_{\ell}} , k_{A_{\ell}}} (\tau_{AB}; N_A) h_{n_{B_{\ell}} , k_{B_{\ell}}} (\tau_{AB}; N_B) , \tag{11}\end{align*}\end{document}

where the summation at a given locus \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\ell \ (\ell = 1 , \ldots , L)$$\end{document} in a given taxon (A or B) ranges from 1 to the number of sampled lineages at that locus in that taxon. The conditional expected value of V _AB for a collection of L loci is obtained by integrating Equation (10). This gives \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} E_{\tau_{AB}} [V_{AB} & \mid {\bf k}_A , {\bf k}_B] \\ & = \int_{t = 0}^{\infty} t f_{V_{AB}} (t \mid {\bf k}_A , {\bf k}_B) dt \\ & = \sum_{m_1 = 1}^{k_1 - 1} \sum_{i_1 = 1}^{m_1} \cdots \sum_{m_L = 1}^{k_L - 1} \sum_{i_L = 1}^{m_L} Pr (M_1 = m_1) \cdots Pr (M_L = m_L) \\ & \qquad \qquad \qquad \qquad \qquad \times c_{i_1 , m_1} \cdots c_{i_L , m_L} \frac {1} {(\gamma_{i_1} + \cdots + \gamma_{i_L})} . & (12) \end{align*}\end{document}

The unconditional expected value \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB}]$$\end{document} can be computed by plugging either Equation (12) or the recursive Equation (5) into Equation (3), thereby completing the derivation of \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB}]$$\end{document} .

Thus, to obtain the iGLASS estimate from the GLASS estimate \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{t}_{AB}$$\end{document} , we evaluate Equation (2), where \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*}g (\tau_{AB}) = \tau_{AB} + \sum_{{\bf k}_A , {\bf k}_B} E_{\tau_{AB}} [V_{AB} \mid {\bf k}_A , {\bf k}_B] \prod_{\ell = 1}^L h_{n_{A_{\ell}} , k_{A_{\ell}}} (\tau_{AB}; N_A) h_{n_{B_{\ell}} , k_{B_{\ell}}} (\tau_{AB}; N_B) , \tag{13}\end{align*}\end{document}

and where \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB} \mid {\bf k}_A , {\bf k}_B]$$\end{document} is given by either Equation (12) or Equation (5). The product \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\prod_{\ell = 1}^L h_{n_{A_{\ell}} , k_{A_{\ell}}} (\tau_{AB}; N_A) h_{n_{B_{\ell}} , k_{B_{\ell}}} (\tau_{AB}; N_B)$$\end{document} is a polynomial in \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$e^{- \tau_{AB}}$$\end{document} and thus, the inverse \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$g^{- 1} (\hat{t}_{AB})$$\end{document} must be evaluated numerically. The iGLASS estimate of the species tree is then constructed by applying an appropriately chosen clustering method to the distance matrix of pairwise iGLASS time estimates. We discuss the choice of clustering method in Section 7.

4.1. An approximation to \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$f_{V_{AB}} (v \mid {\bf k}_A , {\bf k}_B)$$\end{document} for one locus.

As before, consider two taxa A and B. Let k_A and k_B be the numbers of lineages, not necessarily integers, that enter the ancestral population at the divergence time from taxa A and B, respectively. For the remainder of this derivation, it will simplify the notation if we measure time from a reference point at the divergence time τ_AB, rather than from the present. Thus, we take \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{k_A}^t$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{k_B}^t$$\end{document} to be the numbers of lineages remaining at time t from taxa A and B, counting from the divergence time.

Although \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{k_A}^t$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{k_B}^t$$\end{document} are deterministic quantities representing the expected numbers of lineages from taxa A and B, we continue to assume that the interactions between lineages are random. We assume that, in a small time interval [t, t + Δt], a coalescent event occurs with rate \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\varphi_{k_A}^t + \varphi_{k_B}^t \choose 2}$$\end{document} , given that no interspecific coalescence has occurred by time t. In addition, given that a coalescent event occurs in the interval [t, t + Δt], we approximate the probability that it is interspecific by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$2 \frac {\varphi_{k_A}^t} {\varphi_{k_A}^t + \varphi_{k_B}^t} \frac {\varphi_{k_B}^t} {\varphi_{k_A}^t + \varphi_{k_B}^t - 1}$$\end{document} , the conditional probability that a coalescence at time t involves one lineage from taxon A and one lineage from taxon B if the numbers of lineages are integer-valued. Thus, letting \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal I}_{a , b}$$\end{document} be the event that an interspecific coalescence occurs in the interval [a, b], letting \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal I}_{a , b}^c$$\end{document} be the event that an interspecific coalescence does not occur in the interval [a, b], and letting \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal C}_{a , b}$$\end{document} be the event that a coalescence of any kind occurs in the interval [a, b], we find that \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} Pr ({\cal I}_{t , t + \Delta t} \mid {\cal I}_{0 , t}^c) & = Pr ({\cal I}_{t , t + \Delta t} \mid {\cal C}_{t , t + \Delta t} , {\cal I}_{0 , t}^c) Pr ({\cal C}_{t , t + \Delta t} \mid {\cal I}_{0 , t}^c) \\ & \approx 2 \frac {\varphi_{k_A}^t} {\varphi_{k_A}^t + \varphi_{k_B}^t} \frac {\varphi_{k_B}^t} {\varphi_{k_A}^t + \varphi_{k_B}^t - 1} {\varphi_{k_A}^t + \varphi_{k_B}^t \choose 2} \Delta t \\ & = \varphi_{k_A}^t \varphi_{k_B}^t \Delta t. \end{align*}\end{document}

Hence, the approximate probability that an interspecific coalescence does not occur in the interval [t, t + Δt], given that none has occurred more recently than time t, is \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} Pr ({\cal I}_{t , t + \Delta t}^c \mid {\cal I}_{0 , t}^c) \approx 1 - \varphi_{k_A}^t \varphi_{k_B}^t \Delta t \approx e^{- \varphi_{k_A}^t \varphi_{k_B}^t \Delta t}. \end{align*}\end{document}

The probability that no interspecific coalescence occurs in the interval [0, t] can be approximated by the probability that no interspecific coalescence occurs in any of J small intervals of length Δt = t/J: \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*}Pr ({\cal I}_{0 , t}^c) \approx \prod_{j = 0}^{J - 1} e^{- \varphi_{k_A}^{j \Delta t} \varphi_{k_B}^{j \Delta t} \Delta t} = \exp \left\{- \sum_{j = 0}^{J - 1} \varphi_{k_A}^{j \Delta t} \varphi_{k_B}^{j \Delta t} \Delta t \right\} .\end{align*}\end{document}

Thus, as J → ∞ we have Δt → 0, and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*}Pr ({\cal I}_{0 , t}^c) \rightarrow \exp \left\{- \int_{z = 0}^{t} \varphi_{k_A}^{z} \varphi_{k_B}^{z} dz \right\} . \tag{17}\end{align*}\end{document}

We now generalize this result to the case of many loci.

4.2. An approximation to \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${f_{V_{AB}}} (v \mid {\bf k}_{A} , {\bf k}_{B})$$\end{document} for L loci.

Let \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{k_{A_{\ell}}}^t$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{k_{A_{\ell}}}^t$$\end{document} be the deterministic approximations to the numbers of lineages remaining at time t from taxa A and B at locus ℓ. Then the probability that no interspecific coalescence occurs in any one of L independent loci in the interval [0, t] is approximately \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*}Pr (V_{AB} \geq t \mid {\bf k}_A , {\bf k}_B) & \approx \prod_{{\ell} = 1}^L \exp \left\{- \int_{z = 0}^{t} \varphi_{k_{A_{\ell}}}^z \varphi_{k_{B_{\ell}}}^z dz \right\} \\ & = \exp \left\{- \sum_{{\ell} = 1}^L \int_{z = 0}^{t} \varphi_{k_{A_{\ell}}}^z \varphi_{k_{B_{\ell}}}^z dz \right\} . & (18) \end{align*}\end{document}

4.3. The approximate iGLASS correction.

To get the expected time to the first interspecific coalescence at some locus, the approximation to Equation (12), we integrate: \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*}E_{\tau_{AB}} [V_{AB} \mid {\bf k}_A , {\bf k}_B] & = \int_{t = 0}^{\infty} Pr \left(V_{AB} \geq t \mid {\bf k}_A , {\bf k}_B \right) dt \\ & \approx \int_{t = 0}^{\infty} \exp \left\{- \sum_{{\ell} = 1}^L \int_{z = 0}^{t} \varphi_{k_{A_{\ell}}}^z \varphi_{k_{B_{\ell}}}^z dz \right\} dt. & (19) \end{align*}\end{document}

If we assume that the number of lineages remaining at the divergence time is the expected number of lineages at this time, then the approximate iGLASS correction, the approximation to Equation (3), is obtained by making the substitutions \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\bf k}_{A_\ell} = E_{\tau_{AB}} [{\bf K}_{A_\ell}]$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\bf k}_{B_\ell} = E_{\tau_{AB}} [{\bf K}_{B_\ell}]$$\end{document} into Equation (19): \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*}E_{\tau_{AB}} [V_{AB}] & \approx E_{\tau_{AB}} [V_{AB} \mid E_{\tau_{AB}} [{\bf K}_A] , E_{\tau_{AB}} [{\bf K}_B]] \\ & \approx \int_{t = 0}^{\infty} \exp \left\{- \sum_{{\ell} = 1}^L \int_{z = 0}^{t} \varphi_{E_{\tau_{AB}} [K_{A_{\ell}}]}^z \varphi_{E_{\tau_{AB}} [K_{B_{\ell}}]}^z dz \right\} dt \\ & \equiv \tilde{E}_{\tau_{AB}} [V_{AB}] . & (20) \end{align*}\end{document}

This approximate expression is much faster to evaluate than Equation (3) because it does not require a sum over all possible values of \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$k_{k_{A_{\ell}}}$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$k_{k_{B_{\ell}}}$$\end{document} .

Because the values obtained from the approximation (Equation 20) differ from those obtained from the exact solution (Equation 3), we modify our definition of the iGLASS estimator (Equation 2) accordingly. We now define the function \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde{g} (\tau_{AB}) = \tau_{AB} + \tilde{E}_{\tau_{AB}} [V_{AB}]$$\end{document} , and we define the approximate iGLASS estimator \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde{\tau}_{AB}$$\end{document} to be \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} \begin{split} \tilde{\tau}_{AB} = \begin{cases} \tilde{g}^{- 1} (\hat{t}_{AB}^*) , \quad {\rm if} \ \tilde{E}_0 [V_{AB}] \leq \hat{t}_{AB}^* \\ 0 , \qquad \qquad \ {\rm otherwise}.\end{cases}\end{split} \tag{21}\end{align*}\end{document}

As before, the approximate iGLASS estimate \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{S}$$\end{document} of the species tree is then constructed by applying any suitable clustering method to the pairwise approximate iGLASS estimates.

Although Equation (20) is an approximation, it can produce values that are remarkably close to the exact expectations. Figure 3 shows the exact survival function Pr(V _AB  ≥ v) of V _AB (Equation 9) and the approximate survival function (Equation 18) for the case of one locus. From Figure 3, it can be seen that the approximation is exact when one lineage is sampled per taxon, because the expected number of lineages used in the approximation is always equal to one, the true number of lineages.

OPEN IN VIEWER

FIG. 3.

Approximate survival function (Equation 18) (red, dashed) and exact survival function (Equation 9) (blue) of the quantity \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$V_{AB} = \hat{t}_{AB} - \tau_{AB}$$\end{document} for one locus, conditional on the numbers of lineages k_A and k_B remaining at the divergence time from each taxon. Pr(V ≥ v|K_A = k_A, K_B = k_B) is the probability that the GLASS estimate exceeds the divergence time τ_AB by more than v coalescent units. In order from top to bottom, the numbers of lineages that were used to generate the curves are (k_A, k_B) = (1,1), (1,2), (2,2), (2,3), (3,3), (3,5), (5,5), (5,7), where k_A is the number of lineages remaining in taxon A at the divergence time and k_B is the corresponding number of lineages remaining in taxon B. For the top curve, one lineage is sampled from each taxon and the approximation is exact.

For larger numbers of sampled lineages, as the time v is increased the approximation becomes slightly worse and then improves again. This result is a consequence of the behavior of the variability in the number of lineages over time. For small v, with very high probability the number of lineages is close to the number that were initially sampled, and the variance in the number of lineages is small. For intermediate v, greater variation exists in the number of lineages, and the approximation of the stochastic process of coalescence as a deterministic process is less appropriate. Finally, for large v, the number of lineages is equal to one with high probability, and the variance is again small. Thus, the expectation \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_n^v$$\end{document} is a better approximation to the number of lineages for small and large v.

In practice, the approximate iGLASS correction (Equation 20) differs only slightly from the exact iGLASS correction, except in the case of a single locus (Fig. 4). Therefore, in our implementation of the iGLASS correction, we use the approximation (Equation 20), except in the case of a single locus, for which it is fast to compute the exact correction.

OPEN IN VIEWER

FIG. 4.

The difference between the approximate iGLASS correction (Equation 19) and the exact iGLASS correction (Equation 3). Each pixel in the heatmap shows the difference \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde{E}_{\tau_{AB}} [V_{AB}] - E_{\tau_{AB}} [V_{AB}]$$\end{document} for a given divergence time τ_AB, a given number of lineages sampled per taxon, and a given number of loci. Within each block corresponding to a number of loci, the numbers of lineages sampled from each taxon at each locus are, from left to right, (n_A, n_B) = (1,1), (1,2), (1,3), (1,4), (2,2), (2,3), (2,4), (3,3), (3,4), and (4,4), where n_A is the number of lineages sampled from taxon A and n_B is the number sampled from taxon B.

6.1. Exact and approximate iGLASS are consistent estimators of pairwise divergence times.

As we show in Theorem (D.1) in Appendix D, the GLASS method is a consistent estimator of pairwise divergence times. The exact and approximate iGLASS estimators (Equations 2 and 21) approach the GLASS estimator asymptotically in such a way that they are also consistent. We now prove this result.

6.2. Exact and approximate iGLASS are consistent estimators of the species tree topology.

We now show that both the exact and approximate iGLASS methods are consistent estimators of the species tree topology whenever the clustering procedure applied to the estimates of pairwise divergence times has certain desirable properties. Let \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal D}$$\end{document} be a distance matrix whose elements are pairwise distances between taxa in the species tree S computed according to some distance measure. Let \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat {\cal D}$$\end{document} be an estimate of \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal D}$$\end{document} . Let ∥A∥_∞ denote the magnitude of the largest element in a matrix A. Following Atteson (1999), we give the following definition.

Theorem 8.1

Consider two taxa A and B. If a single lineage is sampled from each taxon at each locus ℓ \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$(\ell = 1 , \ldots , L)$$\end{document} , then the estimator defined by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{\tau}_{AB} = g^{- 1} (\hat{t}_{AB})$$\end{document} for all \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{t}_{AB} \in {\mathbb R}$$\end{document} is an unbiased estimator of the divergence time τ_AB.

Proof

Let \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$k_{A_{\ell}}$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$k_{B_{\ell}}$$\end{document} be the numbers of lineages remaining at locus ℓ \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$(\ell = 1 , \ldots , L)$$\end{document} from taxa A and B at the divergence time. When one lineage is sampled from each taxon at each locus, \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$k_{A_{\ell}}$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$k_{B_{\ell}}$$\end{document} equal one for all \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\ell = 1 \ldots , L$$\end{document} . Therefore, letting 1 be the vector of length L with all entries equal to 1, Equation (5) gives E[V _AB |K _A  = 1, K _B  = 1] = 1/L, and Equation (3) simplifies to \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB}] = 1 / L$$\end{document} . The function g(τ_AB) is then given by g(τ_AB) = τ_AB + 1/L, and its inverse by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$g^{- 1} (\hat{t}_{AB}) = \hat t_{AB} - 1 / L$$\end{document} . Hence, g⁻¹(t) is defined for all \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$t \in {\mathbb R}$$\end{document} and it is linear. Thus, by the linearity of the expectation operator, \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [g^{- 1} (\hat{t}_{AB})] = g^{- 1} (E_{\tau_{AB}} [\hat{t}_{AB}]) = g^{- 1} (\tau_{AB} + E_{\tau_{AB}} [V_{AB}]) = g^{- 1} (g (\tau_{AB})) = {\tau_{AB}}.$$\end{document}    ▪

This result implies that the iGLASS estimator defined by Equation (2) is also unbiased for most values of τ_AB whenever one lineage is sampled per taxon. Specifically, as we have assumed that gene trees are inferred with certainty, the GLASS estimate \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{t}_{AB}$$\end{document} always exceeds the true divergence time τ_AB. Therefore, when one lineage is sampled per taxon at each locus and the true divergence time is greater than or equal to 1/L, it follows that \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{t}_{AB} \ge 1 / L$$\end{document} and the iGLASS estimator is defined by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{\tau}_{AB} = g^{- 1} (\hat{t}_{AB}) = \hat{t}_{AB} - 1 / L$$\end{document} . Thus, by Theorem 8.1, the iGLASS estimator will be unbiased in this case.

Note that when more than one lineage is sampled from either taxon, the probability \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\prod_{\ell = 1}^L h_{n_{A_{\ell}} , k_{A_{\ell}}} (\tau_{AB}; N_A) h_{n_{B_{\ell}} , k_{B_{\ell}}} (\tau_{AB}; N_B)$$\end{document} in Equation (3) contains terms of the form \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$e^{- \tau_{AB}}$$\end{document} , and thus, the quantity \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB}]$$\end{document} is no longer linear in τ_AB. In this case, \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$g^{- 1} (\hat{t}_{AB})$$\end{document} is not linear in \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{t}_{AB}$$\end{document} and therefore, we cannot use the relationship \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [g^{- 1} (\hat{t}_{AB})] = \tau_{AB}$$\end{document} when more than one lineage is sampled per taxon. However, as we will see from simulations, the bias is still very small.

9. Comparison of Methods

We used simulations to compare the performance of iGLASS to that of GLASS, evaluating each method on the basis of bias and mean squared error (MSE). We first evaluated the methods for estimating pairwise divergence times, and we then applied them to larger trees.

9.1. Simulations

We simulated gene trees under the multispecies coalescent model for various species trees S, for various numbers of loci, and for various numbers of lineages sampled per taxon. In all simulations, all population sizes were equal to the same value N across the branches of the species tree.

To simulate a gene tree from a given species tree, we used a method similar to that of Rosenberg and Feldman (2002). Let branch i refer to the branch above node i in the species tree. Let t_i be the time at node i, and let \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\bar{t}_{i}$$\end{document} be the time at the node ancestral to node i. Here, we extend our numbering to external branches, with t_i = 0 when i corresponds to a leaf node.

Let n_i be the number of lineages entering branch i at time t_i. If branch i is internal, then n_i is the sum of the numbers of lineages entering from its left and right daughter branches. If branch i is external, then n_i is equal to the number of lineages sampled from the corresponding taxon.

In each branch i, with the enumeration beginning with the external branches and proceeding towards the root in such a way that daughter branches have lower numbers than their parental branches, we first sampled the waiting time T_ni until the first coalescence from an exponential distribution with mean \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$1 / {n_i \choose 2}$$\end{document} . If the sampled time T_ni exceeded \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\bar{t}_i - t_i$$\end{document} , then we let the set of lineages exiting branch i equal the set that entered. Otherwise, we chose two lineages at random without replacement and allowed them to coalesce. We continued in this way, at each coalescence sampling the time to the next coalescence from an exponential distribution with mean \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$1 / {q \choose 2}$$\end{document} , where q was the number of lineages remaining after the previous coalescence, until the sum of waiting times in the branch exceeded \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\bar{t}_i - t_i$$\end{document} . The set of lineages remaining after the last coalescence to occur within branch i was then merged into the set of lineages entering its ancestral branch, along with the set of lineages entering from its sister branch, and the process was repeated in the ancestral branch. Simulations were run until all lineages coalesced to a single lineage. For trees with more than two taxa, the simulations were carried out using the software program ms (Hudson, 2002).

Let \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$n_{X_{\ell}}$$\end{document} denote the number of lineages sampled from taxon X at locus ℓ \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$(\ell = 1 , \ldots , L)$$\end{document} . For a given species tree S together with a set of parameters consisting of a number of loci L and numbers of lineages \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{n_{X_{\ell}} \}_{X \in {{\cal S}}}$$\end{document} , we first sampled r independent sets \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal L}_j$$\end{document} of L gene trees \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$(j = 1 , \ldots , r)$$\end{document} . For each set \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal L}_j$$\end{document} , we computed the GLASS estimate \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{t}^{(j)}_{AB}$$\end{document} for all pairs of species \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$A , B \in {\cal S}$$\end{document} using the GLASS algorithm (Section 1), without applying the single-linkage clustering step. From each observation \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{t}^{(j)}_{AB}$$\end{document} , we then computed an observation \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{\tau}^{(j)}_{AB}$$\end{document} of the exact iGLASS estimate, and an observation \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde{\tau}^{(j)}_{AB}$$\end{document} of the approximate iGLASS estimate. We thus obtained the sets of pairwise estimates \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{\hat{t}^{(j)}_{AB} \}_{A , B \in {\cal S}}$$\end{document} , \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{\hat{\tau}^{(j)}_{AB} \}_{A , B \in {\cal S}}$$\end{document} , and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{\tilde{\tau}^{(j)}_{AB} \}_{A , B \in {\cal S}}$$\end{document} , for each set of gene trees \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal L}_j \ (j = 1 , \ldots , r)$$\end{document} . For species trees with more than two taxa, only \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{\hat{t}^{(j)}_{AB} \}_{A , B \in {\cal S}}$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{\tilde{\tau}^{(j)}_{AB} \}_{A , B \in {\cal S}}$$\end{document} were computed.

For each set \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{\hat{t}^{(j)}_{AB} \}_{A , B \in {\cal S}} \ (j = 1 , \ldots , r)$$\end{document} , we computed the GLASS estimate \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{S}^{(j)}$$\end{document} of the species tree by single-linkage clustering, and for each internal node i in this estimated species tree, we estimated the height \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{t}_i^{(j)}$$\end{document} of the node i by the distance between the two clusters combined on the step of the clustering method that produced the node. The clustering procedure was omitted for trees with two taxa because the estimates \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\hat{t}^{(j)}_{AB}} \ (j = 1 , \ldots , r)$$\end{document} already provide estimates of the divergence time τ_AB. We then compared each estimated node height \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{t}^{(j)}_{i}$$\end{document} to its true value t_i, and we computed the average difference \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${{\cal B}}^{(j)} (\hat{t}) = \sum^{\mid {\cal S} \mid - 1}_{i = 1} (\hat{t}^{(j)}_{i} - t_{i}) / (\mid {{\cal S}} \mid - 1)$$\end{document} and the average squared difference \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${{\cal M}}^{(j)} (\hat{t}) = \sum^{\mid {\cal S} \mid - 1}_{i = 1} (\hat{t}^{(j)}_{i} - t_{i})^2 / (\mid {{\cal S}} \mid - 1)$$\end{document} .

Average bias in the GLASS method was estimated by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat {bias}}_{avg} (\hat {t}) = \frac {1} {r} \sum_{j = 1}^r {{\cal B}}^{(j)} (\hat {t})$$\end{document} , and average MSE by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat {MSE}}_{avg} (\hat {t}) = \frac {1} {r} \sum_{j = 1}^r {{\cal M}}^{(j)} (\hat {t})$$\end{document} . The average bias and MSE in the exact and approximate iGLASS methods were estimated by the same procedure (using single-linkage clustering), but using the times \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{\hat{\tau}^{(j)}_{AB} \}_{A , B \in {{\cal S}};j = 1 , \ldots , r}$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{\tilde{\tau}^{(j)}_{AB} \}_{A , B \in {{\cal S}};j = 1 , \ldots , r}$$\end{document} .

We denote the average bias and MSE in the exact iGLASS method by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat{bias}}_{avg} (\hat{\tau})$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat{MSE}}_{avg} (\hat{\tau})$$\end{document} , and we denote the average bias and MSE in the approximate iGLASS method by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat{bias}}_{avg} (\tilde{\tau})$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat{MSE}}_{avg} (\tilde{\tau})$$\end{document} .

9.2. Estimating pairwise divergence times.

To evaluate the performance of the three methods for estimating pairwise divergence times, we simulated gene trees under the multispecies coalescent from a species tree with two taxa, for various values of the parameters τ_AB, L, \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{n_{A_{\ell}} \}_{\ell = 1}^L$$\end{document} , and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{n_{B_{\ell}} \}_{\ell = 1}^L$$\end{document} , and for r = 50,000 replicates. In varying the parameters \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{n_{A_{\ell}} \}_{\ell = 1}^L$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{n_{B_{\ell}} \}_{\ell = 1}^L$$\end{document} , we maintained the relationships \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$n_{A_{\ell}} = n_A$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$n_{B_{\ell}} = n_B$$\end{document} for all ℓ.

We considered values of 1, 5, 10, and 50 for L. However, because the exact iGLASS estimate is difficult to compute in the case of both multiple loci and large numbers of lineages, only the GLASS estimate and approximate iGLASS estimate were computed when both the number of loci and the number of lineages were large.

9.2.1. Bias

Figure 5 indicates that especially for small divergence times, the bias in the GLASS estimate can be large relative to the divergence time. Whenever a single lineage is sampled from each taxon, the bias in the GLASS method is 1/L in coalescent units of N generations, regardless of the divergence time. One lineage always remains at the divergence time from each taxon at each locus, and therefore, the expected time to the first interspecific coalescence is the expectation of the minimum of L independent exponentially distributed random variables, each with a mean of one coalescent time unit. For example, in a haploid population with an effective size of N = 10, 000, if the GLASS estimate is based on a single lineage sampled from each population at each of 20 loci, then the bias in the GLASS estimate is 10,000/20 = 500 generations.

OPEN IN VIEWER

FIG. 5.

Comparison of bias and mean squared error for the GLASS, exact iGLASS, and approximate iGLASS methods for two taxa and one locus. All values were computed using 50,000 simulation replicates. In each of the fourteen small heatmap panels, the divergence time between two taxa A and B is given in coalescent units on the y-axis. In each heatmap, the divergence times are, from top to bottom, τ_AB = 0, 0.1, 0.5, 1, 2, and 4 coalescent units. In each heatmap, the numbers of lineages sampled from each taxon are given on the x-axis in the format (n_A, n_B), where n_A is the number of lineages sampled from taxon A, and n_B is the number sampled from taxon B. From left to right, the numbers of lineages in each column are (n_A, n_B) = (1,1), (1,3), (1,5), (1,10), (1,15), (1,20), (3,3), (3,5), (3,10), (3,15), (3,20), (5,5), (5,10), (5,15), (5,20), (10,10), (10,15), (10,20), (15,15), (15,20), (20,20).

Although sampling multiple lineages from each population can greatly reduce the bias for low divergence times, it does not reduce the bias for larger divergence times. As noted by Mossel and Roch (2010), when τ_AB is measured in units of generations, the probability that a single lineage remains at the top of the branch corresponding to taxon A is bounded below by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$1 - 3 e^{- \tau_{AB} / N_A}$$\end{document} and the probability that a single lineage remains at the top of the branch corresponding to taxon B is bounded below by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$1 - 3 e^{- \tau_{AB} / N_B}$$\end{document} (Tavaré, 1984). This bound can be made arbitrarily close to one by increasing the divergence time and, as the divergence time increases, the GLASS estimate approaches the value of the GLASS estimate when one lineage is sampled per taxon, or 1/L coalescent units.

To compare the estimated bias in the exact and approximate iGLASS methods to the estimated bias in the GLASS method, we computed the ratios \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat{bias}} (\hat{\tau}_{AB}) / {\widehat{bias}} (\hat{t}_{AB})$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat{bias}} (\tilde{\tau}_{AB}) / {\widehat{bias}} (\hat{t}_{AB})$$\end{document} (Fig. 5). For most values of the divergence time, the bias in the approximate iGLASS method is negligible compared to the bias in the GLASS method; although it is considerably larger in magnitude for small values of τ_AB, the bias ratio continues to be less than 1. The bias is not entirely negligible in this case because we define the exact and approximate iGLASS estimates to be zero whenever the GLASS estimate is lower than its smallest possible expected time (Equations 2 and 21). Thus, when the GLASS estimate is small, instead of subtracting a positive quantity from the GLASS estimate to produce the iGLASS estimate, we estimate the divergence time to be zero, resulting in an iGLASS estimate (exact or approximate) that is biased upwards. This truncation prevents the iGLASS estimators from completely eliminating the bias, but it also leads to a decrease in variance, which ultimately leads to a lower mean squared error at these divergence times. The decrease in MSE due to lower variance can be seen by the yellow bars across the tops of the MSE graphs in Figure 5.

9.2.2. Mean squared error

The ratios \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat{MSE}} (\hat{\tau}_{AB}) / {\widehat{MSE}} (\hat{t}_{AB})$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat{MSE}} (\tilde{\tau}_{AB}) / {\widehat{MSE}} (\hat{t}_{AB})$$\end{document} are shown in Figure 5 for various values of τ_AB, n_A, and n_B. From these plots, we can see that \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat{MSE}} (\hat{\tau}_{AB}) / {\widehat{MSE}} (\hat{t}_{AB})$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat{MSE}} (\tilde{\tau}_{AB}) / {\widehat{MSE}} (\hat{t}_{AB})$$\end{document} are roughly 1/2, and that they appear to approach 1/2 as τ_AB increases.

To see why this is reasonable, consider the case in which a single lineage is sampled per taxon at each locus. In this case, the “overshoot” in the GLASS estimate, \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$V_{AB} = \hat{t}_{AB} - \tau_{AB}$$\end{document} , is distributed exponentially with mean 1/L. Thus, the bias in the GLASS estimator is \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB}] = 1 / L$$\end{document} , its variance is Var(V_AB) = 1/L², and its MSE is \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\rm MSE} (\hat{t}_{AB}) = 2 / L^2$$\end{document} . The variance in the GLASS estimator then accounts for half of the mean squared error when one lineage is sampled per taxon.

When one lineage is sampled per taxon, the iGLASS correction to the GLASS estimator is computed by subtracting a constant quantity 1/L from the GLASS estimate, except when \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{t}_{AB}$$\end{document} is in the region \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\in [0 , 1 / L)$$\end{document} , which decreases in size as L → ∞. Thus, the variance of the (exact or approximate) iGLASS estimator is nearly equal to the variance of the GLASS estimator. As Theorem 8.1 indicates, when a single lineage is sampled per taxon, the iGLASS estimator is almost unbiased. Thus, when a single lineage is sampled per taxon, the MSE in the (exact or approximate) iGLASS estimator is approximately equal to the variance in the GLASS estimator, which is half the MSE in the GLASS estimator. Because \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$k_{A_{\ell}}$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$k_{B_{\ell}}$$\end{document} approach one in probability as τ_AB → ∞, we expect that \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\rm MSE} (\hat{\tau}_{AB})$$\end{document} will approach \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\rm MSE} (\hat{t}_{AB}) / 2$$\end{document} as τ_AB increases to infinity.

9.3. Exact versus approximate iGLASS

In the majority of our simulations, we have used the approximate iGLASS correction rather than the exact method because the exact correction is difficult to compute. However, consider the panels in the first row of Figure 5 that correspond to the case of one locus. It can be seen that the bias and MSE in the approximate iGLASS method are very similar to the bias and MSE in the exact iGLASS method. This result indicates that making the approximation \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB}] \approx {\tilde E}_{\tau_{AB}} [V_{AB}]$$\end{document} (Equation 20) has little effect on the performance of the iGLASS estimator in the case of one locus. Because Figure 4 indicates that the approximation is least accurate in the case of a single locus, the similarity of the bias and MSE for the exact and approximate methods in the case of one locus suggests that making the approximation \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB}] \approx {\tilde E}_{\tau_{AB}} [V_{AB}]$$\end{document} generally has little effect on the performance of the iGLASS method relative to that of GLASS.

9.4. iGLASS for larger trees

Figure 6 shows the ratios \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat{bias}}_{avg} (\tilde{\tau}_{AB}) / {\widehat{bias}}_{avg} (\hat{t}_{AB})$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\widehat{MSE}}_{avg} (\tilde{\tau}_{AB}) / {\widehat{MSE}}_{avg} (\hat{t}_{AB})$$\end{document} computed over r = 50,000 replicates for two different five-taxon species trees similar to those used by Liu et al. (2010) to evaluate the performance of the GLASS method. One internal branch of the tree is short enough that the most likely gene tree given the species tree does not have the topology of the true tree. In other words, the tree is in the anomaly zone of Degnan and Rosenberg (2006).

OPEN IN VIEWER

FIG. 6.

Comparison of mean squared error and bias in the approximate iGLASS and GLASS methods for two five-taxon species trees used in Liu et al. (2010) to evaluate the GLASS method. In Newick format, the tree in (a) and (b), a caterpillar, is given by ((((E:0.5, D:0.5):0.025, C:0.525):0.025, B:0.55):10.0, A:10.55). The tree in (c) and (d), another 5-taxon caterpillar, is ((((E:0.5, D:0.5):0.2, C:0.7):1, B:1.7):10.0, A:11.7). The first tree is in the anomaly zone; the second tree is not. All values were computed using 50,000 replicates. The clustering method applied to the approximate iGLASS estimates was single-linkage. (a) and (c): The ratio \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\widehat{bias}_{avg} (\tilde{\tau}_{AB}) / \widehat{bias}_{avg} (\hat{t}_{AB})$$\end{document} . (b) and (d): The ratio \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\widehat{MSE}_{avg} (\tilde{\tau}_{AB}) / \widehat{MSE}_{avg} (\hat{t}_{AB})$$\end{document} .

From Figure 6, we see that the average bias in the iGLASS estimate is often considerably less than that of the GLASS estimate. The improvement in the bias is best for small numbers of loci and decreases as the number of loci increases. However, the bias in the GLASS method itself decreases quickly as the number of loci is increased.

Note that although the iGLASS correction improves the bias and MSE in the estimates of species tree node heights, it does not improve the accuracy in estimating topologies. For both species trees ((((E:0.5, D:0.5):0.025, C:0.525):0.025, B:0.55):10.0, A:10.55) and ((((E:0.5, D:0.5):0.2, C:0.7):1, B:1.7):10.0, A:11.7) that we considered, the GLASS and iGLASS methods have identical accuracies for estimating the topology. However, for the case in which only one lineage is sampled at only one locus, the GLASS method has slightly higher accuracy for inferring the topology (Fig. 7).

OPEN IN VIEWER

FIG. 7.

The fraction of tree topologies correctly inferred by the approximate iGLASS and GLASS methods for two different five-taxon species trees. The tree in (a) is the same tree considered in Figure 6a,b. The tree in (b) is the same tree considered in Figure 6c,d. Plots show the fraction of 50,000 simulated data sets in which the species tree topology was correctly inferred by GLASS and approximate iGLASS.

The reduction in accuracy for the case of one lineage and one locus was due to the fact that in this case, the iGLASS method estimated more than one pairwise divergence time in the species tree to be zero, resulting in ties that were sometimes resolved to produce a clade that was not on the true species tree. Multiple estimates of zero were produced in this case because the smallest possible expected value \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_0 [\hat t_{AB}]$$\end{document} of the GLASS estimate for a pair of taxa was equal to one, which was greater than at least two of the node heights in each tree that we considered (0.5 and 0.525 for the first tree, and 0.5 and 0.7 for the second tree).

For all other parameter values we considered, \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_0 [\hat t_{AB}]$$\end{document} was smaller than all of the node heights in either tree, and no estimates of zero were produced. For example, when two lineages were sampled per taxon, the smallest possible expected GLASS estimate was E₀[T_AB] = 0.39, which is smaller than 0.5, the smallest node height in either tree. Similarly, when one lineage was sampled per taxon at 5 loci, the smallest expected interspecific coalescence time was E₀[T_AB] = 0.2. Consequently, for all cases we considered except for the case of one sampled lineage per taxon at one locus, the accuracy of the iGLASS method for estimating topologies was the same as that of the GLASS method.

10. Discussion

For two taxa, A and B, we have derived a closed-form expression for the distribution of \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$V_{AB} = \min_{\ell} {\hat{t}}^{(\ell)}_{AB} - \tau_{AB}$$\end{document} , the waiting time to the first interspecific coalescence across L loci, measuring from the divergence time τ_AB. By computing the expectation E_τAB[V_AB], we constructed a correction to the GLASS estimator \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\hat{t}}_{AB} = {\rm min}_{\ell} \ \hat{t}^{(\ell)}_{AB}$$\end{document} of pairwise divergence times, which we call the iGLASS estimator.

Maruvka et al. (2011) have demonstrated that simple functions of time t in a population of constant size can provide useful deterministic continuous approximations of the number of lineages remaining at time t under the standard coalescent model. By approximating the number of lineages at time t by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_x^t$$\end{document} , the expected number of lineages remaining at time t when x lineages are sampled at time t = 0 and when x is not necessarily an integer, we derived an approximation \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde{\tau}_{AB}$$\end{document} to the exact iGLASS estimator \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{\tau}_{AB}$$\end{document} that is faster to compute than the exact value, and that is quite accurate even when the number of lineages is small.

Through simulations, we have shown that the exact and approximate iGLASS estimators reduce the bias in the GLASS estimates of pairwise divergence times. In addition, the exact iGLASS estimator \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{\tau}_{AB}$$\end{document} and its approximation \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde{\tau}_{AB}$$\end{document} generally reduce the mean squared error in the GLASS estimate of pairwise divergence times by approximately one half. This reduction accords with a theoretical prediction in the case in which a single lineage is sampled per taxon.

In our simulations, the accuracy of the iGLASS method for estimating topologies was similar to that of the GLASS method. In the case in which one lineage was sampled per taxon at one locus, iGLASS was slightly poorer, due to the fact that iGLASS produces divergence time estimates of zero whenever the GLASS estimate is smaller than its smallest possible expected value, \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_0 [\hat{t}_{AB}]$$\end{document} . Because \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_0 [\hat{t}_{AB}]$$\end{document} is smaller when the number of sampled lineages or loci is larger, divergence time estimates of zero are less likely when more lineages or loci are sampled. Therefore, the accuracy of the topology estimates produced by iGLASS are likely to be the same as those produced by GLASS whenever sufficiently many lineages or loci are sampled.

We have shown that the exact iGLASS estimator and its approximation are consistent estimators of the pairwise divergence time between a pair of taxa. Further, we have proven that applying any clustering method with nonzero L_∞-radius to the pairwise iGLASS estimates produces a statistically consistent estimator of the species tree topology.

Assuming that gene trees have been correctly inferred, the bias in the GLASS method itself decreases to zero quickly as the number of loci increases. Thus, our correction produces the greatest improvement when information is available for relatively few loci. As we have seen, however, the approximate iGLASS correction is fast to compute even for large numbers of loci, requiring only \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$O (n^2 L \mid {{\cal S}} \mid + L \mid {{\cal S}} \mid^2)$$\end{document} operations for a given level of precision, compared to \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$O (n^2 L \mid {{\cal S}} \mid^2 + \mid {{\cal S}} \mid^3)$$\end{document} operations for GLASS. Consequently, our new estimator provides a method that is reasonable to implement even when information is available at many loci.

11. Appendix A

A recursive formula for \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document}$$E_{\tau_{AB}} [V_{AB} \mid {\rm k}_A , {\rm k}_B ]$$\end{document}

In Appendix A, we derive Equation (5), the expected value of the difference \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$V_{AB} = \hat{t}_{AB} - \tau_{AB}$$\end{document} , conditional on the numbers of lineages that remain at each locus at the divergence time. Let C_ℓ be the event that the first coalescence occurs in locus \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\ell} \ ({\ell} = 1 , \ldots , L)$$\end{document} . We then recursively consider what happens on the next coalescent event: \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} E_{\tau_{AB}} [V_{AB} \mid {\bf k}_A , {\bf k}_B] \\ & = \sum_{{\ell} = 1}^L E_{\tau_{AB}} [V_{AB} \mid {\bf k}_A , {\bf k}_B , C_{\ell}] Pr (C_{\ell}) \\ & = \sum_{{\ell} = 1}^L \Bigg[ \frac {1} {\sum_{\lambda = 1}^L {k_{A_{\lambda}} + k_{B_{\lambda}} \choose 2}} + E_{\tau_{AB}} [V_{AB} \mid {\bf k}_A - {\bf e}_{\ell} , {\bf k}_B] \frac{{k_{A_{\ell}} \choose 2}} {{k_{A_{\ell}} + k_{B_{\ell}} \choose 2}} \Bigg. \\ & \Bigg. \quad + E_{\tau_{AB}} [V_{AB} \mid {\bf k}_A , {\bf k}_B - {\bf e}_{\ell}] \frac{{k_{B_{\ell}} \choose 2}} {{k_{A_{\ell}} + k_{B_{\ell}} \choose 2}} \Bigg] Pr (C_{\ell}) \\ & = \frac {1} {\sum_{{\ell} = 1}^L {k_{A_{\ell}} + k_{B_{\ell}} \choose 2}} + \sum_{{\ell} = 1}^L \Bigg\{ \Bigg[ E_{\tau_{AB}} [V_{AB} \mid {\bf k}_A - {\bf e}_{\ell} , {\bf k}_B] \frac{{k_{A_{\ell}} \choose 2}} {{k_{A_{\ell}} + k_{B_{\ell}} \choose 2}} \Bigg. \Bigg. \\ &\Bigg. \Bigg. \quad + E_{\tau_{AB}} [V_{AB} \mid {\bf k}_A , {\bf k}_B - {\bf e}_{\ell}] \frac{{k_{B_{\ell}} \choose 2}} {{k_{A_{\ell}} + k_{B_{\ell}} \choose 2}} \Bigg] \frac{{k_{A_{\ell}} + k_{B_{\ell}} \choose 2}} {\sum_{\lambda = 1}^L {k_{A_{\lambda}} + k_{B_{\lambda}} \choose 2}} \Bigg\} \\ & = \frac {1} {\sum_{{\ell} = 1}^L {k_{A_{\ell}} + k_{B_{\ell}} \choose 2}} \Bigg[ 1 + \sum_{{\ell} = 1}^L \Big[ E_{\tau_{AB}} [V_{AB} \mid {\bf k}_A - {\bf e}_{\ell} , {\bf k}_B] {k_{A_{\ell}} \choose 2} \Big. \Bigg. \\ &\Bigg. \Big. \quad + E_{\tau_{AB}} [V_{AB} \mid {\bf k}_A , {\bf k}_B - {\bf e}_{\ell}] {k_{B_{\ell}} \choose 2} \Big] \Bigg]. \end{align*}\end{document}

Above, λ is a “dummy” summation variable. The second equality can be understood as follows. Because the time to the first coalescent event at locus ℓ is exponentially distributed with mean \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$1 / {k_{A_{\ell}} + k_{B_{\ell}} \choose 2} (\ell = 1 , \ldots , L)$$\end{document} , the time to the first coalescence at some locus is distributed as the minimum of L such random variables. Therefore, the expected time to the first coalescent event is \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$1 / \sum_{\lambda = 1}^L {k_{A_{\lambda}} + k_{B_{\lambda}} \choose 2}$$\end{document} coalescent units. We must always wait this long on average before the first interspecific coalescent event. Given that the first coalescence occurs at locus ℓ, if the coalescence occurs among lineages from taxon A, an event that occurs with probability \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${k_{A_{\ell}} \choose 2} / {k_{A_{\ell}} + k_{B_{\ell}} \choose 2}$$\end{document} , we must wait on average an additional \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{{\tau}_{AB}} \left[V_{AB} \mid {\bf k}_A - {\bf e}_{\ell} , {\bf k}_B \right]$$\end{document} time units. Similarly, with probability \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${k_{B_{\ell}} \choose 2} / {k_{A_{\ell}} + k_{B_{\ell}} \choose 2}$$\end{document} , we must wait an additional \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} \left[V_{AB} \mid {\bf k}_{A} , {\bf k}_{B} - {\bf e}_{\ell} \right]$$\end{document} time units on average. Finally, if the first coalescence at locus ℓ is interspecific, an event that has probability \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\frac {k_{A_{\ell}} k_{B_{\ell}}} {{k_{A_{\ell}} + k_{B_{\ell}} \choose 2}}$$\end{document} , no further waiting is necessary.

In the third equality, the term \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$1 / \sum_{\lambda = 1}^L {k_{A_{\lambda}} + k_{B_{\lambda}} \choose 2}$$\end{document} does not depend on ℓ and can be brought outside. Additionally, because the time to the first coalescence at locus ℓ is exponentially distributed with mean \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$1 / {k_{A_{\ell}} + k_{B_{\ell}} \choose 2}$$\end{document} , the first coalescence occurs at locus ℓ with probability \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$Pr (C_\ell) = {k_{A_{\ell}} + k_{B_{\ell}} \choose 2} / \sum_{\ell = 1}^L {k_{A_{\ell}} + k_{B_{\ell}} \choose 2}$$\end{document} .

12. Appendix B

Derivation of Equation (7)

In Appendix B, we rely on results from Rosenberg (2003) to derive the probability distribution of M, the number of coalescent events up to and including the first interspecific coalescence, counting backwards in time from the divergence time.

Suppose that k_A and k_B lineages from taxa A and B, respectively, remain at time τ_AB. Equation (A8) of Rosenberg (2003) gives the probability \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$F_w^{A , B} (k_{A} , k_{B})$$\end{document} that an interspecific coalescence occurs among these lineages on or before the (k − w)th coalescence, where k = k_A + k_B. This probability is \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*}F_w^{A , B} (k_A , k_B) = \sum \limits_{x = \max (1 , w + 1 - k_B)}^{k_A} \sum \limits_{y = \max (1 , w + 1 - x)}^{k_B} \frac {I_{k_A , x} I_{k_B , y} {I_{k_A + {k_B , w}}}} W (k_A - x , k_B - y) xyI_{x + y - 1 , w} , \tag {{\rm B} .1} \end{align*}\end{document}

where I_n,k = [n!(n − 1)!]/[2^n−kk!(k − 1)!] (Rosenberg, 2003) is the number of ways in which n lineages can coalesce down to k lineages, and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$W (i , j) = {i + j \choose i}$$\end{document} (Rosenberg, 2003) is the number of ways of “interweaving” the coalescent events among lineages only from taxon A with the coalescent events among lineages only from taxon B.

Each term in the summation (B.1) is the joint probability that the first interspecific coalescence occurs when the k_A and k_B lineages have x and y ancestors, respectively, and that the first interspecific coalescence occurs on or before the (k − w)th coalescence. If w = 1, then each term is just the probability that the first interspecific coalescence occurs when the k_A and k_B lineages have x and y ancestors.

Since different choices of x and y (say, (x₁, y₁) and (x₂, y₂) where x₁ ≠ x₂ or y₁ ≠ y₂, or both) correspond to mutually exclusive events, and since the sum x + y specifies M through the relationship x + y = k − M + 1, to derive the probability that the first interspecific coalescence is the Mth coalescence (Equation 7), we can set w = 1 and sum over all x and y such that x + y = k − M + 1, i.e., over all mutually exclusive events corresponding to the case in which the first interspecific coalescence is the Mth coalescence.

To determine the values of x and y corresponding to the case M = m, we can write x = k − m + 1 − y. Note that x is at most k_A and at least 1, and thus, 1 ≤ x ≤ min{k − m, k_A}. Similarly, by symmetry in x and y, 1 ≤ y ≤ min{k − m, k_B}, giving x = k − m + 1 − y ≥ k − m + 1 − min{k − m, k_B} = max{1, k_A − m + 1}. This inequality yields the constraint max{1, k_A − m + 1} ≤ x ≤ min{k − m, k_A}. Thus, we obtain \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} Pr (& M = m) \\ & = \sum_{x = \max \{1 , k_A - m + 1\}}^{\min \{k - m , k_A\}} {I_{k_A , x} \frac{I_{k_B , k - m + 1 - x}} {I_{k_A + k_B , 1}}} W (k_A - x , k_B - y) x (k - m + 1 - x) I_{k_A + k_B - m , 1.} \end{align*}\end{document}

Making the change of variables η = k_A − x and noting that k_B − y = m − 1 − η because k_A + k_B − m + 1 = x + y, we get \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} Pr (M & = m) \\ & = \frac{I_ {k_A + k_B - m , 1} {{I_{k_A + k_B , 1}}}} \sum_{\eta = \max \{0 , m - k_B \}}^{\min \{m - 1 , k_A - 1 \}} W (\eta , m - 1 - \eta) I_{k_{A,}} \,_{k_{A^{- \eta}}} I_{k_{B,}} \,_{k_{B^{- m + 1 + \eta}}} \\ & \quad \times (k_A - \eta) (k_B - m + 1 + \eta). \end{align*}\end{document}

Using I _n,k  = [n!(n − 1)!]/[2^n−kk!(k − 1)!] and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$W (i , j) = {i + j \choose i}$$\end{document} , we get \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} Pr (M = m) & \\ & = \frac{I_ {k_A + k_B - m , 1} {{I_{k_A + k_B , 1}}}} \sum_{\eta = \max \{0 , m - k_B \}}^{\min \{m - 1 , k_A - 1 \}} {m - 1 \choose \eta} \frac {k_{A} ! (k_{A} - 1) !} {2^{k_A - (k_A - \eta)} (k_{A} - \eta) ! (k_{A} - \eta - 1) !} \\ & \times \frac {k_{B} ! (k_{B} - 1) !} {2^{k_B - (k_B - m + 1 + \eta)} (k_{B} - m + 1 + \eta) ! (k_{B} - m + \eta) !} \\ & \times (k_{A} - \eta) (k_{B} - m + 1 + \eta) \\ & = \frac {I_{k_A + k_B - m , 1}} {I_{k_A + k_B , 1}} \sum_{\eta = \max \{0 , m - k_B \}}^{\min \{m - 1 , k_A - 1 \}} {m - 1 \choose \eta} \frac {k_{A} !^2} {k_{A} (k_{A} - \eta - 1) !^2} \frac {k_{B} !^2} {k_{B} (k_{B} - m + \eta) !^2} \frac {1} {2^{m - 1}} \\ & = \frac {I_{k_{A} + k_{B} - m , 1}} {2^{m - 1} k_A k_B I_{k_A + k_B , 1}} \sum_{\eta = \max \{0 , m - k_B \}}^{\min \{m - 1 , k_A - 1 \}} {m - 1 \choose \eta} k_{A_{[\eta + 1]}}^2 k_{B_{[m - \eta]}}^2 , \qquad \qquad \qquad ({\rm B} .2) \end{align*}\end{document}

where k_[i] = k!/(k − i)!.

When either k_A = 1 or k_B = 1, Equation (B.2) has a particularly simple form. Without loss of generality, suppose that k_B = 1. Then max{0, m − k_B} = m − 1 because m ≥ 1, and min{m − 1, k_A − 1} = m − 1 because m ≤ k_A + k_B − 1 = k_A. Therefore, using k = k_A + 1, Equation (B.2) simplifies as follows: \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} \begin{split} Pr (M = m) & = \frac {1} {k_{A} 2^{m - 1}} \frac {2^{k_A + 1 - 1}} {(k_{A} + 1) !k_{A} !} \frac {(k_{A} + 1 - m) ! (k_{A} + 1 - m - 1) !} {2^{k_{A} + 1 - m - 1}} k^{2}_{A {[m]}} 1^2_{[1]} \\ & = \frac {2 (k_{A} + 1 - m)} {k_A (k_{A} + 1)} \frac {(k_{A} - m) ! (k_{A} - m) !} {k_{A} !k_{A} !} \frac {k_{A} !^2} {(k_{A} - m) !^2} \\ & = \frac {2 (k - m)} {k (k - 1)}. \end{split} \tag {{\rm B} .3} \end{align*}\end{document}

13. Appendix C.

Computational complexity of approximate iGLASS

We now compute the computational complexity of the approximate iGLASS method, Equation (21). To compute the iGLASS correction for each pair of taxa X and Y in \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal S}$$\end{document} , we first evaluate Equation (20) for many different values of τ_XY. In particular, to numerically obtain the inverse in Equation (21), we compute Equation (20) for each divergence time estimate τ_XY in the set \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{i \Delta t \}_{i = 0}^{\Gamma_{XY}}$$\end{document} , where Δt is a fixed time-step and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\Gamma_{XY} = \lceil \hat{t}_{XY} / \Delta t \rceil$$\end{document} . We then estimate τ_XY by the value \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tau \in \{i \Delta t \}_{i = 0}^{\Gamma_{XY}}$$\end{document} that minimizes the quantity \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\mid \tau + \tilde{E}_{\tau} [V_{XY}] - \hat{t}_{XY} \mid$$\end{document} .

To evaluate the integral in Equation (20), we assume that numerical integration is carried out by computing the Riemann sum with fixed step-size Δt. We truncate the outer integral at PΔt, where P is large enough that the tail of the outer integral in Equation (20) is smaller than some predefined value ε > 0. For a given value of ε, a sufficiently-large value of P can be found by bounding the integral in Equation (20). The bound can be obtained by noting that \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_n^z \geq 1$$\end{document} for all n and z in \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\mathbb R}^+$$\end{document} , and thus, the integrand in Equation (20) is smaller than exp{ − Lt}, which is easily integrated. Converting the integrals in Equation (20) to summations gives \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} \tilde{E}_{\tau_{XY}} [V_{XY}] \approx \sum_{\alpha = 0}^{P} \Delta t \exp \left\{- \sum_{{\ell} = 1}^L \int_{z = 0}^{\alpha \Delta t} \varphi_{E_{\tau_{XY}} [K_{X_{\ell}}]}^z \varphi_{E_{\tau_{XY}} [K_{Y_{\ell}}]}^z dz \right\} \\ \approx \sum_{\alpha = 0}^{P} \Delta t \exp \left\{- \sum_{{\ell} = 1}^L \sum_{\beta = 0}^{\alpha} \Delta t \ \varphi_{E_{\tau_{XY}} [K_{X_{\ell}}]}^{\beta \Delta t} \varphi_{E_{\tau_{XY}} [K_{Y_{\ell}}]}^{\beta \Delta t} \right\} . & {(\rm C}.1) \end{align*}\end{document}

Once \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{E_{\tau_{XY}} [K_{X_{\ell}}]}^{\beta \Delta t}$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{E_{\tau_{XY}} [K_{Y_{\ell}}]}^{\beta \Delta t}$$\end{document} have been pre-computed and stored for all values at which they are evaluated in the summation, the exponent in Equation (C.1) requires O(Lα) operations, where α is the index in the outermost summation. Thus, we have the following result:

After pre-computing the terms in the summand, the summation (C.1) requires O(LP²) operations.

For each taxon \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$X \in {\cal S}$$\end{document} , let \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\Gamma_X = \max_{Y \in {\cal S} , Y \neq X} \Gamma_{XY}$$\end{document} ; in other words, Γ _X Δt is, to precision Δt, the maximum pairwise divergence time between taxon X and any other taxon. For each \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\ell = 1 , \ldots , L$$\end{document} and for each \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$X \in {\cal S}$$\end{document} , we must ultimately compute \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau} [K_{X_{\ell}}]$$\end{document} for each \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tau \in \{i \Delta t \}_{i = 0}^{\Gamma_X}$$\end{document} , and we must compute \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{E_{\tau} [K_{X_{\ell}}]}^{t}$$\end{document} for each \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$t \in \{i \Delta t \}_{i = 0}^{P}$$\end{document} and for each \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tau \in \{j \Delta t \}_{j = 0}^{\Gamma_X}$$\end{document} . However, note that \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_n^{t + \Delta t} = \varphi_{\varphi_n^t}^{\Delta t}$$\end{document} , and note that \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi^\tau_{n_{X_\ell}} = E_{\tau} [K_{X_\ell}]$$\end{document} by definition (Equations (15) and (16)). Therefore, we have \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{E_{\tau} [K_{X_{\ell}}]}^{i \Delta t} = \varphi_{\varphi^\tau_{n_{X_\ell}}}^{i \Delta t} = \varphi_{n_{X_{\ell}}}^{\tau + i \Delta t}$$\end{document} for all \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$i \in {\rm Z}^+$$\end{document} , and thus, it suffices to pre-compute \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{n_{X_{\ell}}}^\tau$$\end{document} for all \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tau \in \{j \Delta t \}_{j = 0}^{\Gamma_X + P}$$\end{document} for each \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$X \in {\cal S}$$\end{document} and for each \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\ell \ (\ell = 1 , \ldots , L)$$\end{document} .

Let \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$n = \max_{X \in {\cal S}} \max_{\ell \in \{1 , \ldots , L \}} n_{X_{\ell}}$$\end{document} be the maximal number of lineages sampled from any taxon, and let \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$Q = P + \max_{X \in {\cal S}} \Gamma_{X}$$\end{document} . Then for a given taxon \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$X \in {\cal S}$$\end{document} and for a given \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\ell \ (\ell = 1 , \ldots , L)$$\end{document} , the amount of time needed to compute \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{n_{X_{\ell}}}^\tau$$\end{document} for all \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tau \in \{j \Delta t \}_{j = 0}^{\Gamma_X + P}$$\end{document} is bounded by the time needed to compute \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_n^{t}$$\end{document} for all \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$t \in \{i \Delta t \}_{i = 0}^Q$$\end{document} .

Because the summand in Equation (15) requires O(n) operations, (a rising factorial and a falling factorial totaling n multiplications), computing Equation (15) for a given value of t requires O(n²) operations. Therefore, evaluating (15) for each time t in \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\{i \Delta t \}_{i = 0}^Q$$\end{document} requires O(n²Q) operations, and pre-computing \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_n^{t}$$\end{document} for all \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$t \in \{i \Delta t \}_{i = 0}^Q$$\end{document} also requires O(n²Q) operations. This gives the following result:

Pre-computing \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{n_{X_{\ell}}}^\tau$$\end{document} for all \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tau \in \{j \Delta t \}_{j = 0}^{\Gamma_X + P}$$\end{document} for all \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$X \in {\cal S}$$\end{document} and for all ℓ \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$(\ell = 1 , \ldots , L)$$\end{document} requires \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$O (n^2 Q L \mid {\cal S} \mid)$$\end{document} operations.

Once all values of \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{n_{X_{\ell}}}^\tau$$\end{document} have been pre-computed and stored, Equation (C.1) must be computed for each \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tau \in \{0 , 1 , \ldots , \Gamma_{XY} \}$$\end{document} for each pair of taxa \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$X , Y \in {\cal S}$$\end{document} . Equation (C.1) requires O(LP²) operations for each value of τ. Therefore, because Γ _XY  ≤ Q, computing (C.1) for a pair of taxa requires O(LP²Q) operations. Because P ≤ Q, this simplifies to O(LQ³) operations. Therefore, computing (C.1) for all \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\mid {\cal S} \mid \choose 2}$$\end{document} pairs of taxa requires \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$O (LQ^3 \mid {\cal S} \mid^{2})$$\end{document} operations. Combining this quantity with the number of operations necessary to pre-compute the values of \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{n_{X_{\ell}}}^\tau$$\end{document} gives the following result:

Including all pre-computations, the total number of operations required to compute Equation (C.1) for all \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\mid {\cal S} \mid \choose 2}$$\end{document} pairs of taxa is \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$O (n^2QL \mid {\cal S} \mid + LQ^3 \mid {\cal S} \mid^2)$$\end{document} .

Note that once all values of \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_{n_{X_{\ell}}}^\tau$$\end{document} have been pre-computed and stored, the cost of computing (C.1) does not depend on the magnitude of the \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$n_{X_\ell}$$\end{document} , only on the number of terms in the summation. Thus, the complexity only depends on n through the pre-computation step.

The only other computations needed to compute the approximate iGLASS correction are those associated with finding arg \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\min_{\tau} \mid \tau + \tilde{E}_{\tau} [V_{AB}] - \hat{t}_{AB} \mid$$\end{document} and those associated with the single-linkage clustering step. We must perform \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\mid {\cal S} \mid \choose 2}$$\end{document} searches to find the value of τ that minimizes \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\mid \tau + \tilde{E}_{\tau} [V_{AB}] - \hat{t}_{AB} \mid$$\end{document} for each of the \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\mid {\cal S} \mid \choose 2}$$\end{document} pairs of taxa. An exhaustive search is bounded by the number of values of τ, which is always less than or equal to Q. Thus, correcting the GLASS method requires \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$O (Q \mid {\cal S} \mid^2)$$\end{document} operations. Finally, single-linkage clustering requires at most \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$O (\mid {\cal S} \mid^2)$$\end{document} operations (Gordon, 1996). Thus, the entire correction procedure requires \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$O (n^2QL \mid {\cal S} \mid + LQ^3 \mid {\cal S} \mid^2 + Q \mid {\cal S} \mid^2 + \mid {\cal S} \mid^2)$$\end{document} operations. Terms can be combined to get the following result:

The entire approximate iGLASS correction procedure requires \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$O (n^2QL \mid {\cal S} \mid + LQ^3 \mid {\cal S} \mid^2)$$\end{document} operations.

It is useful to compare the complexity of approximate iGLASS to the complexity of GLASS for a given precision. The choices of Δt and P determine the precision in computing the approximate iGLASS correction, in other words, the error between the outcome of the numerical steps that we have just outlined, and the outcome of exactly computing Equation (20) and exactly solving Equation (21). Together, Δt and P determine \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$Q = P + \lceil \max_{XY \in {\cal S}} \hat{t}_{XY} / \Delta t \rceil$$\end{document} . Thus, Q is a tuning parameter that affects the precision in our numerical steps. For fixed Q, the complexity of approximate iGLASS is \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$O (n^2 L \mid {\cal S} \mid + L \mid {\cal S} \mid^2)$$\end{document} . In comparison, a similar analysis demonstrates that the GLASS method requires \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$O (n^2L \mid {\cal S} \mid^2 + \mid {\cal S} \mid^3)$$\end{document} operations.

14. Appendix D

Consistency of GLASS for divergence times

Mossel and Roch (2010) proved that the GLASS method is a consistent estimator of the species tree topology as the number of loci approaches infinity. Liu et al. (2010) proved that the GLASS estimator is consistent for pairwise divergence times in the case in which a single lineage is sampled per taxon.

Here, we prove that GLASS is a consistent estimator of pairwise divergence times in the case in which arbitrarily many lineages are sampled per taxon. Our argument is a minor extension of the consistency proof in Liu et al. (2010).

Theorem D.1

Consider two taxa, A and B, with divergence time τ_AB. The GLASS estimator \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\hat t}_{AB}$$\end{document} is a consistent estimator of τ_AB.

Proof

At each locus ℓ \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$(\ell = 1 , \ldots , L)$$\end{document} , consider a lineage a_ℓ sampled at random from taxon A and a lineage b_ℓ sampled at random from taxon B. The time \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$V_{AB}^{\ell}$$\end{document} to the first interspecific coalescence at locus ℓ is less than or equal to the coalescence time between a_ℓ and b_ℓ, which we denote by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$V_{a_\ell , b_\ell}$$\end{document} . Therefore, using the fact that the GLASS estimate is given by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\hat t}_{AB} = \tau_{AB} + V_{AB}$$\end{document} , and following Liu et al. (2010), we obtain \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$Pr (\mid \hat t_{AB} - \tau_{AB} \mid > \epsilon) = Pr (V_{AB} > \epsilon) = Pr (\min_\ell V_{AB}^\ell > \epsilon) \leq Pr (\min_\ell V_{a^\ell , b^\ell} > \epsilon) = \prod_{\ell = 1}^L Pr (V_{a^\ell , b^\ell} > \epsilon) = e^{- L \epsilon}$$\end{document} . Here, to obtain the last equality, we have used the fact that \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$V_{a^\ell , b^\ell}$$\end{document} is exponentially distributed with mean 1 coalescent unit of N generations. Thus, we have \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} 0 \leq Pr (\mid {\hat t}_{AB} - \tau_{AB} \mid > \epsilon) \leq e^{- L \epsilon} , \end{align*}\end{document}

from which it follows that \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$Pr (\mid \hat{t}_{AB} - \tau_{AB} \mid > \ \epsilon) \rightarrow 0$$\end{document} as L → ∞ by the “squeeze theorem.”   ▪

15. Appendix E

iGLASS and approximate iGLASS are consistent estimators of pairwise divergence times

Here, we prove that the expectation \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB}]$$\end{document} of the difference \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$V_{AB} = \hat{t}_{AB} - \tau_{AB}$$\end{document} between the GLASS estimate \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{t}_{AB}$$\end{document} and the divergence time τ_AB is bounded above by 1/L. Thus, \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB}] \rightarrow 0$$\end{document} as L → ∞. Using Equation (1), we then show that the difference between the GLASS estimator and the iGLASS estimator is bounded above by 1/L. Thus, the difference goes to 0 as L → ∞. A similar result is proven for the expectation \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\tilde E}_{\tau_{AB}} [V_{AB}]$$\end{document} used in the approximate iGLASS correction (Equation 21).

Since GLASS is a consistent estimator of pairwise divergence times, these results can be used to show that exact iGLASS and approximate iGLASS are consistent estimators of pairwise divergence times, as they converge to the same limit as the GLASS estimator in the limit L → ∞.

Lemma E.1

For taxa A and B, let \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB}]$$\end{document} be the expectation of the difference \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$V_{AB} = \hat{t}_{AB} - \tau_{AB}$$\end{document} between the GLASS estimate \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{t}_{AB}$$\end{document} and the divergence time τ_AB. Then \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_{\tau_{AB}} [V_{AB}] \leq 1 / L$$\end{document} .

Proof

In Theorem D.1, we saw that \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$0 \leq Pr (\mid \hat{t}_{AB} - \tau_{AB} \mid > \epsilon) \leq e^{- L \epsilon}$$\end{document} for all \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\epsilon > 0$$\end{document} . Thus, \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} E_{\tau_{AB}} [V_{AB}] & = \int_{v = 0}^\infty Pr (V_{AB} > v) dv \\ & = \int_{v = 0}^{\infty} Pr (\mid \hat {t}_{AB} - \tau_{AB} \mid > v) dv \\ & \leq \int_{v = 0}^{\infty} e^{- Lv} dv \\ & = \frac {1} {L} , \end{align*}\end{document}

proving the result.   ▪

Lemma E.2

The approximation \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde{E}_{\tau_{AB}} [V_{AB}]$$\end{document} satisfies \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde{E}_{\tau_{AB}} [V_{AB}] \leq 1 / L$$\end{document} .

Proof

For any n and t, the expected number of lineages \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\varphi_n^t$$\end{document} remaining at any given time t is at least 1. Therefore, for any ℓ \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$(\ell = 1 , \ldots , L)$$\end{document} \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*}\int_{z = 0}^{t} \varphi_{E_{\tau_{AB}} [K_{A_{\ell}}]}^z \varphi_{E_{\tau_{AB}} [K_{B_{\ell}}]}^z dz \ge t.\end{align*}\end{document}

Consequently, \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} \begin{align*} \tilde {E}_{\tau_{AB}} [V_{AB}] & = \int_{t = 0}^{\infty} \exp \left\{- \sum_{\ell = 1}^L \int_{z = 0}^{t} \varphi_{E_{\tau_{AB}} [K_{A_{\ell}}]}^z \varphi_{E_{\tau_{AB}} [K_{B_{\ell}}]}^z dz \right\} dt \\ & \leq \int_{t = 0}^{\infty} \exp \left\{- \sum_{\ell = 1}^L t \right\} dt \\ & = \int_{t = 0}^{\infty} e^{- L t} dt \\ & = \frac {1} {L} , \end{align*}\end{document}

proving the result.   ▪

The following corollary proves that after the correction procedure (Equation 2), both the exact and approximate iGLASS estimates differ from the GLASS estimate by at most 1/L coalescent units.

Corollary E.3

For two taxa A and B, let \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$C_{AB} = \hat{t}_{AB} - \hat{\tau}_{AB}$$\end{document} and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde{C}_{AB} = \hat{t}_{AB} - \tilde{\tau}_{AB}$$\end{document} be the differences between the GLASS estimate and the exact and approximate iGLASS estimates, respectively. Then C_AB ≤1/L and \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde{C}_{AB} \leq 1 / L$$\end{document} .

Proof

Using Equation (2), if \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E_0 [V_{AB}] \leq \hat{t}_{AB}$$\end{document} , then the iGLASS estimate is obtained by solving \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\tau}_{AB} = \hat{t}_{AB} - E_{\tau_{AB}} [V_{AB}]$$\end{document} for τ_AB. In this case, the difference C_AB is at most 1/L by Lemma E.1. On the other hand, if \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{t}_{AB} \in [0 , E_0 [V_{AB}])$$\end{document} , then the iGLASS estimate is given by \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\hat{\tau}_{AB} = 0$$\end{document} . Since [0, E₀[V _AB ]) ⊆ [0, 1/L) by Lemma E.1, we have \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\mid \hat{t}_{AB} - \hat{\tau}_{AB} \mid < 1 / L$$\end{document} . Thus, in both cases, C_AB ≤ 1/L.

The same argument using Lemma E.2 and Equation (21) rather than Lemma E.1 and Equation (2) establishes \documentclass{aastex} \usepackage{amsbsy} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{bm} \usepackage{mathrsfs} \usepackage{pifont} \usepackage{stmaryrd} \usepackage{textcomp} \usepackage{portland, xspace} \usepackage{amsmath, amsxtra} \pagestyle{empty} \DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde{C}_{AB} \leq 1 / L$$\end{document} , proving the result.   ▪