Inferring phenotypic trait evolution on large trees with many incomplete measurements

TL;DR

This paper introduces a scalable inference method for analyzing phenotypic trait evolution across large phylogenetic trees with many missing measurements, improving computational efficiency significantly.

Contribution

It presents an analytical integration technique for missing data under a multivariate Brownian diffusion model that scales linearly with the number of taxa.

Findings

Achieves up to two orders-of-magnitude speedup over existing methods.

Successfully applied to diverse biological datasets including mammals, prokaryotes, and HIV.

Extends the model to incorporate sampling error and residual variance.

Abstract

Comparative biologists are often interested in inferring covariation between multiple biological traits sampled across numerous related taxa. To properly study these relationships, we must control for the shared evolutionary history of the taxa to avoid spurious inference. Existing control techniques almost universally scale poorly as the number of taxa increases. An additional challenge arises as obtaining a full suite of measurements becomes increasingly difficult with increasing taxa. This typically necessitates data imputation or integration that further exacerbates scalability. We propose an inference technique that integrates out missing measurements analytically and scales linearly with the number of taxa by using a post-order traversal algorithm under a multivariate Brownian diffusion (MBD) model to characterize trait evolution. We further exploit this technique to extend the MBD model to account for sampling error or non-heritable residual variance. We test these methods to examine mammalian life history traits, prokaryotic genomic and phenotypic traits, and HIV infection traits. We find computational efficiency increases that top two orders-of-magnitude over current best practices. While we focus on the utility of this algorithm in phylogenetic comparative methods, our approach generalizes to solve long-standing challenges in computing the likelihood for matrix-normal and multivariate normal distributions with missing data at scale.

Figures

Tables

Table 1. Table 1: Algorithmic improvement. We report MCMC sampling efficiency through effective sample size (ESS) that shows both a decrease in autocorrelation (as shows by ESS / Sample) and in the overall work required per sample (as shown by Samples / Hour).

Data set	Model	Integration method	ESS/hour			ESS/sample		Samples /hour
			minimum	median		minimum	median
Mammals	Diffusion only	Analytic	1,200	3,600		0.044	0.13	27,000
		Sampling	3	9.8		0.0043	0.014	700
		Speed-up	410$\times$	360$\times$		10$\times$	9.3$\times$	39$\times$
	Diffusion with residual	Analytic	140	320		0.0062	0.015	22,000
		Sampling	0.38	3		0.00024	0.0018	1,600
		Speed-up	350$\times$	110$\times$		26$\times$	7.9$\times$	14$\times$
HIV	Diffusion only	Analytic	100,000	220,000		0.31	0.66	320,000
		Sampling	1,500	8,500		0.01	0.057	150,000
		Speed-up	65$\times$	25$\times$		30$\times$	12$\times$	2.2$\times$
	Diffusion with residual	Analytic	1,600	2,500		0.0061	0.0096	260,000
		Sampling	5.1	8.7		9.9e-5	0.00017	51,000
		Speed-up	320$\times$	290$\times$		61$\times$	56$\times$	5.2$\times$

Table 2. Table 2: Missing data summary for all three examples.

Data set	Trait	Number observed	Percent missing
Mammals $N=3690$	Body mass	3508	4.9%
	Litter size	2538	31.2%
	Gestation length	1427	61.3%
	Weaning age	1253	66.0%
	Litter frequency	1231	66.6%
	Neonatal body mass	1108	70.0%
	Age at first birth	945	74.4%
	Reproductive lifespan	748	79.7%
	Total	12758	56.8%
Prokaryotes $N=705$	Cell diameter	690	2.1%
	Cell length	657	6.8 %
	Genome length	563	20.1%
	GC content	563	20.1%
	Coding sequence length	558	20.9%
	Optimal temperature	548	22.2%
	Optimal pH	487	30.9%
	Total	4066	17.6%
HIV $N=1536$	GSVL	1536	0.0%
	SPVL	1536	0.0%
	CD4 slope	1102	28.3%
	Total	4174	9.4%

Equations

$${\mathbf{X}}_k\sim \text{MVN}\left({\mathbf{X}}_{\text{pa}(k)},t_k\mathbf{\Sigma }\right)\text{ for }k=1,\dots ,2N-2,$$

$$p({\mathbf{Y}}^{\text{obs}}|\mathbf{\Sigma },\mathcal{F},{\mathbf{\mu }}_0,{\kappa }_0)=\int p({\mathbf{Y}}^{\text{obs}},{\mathbf{Y}}^{\text{mis}}|\mathbf{\Sigma },\mathcal{F},{\mathbf{\mu }}_0,{\kappa }_0)\text{d}{\mathbf{Y}}^{\text{mis}},$$

$$\text{mbox}vec\left[\mathbf{Y}|\mathbf{\Sigma },\mathcal{F},{\mathbf{\mu }}_0,{\kappa }_0\right]\sim \text{MVN}\left(\text{mbox}vec\left[1_N{\mathbf{\mu }}_0^t\right],\mathbf{\Sigma }\otimes \mathbf{{\rm Y}}\right),$$

$${\mathbf{P}}_i={\mathbf{C}}_i\text{diag}\infty {\mathbf{C}}_i^t,$$

$${\mathbf{P}}_i^{-}={\mathbf{C}}_i\text{diag}\infty {\mathbf{C}}_i^t.$$

$$\log \hat{\varphi }\left(\mathbf{z};\mathbf{\mu },\mathbf{P}\right)=\frac{1}{2}\log \widehat{\text{det}}\left(\mathbf{P}\right)-\frac{\text{rank}\left(\mathbf{P}\right)}{2}\log 2\pi -\frac{1}{2}{\left(\mathbf{z}-\mathbf{\mu }\right)}^t\mathbf{P}\left(\mathbf{z}-\mathbf{\mu }\right),$$

$$p({\mathbf{Y}}_{\lfloor k\rfloor }^{\text{obs}}|{\mathbf{X}}_k,\mathbf{\Sigma },\mathcal{F})=r_k\hat{\varphi }\left({\mathbf{X}}_k;{\mathbf{m}}_k,{\mathbf{P}}_k\right),$$

$$p({\mathbf{Y}}_{\lfloor \ell \rfloor }^{\text{obs}}|{\mathbf{X}}_{\ell },\mathbf{\Sigma },\mathcal{F})=p({\mathbf{Y}}_{\lfloor j\rfloor }^{\text{obs}}|{\mathbf{X}}_{\ell },\mathbf{\Sigma },\mathcal{F})p({\mathbf{Y}}_{\lfloor k\rfloor }^{\text{obs}}|{\mathbf{X}}_{\ell },\mathbf{\Sigma },\mathcal{F})$$

$$\begin{array}{rl}p({\mathbf{Y}}_{\lfloor j\rfloor }^{\text{obs}}|{\mathbf{X}}_{\ell },\mathbf{\Sigma },\mathcal{F})& =\int p({\mathbf{Y}}_{\lfloor j\rfloor }^{\text{obs}}|{\mathbf{X}}_j,\mathbf{\Sigma },\mathcal{F})p({\mathbf{X}}_j|{\mathbf{X}}_{\ell },\mathbf{\Sigma },\mathcal{F})\text{d}{\mathbf{X}}_j=r_j\hat{\varphi }\left({\mathbf{X}}_{\ell };{\mathbf{m}}_j,{\mathbf{P}}_j^{\star }\right),\end{array}$$

$$\begin{array}{rl}\log p({\mathbf{Y}}_{\lfloor \ell \rfloor }^{\text{obs}}|{\mathbf{X}}_{\ell },\mathbf{\Sigma },\mathcal{F})& =\log r_j+\log r_k+\log \hat{\varphi }\left({\mathbf{X}}_{\ell };{\mathbf{m}}_j,{\mathbf{P}}_j^{\star }\right)+\log \hat{\varphi }\left({\mathbf{X}}_{\ell };{\mathbf{m}}_k,{\mathbf{P}}_k^{\star }\right)\\ & =\log r_{\ell }+\log \hat{\varphi }\left({\mathbf{X}}_{\ell };{\mathbf{m}}_{\ell },{\mathbf{P}}_{\ell }\right),\end{array}$$

$$\begin{array}{r}\log r_{\ell }=\log r_j+\log r_k+\frac{1}{2}\log \widehat{\text{det}}\left({\mathbf{P}}_j^{\star }\right)+\frac{1}{2}\log \widehat{\text{det}}\left({\mathbf{P}}_k^{\star }\right)-\frac{{\Delta }_{jk\ell }}{2}\log 2\pi \\ -\frac{1}{2}\log \widehat{\text{det}}\left({\mathbf{P}}_{\ell }\right)-\frac{1}{2}\left({\mathbf{m}}_j^t{\mathbf{P}}_j^{\star }{\mathbf{m}}_j+{\mathbf{m}}_k^t{\mathbf{P}}_k^{\star }{\mathbf{m}}_k-{\mathbf{m}}_{\ell }^t{\mathbf{P}}_{\ell }{\mathbf{m}}_{\ell }\right).\end{array}$$

$$\begin{array}{rl}p({\mathbf{Y}}^{\text{obs}}|\mathbf{\Sigma },\mathcal{F},{\mathbf{\mu }}_0,{\kappa }_0)& =\int p({\mathbf{Y}}^{\text{obs}}|{\mathbf{X}}_{2N-1},\mathbf{\Sigma },\mathcal{F})p({\mathbf{X}}_{2N-1}|\mathbf{\Sigma },{\mathbf{\mu }}_0,{\kappa }_0)\text{d}{\mathbf{X}}_{2N-1}\\ & =\int r_{2N-1}\hat{\varphi }\left({\mathbf{X}}_{2N-1};{\mathbf{m}}_{2N-1},{\mathbf{P}}_{2N-1}\right)\hat{\varphi }\left({\mathbf{X}}_{2N-1};{\mathbf{\mu }}_0,{\kappa }_0{\mathbf{\Sigma }}^{-}\right)\text{d}{\mathbf{X}}_{2N-1}\\ & =r_{\text{full}}\int \hat{\varphi }\left({\mathbf{X}}_{2N-1};{\mathbf{m}}_{\text{full}},{\mathbf{P}}_{\text{full}}\right)\text{d}{\mathbf{X}}_{2N-1},\end{array}$$

$$\begin{array}{rl}\log p({\mathbf{Y}}^{\text{obs}}|\mathbf{\Sigma },\mathcal{F},{\mathbf{\mu }}_0,{\kappa }_0)& =\log r_{\text{full}}\\ & =\log r_{2N-1}-\frac{\text{rank}\left({\mathbf{P}}_{2N-1}\right)}{2}\log 2\pi \\ & +\frac{1}{2}\log \widehat{\text{det}}\left({\mathbf{P}}_{2N-1}\right)+\frac{1}{2}\log \widehat{\text{det}}\left({\kappa }_0{\mathbf{\Sigma }}^{-1}\right)-\frac{1}{2}\log \widehat{\text{det}}\left({\mathbf{P}}_{\text{full}}\right)\\ & -\frac{1}{2}\left({\mathbf{m}}_{2N-1}^t{\mathbf{P}}_{2N-1}{\mathbf{m}}_{2N-1}+{\kappa }_0{\mathbf{\mu }}_0^t{\mathbf{\Sigma }}^{-1}{\mathbf{\mu }}_0-{\mathbf{m}}_{\text{full}}^t{\mathbf{P}}_{\text{full}}{\mathbf{m}}_{\text{full}}\right).\end{array}$$

$$p\left(\mathbf{\Sigma },\mathcal{F},\mathbf{\theta }|{\mathbf{Y}}^{\text{obs}},\mathbf{S}\right)\propto p\left({\mathbf{Y}}^{\text{obs}}|\mathbf{\Sigma },\mathcal{F},\mathbf{\theta }\right)p(\mathcal{F},\mathbf{S})p(\mathbf{\Sigma })p(\mathbf{\theta }),$$

$$\begin{gathered} {\mathbf{\Sigma }}^{-1}\mid \mathbf{X},\mathcal{F},{\mathbf{\mu }}_0,{\kappa }_0,\nu ,{\mathbf{\Lambda }}_0\sim \\ {\text{Wishart}}_P\left[{\mathbf{\Lambda }}_0+{\left(\mathbf{X}-{\mathbf{J}}_N{\mathbf{\mu }}_0^t\right)}^t{\left(\mathbf{{\rm Y}}+\frac{1}{{\kappa }_0}{\mathbf{J}}_N\right)}^{-1}\left(\mathbf{X}-{\mathbf{J}}_N{\mathbf{\mu }}_0^t\right),\nu +N\right]. \end{gathered}$$

$$\begin{array}{rl}p({\mathbf{X}}_{2N-1}|{\mathbf{Y}}^{\text{obs}},{\mathbf{\mu }}_0,{\kappa }_0)& \propto p({\mathbf{Y}}^{\text{obs}}|{\mathbf{X}}_{2N-1})p({\mathbf{X}}_{2N-1}|{\mathbf{\mu }}_0,{\kappa }_0)\\ & \propto \hat{\varphi }\left({\mathbf{X}}_{2N-1};{\mathbf{m}}_{\text{full}},{\mathbf{P}}_{\text{full}}\right)\text{, which implies that}\\ {\mathbf{X}}_{2N-1}|{\mathbf{Y}}^{\text{obs}},{\mathbf{\mu }}_0,{\kappa }_0& \sim \text{MVN}\left({\mathbf{m}}_{\text{full}},{\mathbf{P}}_{\text{full}}\right).\end{array}$$

$$\begin{array}{rl}p({\mathbf{X}}_j|{\mathbf{Y}}_{\lfloor j\rfloor }^{\text{obs}},{\mathbf{X}}_{\text{pa}(j)})& \propto p({\mathbf{Y}}_{\lfloor j\rfloor }^{\text{obs}}|{\mathbf{X}}_j)p({\mathbf{X}}_j|{\mathbf{X}}_{\text{pa}(j)})\\ & \propto \hat{\varphi }\left({\mathbf{X}}_j;{\mathbf{m}}_j,{\mathbf{P}}_j\right)\hat{\varphi }\left({\mathbf{X}}_j;{\mathbf{X}}_{\text{pa}(j)},{\left(t_j\mathbf{\Sigma }\right)}^{-1}\right)\\ & \propto \hat{\varphi }\left({\mathbf{X}}_j;{\mathbf{n}}_j,{\mathbf{Q}}_j\right)\end{array}$$

$${\mathbf{\Sigma }}^{-1}={\mathbf{C}}_i\left(\begin{array}{cc}{\mathbf{S}}_i^{\text{obs}}& {\mathbf{S}}_i^{\text{om}}\\ {\mathbf{S}}_i^{\text{mo}}& {\mathbf{S}}_i^{\text{mis}}\end{array}\right){\mathbf{C}}_i^t,{\mathbf{X}}_i={\mathbf{C}}_i\left(\begin{array}{c}{\mathbf{X}}_i^{\text{obs}}\\ {\mathbf{X}}_i^{\text{mis}}\end{array}\right),\text{ and}{\mathbf{X}}_{\text{pa}(i)}={\mathbf{C}}_i\left(\begin{array}{c}{\mathbf{X}}_{\text{pa}(i)}^{\text{obs}}\\ {\mathbf{X}}_{\text{pa}(i)}^{\text{mis}}\end{array}\right)$$

$$p\left({\mathbf{Y}}_i|{\mathbf{X}}_i\right)=\hat{\varphi }\left({\mathbf{Y}}_i;{\mathbf{X}}_i,\mathbf{\Gamma }\right)\text{ for }i=1,\dots ,N,$$

$$\mathbf{\Gamma }|\mathbf{Y},\mathbf{X}\sim {\text{Wishart}}_P\left({\mathbf{\Lambda }}_s+{\left(\mathbf{Y}-\mathbf{X}\right)}^t\left(\mathbf{Y}-\mathbf{X}\right),{\nu }_s+N\right).$$

$$\mathbf{\Gamma }={\mathbf{C}}_i\left(\begin{array}{cc}{\mathbf{\Gamma }}_i^{\text{obs}}& {{\mathbf{\Gamma }}_i^{\text{mo}}}^t\\ {\mathbf{\Gamma }}_i^{\text{mo}}& {\mathbf{\Gamma }}_i^{\text{mis}}\end{array}\right){\mathbf{C}}_i^t\text{ for }i=1,\dots ,N.$$

$$p\left({\mathbf{Y}}_i^{\text{mis}}|{\mathbf{Y}}_i^{\text{obs}},{\mathbf{X}}_i,\mathbf{\Gamma }\right)=\hat{\varphi }\left({\mathbf{Y}}_i^{\text{mis}};{\mathbf{X}}_i^{\text{mis}}+{{\mathbf{\Gamma }}_i^{\text{mis}}}^{-1}{\mathbf{\Gamma }}_i^{\text{mo}}\left({\mathbf{X}}_i^{\text{obs}}-{\mathbf{Y}}_i^{\text{obs}}\right),{\mathbf{\Gamma }}_i^{\text{mis}}\right).$$

$$\begin{array}{rl}{\mathbf{S}}^2\left(\mathbf{Y}\right)& =\frac{1}{N}\sum _{i=1}^N\left({\mathbf{Y}}_i-\bar{\mathbf{y}}\right){\left({\mathbf{Y}}_i-\bar{\mathbf{y}}\right)}^t=\frac{1}{N}{\left(\mathbf{Y}-\bar{\mathbf{Y}}\right)}^t\left(\mathbf{Y}-\bar{\mathbf{Y}}\right),\end{array}$$

$$\mathbb{E}\left[{\mathbf{S}}^2\left(\mathbf{Y}\right)\right]=\frac{N-1}{N}{\mathbf{\Gamma }}^{-1}+\left(\frac{1}{N}\text{mbox}tr\left[\mathbf{{\rm Y}}\right]-\frac{1}{N^2}{1}_N^t\mathbf{{\rm Y}}{1}_N\right)\mathbf{\Sigma }.$$

$$h_{kl}=\frac{c_{\sigma }{\Sigma }_{kl}}{\sqrt{\left(c_{\sigma }{\Sigma }_{kk}+c_{\gamma }{\Gamma }_{kk}^{-1}\right)\left(c_{\sigma }{\Sigma }_{ll}+c_{\gamma }{\Gamma }_{ll}^{-1}\right)}},$$

$$\mathbf{\Sigma }$$

$${\mathbf{P}}_i$$

$${\mathbf{P}}_i^{\star }$$

$$\displaystyle=\mathbf{C}_{i}\,\left(\left(\text{diag}\left[\raisebox{1.00006pt}{$\infty$}\mathbf{I},\tilde{\mathbf{P}_{i}},0\mathbf{I}\right]\right)^{-}+\text{diag}\left[t_{i}\left(\begin{array}[]{ll}\boldsymbol{\Sigma}^{\text{\tiny{{obs}}}}_{i}&\boldsymbol{\Sigma}^{\text{\tiny{{ol}}}}_{i}\\ -&\boldsymbol{\Sigma}^{\text{\tiny{{lat}}}}_{i}\\ \end{array}\right),0\mathbf{I}\right]\right)^{-}\mathbf{C}_{i}^{t}$$

$$\displaystyle=\mathbf{C}_{i}\,\left(\text{diag}\left[0\mathbf{I},\tilde{\mathbf{P}_{i}}^{-1},\raisebox{1.00006pt}{$\infty$}\mathbf{I}\right]+\text{diag}\left[t_{i}\left(\begin{array}[]{ll}\boldsymbol{\Sigma}^{\text{\tiny{{obs}}}}_{i}&\boldsymbol{\Sigma}^{\text{\tiny{{ol}}}}_{i}\\ -&\boldsymbol{\Sigma}^{\text{\tiny{{lat}}}}_{i}\\ \end{array}\right),0\mathbf{I}\right]\right)^{-}\mathbf{C}_{i}^{t}$$

$$={\mathbf{C}}_i{\left(\text{diag}\infty \right)}^{-}{\mathbf{C}}_i^t$$

$$={\mathbf{C}}_i\text{diag}\left[{\mathbf{T}}^{-1},0\mathbf{I}\right]{\mathbf{C}}_i^t,$$

$$\mathbf{T}=\text{diag}\left[0\mathbf{I},\tilde{\mathbf{P}_{i}}^{-1}\right]+t_{i}\left(\begin{array}[]{ll}\boldsymbol{\Sigma}^{\text{\tiny{{obs}}}}_{i}&\boldsymbol{\Sigma}^{\text{\tiny{{ol}}}}_{i}\\ -&\boldsymbol{\Sigma}^{\text{\tiny{{lat}}}}_{i}\\ \end{array}\right)=\left(\begin{array}[]{cc}t_{i}\boldsymbol{\Sigma}^{\text{\tiny{{obs}}}}_{i}&t_{i}\boldsymbol{\Sigma}^{\text{\tiny{{ol}}}}_{i}\\ -&\tilde{\mathbf{P}_{i}}^{-1}+t_{i}\boldsymbol{\Sigma}^{\text{\tiny{{lat}}}}_{i}\end{array}\right).$$

Full text

Version dated:

Inferring Phenotypic Trait Evolution on Large Trees

With Many Incomplete Measurements

Gabriel Hassler1,

Max R. Tolkoff2,

William L. Allen3,

Lam Si Tung Ho4,

Philippe Lemey5,

and Marc A. Suchard1,2,6

1*Department of Biomathematics, David Geffen School of Medicine at UCLA, University of California, Los Angeles, United States

2Department of Biostatistics, Jonathan and Karin Fielding School of Public Health, University of California, Los Angeles, United States*

3*Department of Biosciences, Swansea University, Swansea, United Kingdom

4Department of Mathematics and Statistics, Dalhousie University, Halifax, Nova Scotia, Canada

5Department of Microbiology and Immunology, Rega Institute, KU Leuven, Leuven, Belgium

6Department of Human Genetics, David Geffen School of Medicine at UCLA, Universtiy of California, Los Angeles, United States*

Corresponding author: Marc A. Suchard, Departments of Biostatistics, Biomathematics, and Human Genetics, University of California, Los Angeles, 695 Charles E. Young Dr., South, Los Angeles, CA 90095-7088, USA; E-mail: msuchard@ucla.edu

Abstract

Comparative biologists are often interested in inferring covariation between multiple biological traits sampled across numerous related taxa. To properly study these relationships, we must control for the shared evolutionary history of the taxa to avoid spurious inference. Existing control techniques almost universally scale poorly as the number of taxa increases. An additional challenge arises as obtaining a full suite of measurements becomes increasingly difficult with increasing taxa. This typically necessitates data imputation or integration that further exacerbates scalability. We propose an inference technique that integrates out missing measurements analytically and scales linearly with the number of taxa by using a post-order traversal algorithm under a multivariate Brownian diffusion (MBD) model to characterize trait evolution. We further exploit this technique to extend the MBD model to account for sampling error or non-heritable residual variance. We test these methods to examine mammalian life history traits, prokaryotic genomic and phenotypic traits, and HIV infection traits. We find computational efficiency increases that top two orders-of-magnitude over current best practices. While we focus on the utility of this algorithm in phylogenetic comparative methods, our approach generalizes to solve long-standing challenges in computing the likelihood for matrix-normal and multivariate normal distributions with missing data at scale.

Keywords:

Bayesian inference, matrix-normal, missing data, phylogenetics

.   ntroduction

Phylogenetic comparative methods explore the relationships between different biological phenotypes across sets of organisms. To properly understand these phenotypic trait relationships, methods must adjust for the shared evolutionary history of the taxa (Felsenstein 1985). Molecular sequences from emerging sequencing technology and high-throughput biological experimentation enable such phylogenetic adjustment for rapidly growing numbers of taxa and increasing numbers of trait measurements. Comparative studies incorporating dense taxonomic sampling create the potential for new research into general patterns in phenotypic evolution, key differences between subgroups and the relationship between phenotypic and genetic evolutionary dynamics. Unfortunately, many phylogenetic comparative methods remain poorly equipped to handle these research questions at scale.

Popular methods often assume an underlying Brownian diffusion process acts along each branch of a phylogenetic tree, such that the traits are multivariate normally distributed. Revell (2012) and Adams (2014), for example, parameterize this distribution in terms of a highly-structured variance-covariance matrix that characterizes the tree and trait covariation. Computational work to invert this matrix to evaluate the multivariate normal likelihood scales cubically with the number of taxa. This work stands even more troublesome when the phylogenetic tree remains unknown and requires joint inference with the trait process, necessitating repeated inversion. Freckleton (2012), Pybus et al. (2012), and Ho and Ané (2014) all independently develop algorithms that take advantage of the matrix-normal structure of the data under the MBD model to evaluate the likelihood. Using the tree structure, these algorithms then scale linearly with the number of taxa with complete data, but this ideal run-time currently stumbles when trait measurements are missing.

As the number of taxa grows large, measuring a complete suite of traits for all taxa becomes increasingly challenging. Recent solutions to this problem in phylogenetics include those proposed by Tolkoff et al. (2017), who use Markov chain Monte Carlo (MCMC) to numerically integrate via sampling the missing data, and Bastide et al. (2018), who exploit expectation-maximization (EM) to impute missing values and all unobserved ancestral node traits. Both lines of work, however, require iterative manipulation of the likelihood function on a per-taxon basis. Current sampling methods scale at best quadratically with the number of taxa, and likewise the number of ancestral node trait values that require imputation is directly related to the number of taxa. As a consequence, these methods remain computationally prohibitive for large trees.

In this paper, we reformulate evaluation of the data likelihood function under a Brownian diffusion process on a tree such that we achieve the marginalized likelihood of the observed trait measurements only. This innovation arises from thinking about observed tip traits as multivariate normally distributed with infinite precision in their sampling, while missing traits have zero precision, and appropriately propagating these precisions up the tree through dynamic programming involving an unusual matrix pseudo-inverse definition. This pseudo-inverse finds similar use, but independent discovery, in Bastide et al. (2018). Unlike previous approaches, the integration avoids EM iteration making simultaneous inference with the phylogeny practical and enables researchers to analyze all available measurements when inferring the trait relationships. Surprisingly, we can still evaluate the observed-data likelihood in linear time with respect to the number of taxa. The price to be paid is that computation now scales cubically, rather than quadratically, in the number of traits. This remains a small price since the number of taxa is often orders-of-magnitude larger than the number of traits. It is also notable that this method has applications beyond phylogenetic comparative methods and can be used more generally in a special class of matrix-normal and multivariate normal distributions with missing data. This has been a long standing problem in statistics since at least the 1930’s (Wilks 1932), with more recent work by Dominici et al. (2000); Cantet et al. (2004); Allen and Tibshirani (2010); and Glanz and Carvalho (2018). One important limitation to our approach is that it assumes data are missing at random (Little and Rubin 1987) which is inappropriate for many data sets.

We also demonstrate how this framework can be easily extended to incorporate residual variance in the MBD model, which is only one of many possible model extensions. Our strategy of analytically marginalizing the observed data likelihood extends seamlessly to this and other model extensions and allows for efficient inference on these models while maintaining likelihood computations that scale linearly with the number of taxa. These extensions open up lines of inquiry not available in the simple MBD model. In particular, including residual variance in the model enables inference of phylogenetic heritability.

We demonstrate the broad utility of our algorithm to compute the marginalized likelihood through three examples. First, we examine covariation in mammalian life history traits using data on 3690 taxa from the PanTHERIA ecological database (Jones et al. 2009). Second, we use our new efficient algorithm to simultaneously evaluate several theories regarding prokaryotic evolutionary theory. We use data from NCBI Genome and a recent study by Goberna and Verdú (2016), along with matching 16S sequences from the ARB Silva Database (Ludwig et al. 2004), to jointly infer both the phylogenetic tree and evolutionary correlation between several prokaryotic genotypic and phenotypic traits. Finally, we apply our multivariate residual variance model extension to data presented by Blanquart et al. (2017) concerning HIV virulence to evaluate the heritability of HIV viral load and CD4 T-cell decline. We compare the computation speed of our analytical integration method against current best-practice methods and observed increases in speed that top two orders-of-magnitude.

.   henotypic diffusion on trees

Consider a data-complete collection $\mathbf{Y}=\left(\mathbf{Y}_{1},\ldots,\mathbf{Y}_{N}\right)^{t}$ where $\mathbf{Y}_{i}=\left(Y_{i1},\ldots,Y_{iP}\right)^{t}$ of $P$ real-valued phenotypic traits measured across $N$ biological taxa. Relating the taxa stands a known and fixed or unknown and random phylogeny $\mathscr{F}$ that is a bifurcating, directed acyclic graph whose $2N-1$ vertices originate with a degree-2 root node $\nu_{2N-1}$ and terminate with degree-1 tip nodes $\left(\nu_{1},\ldots,\nu_{N}\right)$ that correspond to the $N$ taxa. Linking vertices are edge weights or branch lengths $(t_{1},\ldots,t_{2N-2})$. Let $\mathbf{X}_{k}=\left(X_{k1},\ldots,X_{kP}\right)$ be latent values of the traits at node $\nu_{k}$ on the tree for $k=1,\ldots,2N-1$. For tip nodes $i=1,\ldots,N$, we posit a stochastic link $p\left(\mathbf{Y}_{i}\left|\mathbf{X}_{i}\right.\right)$ where $\mathbf{Y}_{i}$ is drawn from some distribution parameterized by $\mathbf{X}_{i}$ and other hyperparameters (see Figure 1). Comparative methods standardly assume that the density $p\left(\mathbf{Y}_{i}\left|\mathbf{X}_{i}\right.\right)$ is degenerate at $\mathbf{X}_{i}$ (i.e. $\mathbf{Y}_{i}=\mathbf{X}_{i}$ with probability 1), but we relax this assumption in future sections.

The most common phenotypic model of evolution (Felsenstein 1985) assumes a multivariate Brownian diffusion process acts conditionally independently along each branch generating a multivariate normal (MVN) increment,

[TABLE]

centered around the realized value $\mathbf{X}_{\text{pa}(k)}$ at its parent node and variance proportional to an estimable $P\times P$ positive-definite matrix $\boldsymbol{\Sigma}$. Since the trait values at the root are also unknown, Pybus et al. (2012) suggest further assuming $\mathbf{X}_{2N-1}\sim\text{MVN}\left(\boldsymbol{\mu}_{0},\kappa_{0}^{-1}\boldsymbol{\Sigma}\right)$ with fixed prior mean $\boldsymbol{\mu}_{0}$ and sample-size $\kappa_{0}$.

2.1   Computation of Observed Data Likelihood

When there are no missing data and under our standard assumption that $p\left(\mathbf{Y}_{i}\left|\mathbf{X}_{i}\right.\right)$ is degenerate, integrating out unobserved internal and root node traits leads to a seemingly simple expression for the data likelihood $p(\mathbf{Y}\,|\,\boldsymbol{\Sigma},\mathscr{F},\boldsymbol{\mu}_{0},\kappa_{0})$ (Freckleton 2012; Vrancken et al. 2015). Namely, $\mathbf{Y}$ is matrix-normal (MN) distributed around mean $\mathbf{1}_{N}\boldsymbol{\mu}_{0}^{t}$, with across-row variance $\boldsymbol{\Upsilon}+\kappa_{0}^{-1}\mathbf{J}_{N}$ and across-column variance $\boldsymbol{\Sigma}$, where $\mathbf{1}_{N}$ is a vector of length $N$ populated by ones, $\mathbf{J}_{N}=\mathbf{1}_{N}\mathbf{1}_{N}^{t}$, and $\boldsymbol{\Upsilon}$ is a deterministic function of $\mathscr{F}$. Specifically, element $\Upsilon_{ii^{\prime}}$ measures shared evolutionary history and equals the sum of the branch lengths from the root to the most recent common ancestral node of taxa $i$ and $i^{\prime}$ when $i\neq i^{\prime}$ or the sum of the branch lengths from the root to taxon $i$ otherwise. For example, in Figure 1, $\Upsilon_{12}=t_{4}$ and $\Upsilon_{11}=t_{1}+t_{4}$. One can evaluate this highly structured matrix-normal likelihood function with computational complexity $\mathcal{O}\left(NP^{2}\right)$ given the acyclic nature of $\mathscr{F}$. When some data points are missing, however, the observed-data likelihood is no longer matrix-normal and new approaches are needed. This becomes increasingly urgent as the prevalence of missing observations grows with the size of trait data sets. In this context we wish to compute

[TABLE]

where $\mathbf{Y}^{\text{\tiny obs}}$ and $\mathbf{Y}^{\text{\tiny mis}}$ contain the observed and missing trait values, respectively.

The two simplest strategies for calculating the observed-data likelihood are, unfortunately, computationally prohibitive for most large problems. One such solution forfeits the MN structure of the data in favor a simple expression of the observed-data likelihood. This strategy uses the fact that the matrix-normal distribution of $\mathbf{Y}$ can also be expressed as

[TABLE]

using the Kronecker product $\otimes$. Assuming data are missing at random (Little and Rubin 1987), one can simply remove the rows and columns of $\mbox{vec}\left[\mathbf{1}_{N}\boldsymbol{\mu}_{0}^{t}\right]$ and $\boldsymbol{\Sigma}\otimes\boldsymbol{\Upsilon}$ corresponding to the missing data and compute the likelihood for this $NP-M^{\prime}$ dimensional MVN distribution, where $M^{\prime}$ is the number of missing measurements. This likelihood calculation carries the onerous computational complexity $\mathcal{O}\left(\left(NP-M^{\prime}\right)^{3}\right)$. Alternatively, from a Bayesian perspective, one could numerically integrate out the missing data by treating each missing data point as an unknown model parameter and employing MCMC to sample each value. This strategy restores the matrix-normal structure, but requires the likelihood be evaluated each time one samples a missing data point. This results in computation complexity of at least $\mathcal{O}\left(NP^{2}M\right)$, where $M$ is the number of taxa with missing measurements. Because $M$ often scales with $N$, this method remains prohibitively slow for many data sets with large $N$. Our goal is to integrate out these missing values analytically using a dynamic programming algorithm in order to bring run time down to a much more manageable $\mathcal{O}\left(NP^{3}\right)$.

2.1.1   Missing Data Definitions and Operations

To develop our algorithm, we first introduce some useful abstractions and notation. At each tip in $\mathscr{F}$, information about each of the $P$ traits comes in one of three forms: a trait value may be directly observed, latent, or completely missing. When directly observed, we posit without loss of generality that the value arises from a normal distribution centered at the observed value with infinite precision. We assume that trait data that arise from latent values are jointly multivariate normally distributed about the unknown latent values with known or estimable precision. Finally, a completely missing value arises also without loss of generality from a normal distribution centered at 0 with zero precision. To formalize this, for tip $i=1,\ldots,N$, we construct a permutation matrix $\mathbf{C}_{i}$ that groups traits in directly observed, latent, and completely missing order and populate a pseudo-precision matrix

[TABLE]

where $\text{diag}\left[\cdot\right]$ is a function that arranges its constituent elements into block-diagonal form and $\mathbf{R}_{i}$ is the latent block precision. Note that any block may be 0-dimensional. This construction arbitrarily forces off-diagonal elements of $\mathbf{P}_{i}$ involving directly observed and completely missing traits to equal 0 and plays an important role in simplifying computations.

We additionally define a series of operations that we will find useful for defining this algorithm. We define the pseudo-inverse

[TABLE]

We define the pseudo-determinant $\hat{\text{det}}\left(\right)$ as the product of the non-zero singular values. We also define the matrix $\boldsymbol{\delta}_{i}=\text{diag}\left[\delta_{i1},\ldots,\delta_{iP}\right]$ for $i=1,\ldots,N$, where $\delta_{ij}$ is an indicator variable which takes a value of $1$ if trait $Y_{ij}$ is observed or latent and [math] if it is missing. Lastly, we define the possibly degenerate multivariate normal density function

[TABLE]

for some argument $\mathbf{z}$, mean $\boldsymbol{\mu}$ and precision $\mathbf{P}$ of appropriate dimensions.

2.1.2   Post-Order Observed Data Likelihood Algorithm

Our goal is to efficiently compute the likelihood $p(\mathbf{Y}^{\text{\tiny obs}}\,|\,\boldsymbol{\Sigma},\mathscr{F},\boldsymbol{\mu}_{0},\kappa_{0})$. Following from Pybus et al. (2012), we perform a post-order traversal where we calculate the observed-data partial likelihood $p(\mathbf{Y}^{\text{\tiny obs}}_{\lfloor{k}\rfloor}\,|\,\mathbf{X}_{k},\boldsymbol{\Sigma},\mathscr{F})$ at each node $\nu_{k}$ where $\mathbf{Y}^{\text{\tiny obs}}_{\lfloor{k}\rfloor}$ is the observed data restricted to all descendants of node $k$ on the tree. For example, in Figure 1, $\mathbf{Y}^{\text{\tiny obs}}_{\lfloor{4}\rfloor}=\{\mathbf{Y}^{\text{\tiny obs}}_{1},\mathbf{Y}^{\text{\tiny obs}}_{2}\}$.

We posit that, given an appropriate stochastic link function $p\left(\mathbf{Y}_{i}\left|\mathbf{X}_{i}\right.\right)$, we can express the observed-data partial likelihood as

[TABLE]

for all nodes $k=1,\ldots,2N-1$ and some remainder $r_{k}$, mean $\mathbf{m}_{k}$, and precision $\mathbf{P}_{k}$. Given a parent node $\ell$ with children $j$ and $k$, let us assume we can express the observed-data likelihood of $\mathbf{Y}^{\text{\tiny obs}}_{\lfloor{j}\rfloor}$ and $\mathbf{Y}^{\text{\tiny obs}}_{\lfloor{k}\rfloor}$ as in Equation 6. Conditioning on $\mathbf{X}_{\ell}$, we can compute

[TABLE]

as $\mathbf{Y}^{\text{\tiny obs}}_{\lfloor{j}\rfloor}$ and $\mathbf{Y}^{\text{\tiny obs}}_{\lfloor{k}\rfloor}$ are conditionally independent given $\mathbf{X}_{\ell}$. Using Equations 1 and 6, we form

[TABLE]

where the branch-deflated pseudo-precision $\mathbf{P}_{j}^{\star}=\left(\mathbf{P}_{j}^{-}+t_{j}\boldsymbol{\delta}_{j}\boldsymbol{\Sigma}\boldsymbol{\delta}_{j}\right)^{-}$. See Supplementary Information (SI) Section 1 for details on computing this pseudo-inverse. We use the results of Equation 8 in Equation 7 to compute the partial log-likelihood

[TABLE]

where $\mathbf{P}_{\ell}=\mathbf{P}_{j}^{\star}+\mathbf{P}_{k}^{\star}$, $\mathbf{m}_{\ell}$ is a solution to $\mathbf{P}_{\ell}\mathbf{m}_{\ell}=\mathbf{P}_{j}^{\star}\mathbf{m}_{j}+\mathbf{P}_{k}^{\star}\mathbf{m}_{k}$, and

[TABLE]

Note that the change of informative dimensions $\Delta_{jk\ell}=\text{rank}\left(\mathbf{P}_{j}^{\star}\right)+\text{rank}\left(\mathbf{P}_{k}^{\star}\right)-\text{rank}\left(\mathbf{P}_{\ell}\right)$. We update $\boldsymbol{\delta}_{\ell}=\boldsymbol{\delta}_{j}\lor\boldsymbol{\delta}_{k},$ where $\lor$ is the element-wise “logical or” operation.

Our algorithm initializes $r_{i}$, $\mathbf{m}_{i}$, and $\mathbf{P}_{i}$ such that $p(\mathbf{Y}^{\text{\tiny obs}}_{i}\,|\,\mathbf{X}_{i})=r_{i}\hat{\phi}\left(\mathbf{X}_{i};\mathbf{m}_{i},\mathbf{P}_{i}\right)$ at the tips of the tree. For the standard assumption that $\mathbf{Y}_{i}=\mathbf{X}_{i}$, we have $r_{i}=1$, $\mathbf{m}_{i}=\mathbf{C}_{i}\left[\mathbf{Y}^{\text{\tiny obs}}_{i},\mathbf{0}\right]$, and \mathbf{P}_{i}=\mathbf{C}_{i}\,\text{diag}\left[\raisebox{1.00006pt}{\infty}\mathbf{I},0\mathbf{I}\right]\mathbf{C}_{i}^{t}. We perform a post-order traversal of the tree computing $\mathbf{m}_{\ell},\mathbf{P}_{\ell},$ and $r_{\ell}$ for internal nodes $\ell=N+1,\ldots,2N-2$ using the already-computed node remainders, means, and precisions for their respective child nodes. At the root, $\mathbf{Y}^{\text{\tiny obs}}_{\lfloor{2N-1}\rfloor}=\mathbf{Y}^{\text{\tiny obs}}$ and we return the observed-data log-likelihood

[TABLE]

where $\mathbf{P}_{\text{full}}=\mathbf{P}_{2N-1}+\kappa_{0}\boldsymbol{\Sigma}^{-1}$ and $\mathbf{m}_{\text{full}}=\mathbf{P}_{\text{full}}^{-1}\left(\mathbf{P}_{2N-1}\mathbf{m}_{2N-1}+\kappa_{0}\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}_{0}\right)$. The integral evaluates to one, leaving the observed-data log-likelihood

[TABLE]

This tree traversal visits each node in $\mathscr{F}$ exactly once and inverts a $P\times P$ matrix each time, resulting in an overall computational complexity of $\mathcal{O}\left(NP^{3}\right)$ for each likelihood evaluation.

2.2   Inference

The primary parameter of scientific interest is the diffusion variance $\boldsymbol{\Sigma}$. We are also often interested in additional hyper-parameters $\boldsymbol{\theta}$ related to the stochastic link function $p\left(\mathbf{Y}_{i}\left|\mathbf{X}_{i}\right.\right)$. In cases where the tree structure is unknown, we use sequence data $\mathbf{S}$ to simultaneously infer $\mathscr{F}$. As such, from a Bayesian perspective, we are interested in approximating

[TABLE]

for inference. We place a $\text{Wishart}_{P}\left(\boldsymbol{\Lambda}_{0},\nu\right)$ prior on $\boldsymbol{\Sigma}^{-1}$, where $\boldsymbol{\Lambda}_{0}$ is a $P\times P$ rate matrix. The prior on $\boldsymbol{\theta}$ depends the problem of interest, and there are many ways to specify $p(\mathscr{F},\mathbf{S})$ (see Suchard et al. 2018). To approximate the posterior distributions via MCMC simulation, we apply a random scan Metropolis-within-Gibbs (Liu et al. 1995) approach by which we sample parameter blocks one at a time at random from their full conditional distribution.

Let $\mathbf{X}=\left(\mathbf{X}_{1}^{t},\cdots,\mathbf{X}_{N}^{t}\right)^{t}$ be the latent trait values at the tips of the phylogeny. The conjugate $\text{Wishart}_{P}\left(\boldsymbol{\Lambda}_{0},\nu\right)$ prior on $\boldsymbol{\Sigma}^{-1}$ implies that

[TABLE]

We apply the post-order computation method proposed by Ho and Ané (2014) to compute $\left(\mathbf{X}-\mathbf{J}_{N}\boldsymbol{\mu}_{0}^{t}\right)^{t}\left(\boldsymbol{\Upsilon}+\frac{1}{\kappa_{0}}\mathbf{J}\right)^{-1}\left(\mathbf{X}-\mathbf{J}_{N}\boldsymbol{\mu}_{0}^{t}\right)$, which has computational complexity $\mathcal{O}\left(NP^{2}\right)$. When $\mathbf{X}$ are known (i.e. when there are no missing values and $p\left(\mathbf{Y}_{i}\left|\mathbf{X}_{i}\right.\right)$ is degenerate at $\mathbf{X}_{i}$), we can sample from the distribution in Equation 14 immediately without any additional steps. However, if either assumption is violated, we must first draw from the full conditional distribution of $\mathbf{X}$ via the data augmentation algorithm described below.

2.2.1   Pre-Order Missing Data Augmentation Algorithm

To sample jointly from the full conditional of $\mathbf{X}=\left(\mathbf{X}_{1},\ldots,\mathbf{X}_{N}\right)^{t}$, we draw on the calculations made in Section 2.1.2 and perform a pre-order traversal of the tree. Note that we omit explicit dependence on the parameters $\boldsymbol{\Sigma},\mathscr{F}$, and $\boldsymbol{\theta}$ in all calculations below for clarity. Starting at the root, $\mathbf{X}_{2N-1}$, we draw from $\left.{\mathbf{X}_{2N-1}}\right|{\mathbf{Y}^{\text{\tiny obs}},\boldsymbol{\mu}_{0},\kappa_{0}}$. Using Bayes’ rule and Equation 11, we see that

[TABLE]

After sampling the root traits from their full conditional, we continue the traversal of the tree where we sample each node $\mathbf{X}_{j}$ conditional on its (previously sampled) parent node $\mathbf{X}_{\text{pa}(j)}$ and the observed data below node $j$ $\mathbf{Y}^{\text{\tiny obs}}_{\lfloor{j}\rfloor}$ for $j=1,\ldots,2N-2$. For the internal nodes, we compute $p(\mathbf{X}_{j}\,|\,\mathbf{Y}^{\text{\tiny obs}}_{\lfloor{j}\rfloor},\mathbf{X}_{\text{pa}(j)})$ as follows:

[TABLE]

where $\mathbf{Q}_{j}=\mathbf{P}_{j}+\left(t_{j}\boldsymbol{\Sigma}\right)^{-1}$ and $\mathbf{n}_{j}=\mathbf{Q}_{j}^{-1}\left(\mathbf{P}_{j}\mathbf{m}_{j}+\left(t_{j}\boldsymbol{\Sigma}\right)^{-1}\mathbf{X}_{\text{pa}(j)}\right)$. This implies $\left.{\mathbf{X}_{j}}\right|{\mathbf{Y}^{\text{\tiny obs}}_{\lfloor{j}\rfloor},\mathbf{X}_{\text{pa}(j)}}\sim\text{MVN}\left(\mathbf{n}_{j},\mathbf{Q}_{j}\right)$, and we sample $\mathbf{X}_{j}$ from this distribution.

At the tips, we employ one of two techniques depending on the specific model. Under our standard assumption (i.e. $\mathbf{X}_{i}=\mathbf{Y}_{i}$ with probability 1), we partition the precision $\boldsymbol{\Sigma}^{-1}$ and trait values $\mathbf{X}_{i}$ and $\mathbf{X}_{\text{pa}(i)}$ such that

[TABLE]

and draw from $\left.{\mathbf{X}^{\text{\tiny mis}}_{i}}\right|{\mathbf{Y}^{\text{\tiny obs}}_{i},\mathbf{X}_{\text{pa}(i)}}\sim\text{MVN}\left(\mathbf{X}^{\text{\tiny mis}}_{\text{pa}(i)}+{\mathbf{S}^{\text{\tiny{{mis}}}}_{i}}^{-1}\mathbf{S}^{\text{\tiny{{mo}}}}_{i}\left(\mathbf{X}^{\text{\tiny obs}}_{\text{pa}(i)}-\mathbf{X}^{\text{\tiny obs}}_{i}\right),\frac{1}{t_{i}}\mathbf{S}^{\text{\tiny{{mis}}}}_{i}\right)$ for $i=1,\ldots,N$. For cases where $p\left(\mathbf{Y}_{i}\left|\mathbf{X}_{i}\right.\right)$ is non-degenerate, we simply use Equation 16 to sample from $\left.{\mathbf{X}_{i}}\right|{\mathbf{Y}^{\text{\tiny obs}}_{i},\mathbf{X}_{\text{pa}(i)}}$. Once we have sampled $\mathbf{X}\left|\mathbf{Y}^{\text{\tiny obs}},\boldsymbol{\Sigma},\mathscr{F},\boldsymbol{\theta}\right.$, we can draw from the full conditional distribution of $\boldsymbol{\Sigma}^{-1}$ via Equation 14. This pre-order data augmentation procedure requires a single $P\times P$ matrix inversion at each of the ${2N-1}$ nodes in the tree, resulting in overall computational complexity $\mathcal{O}\left(NP^{3}\right)$.

.   odel Extension: Residual Variance

We extend the MBD model of phenotypic evolution to include multivariate normal residual variance at each of the tips. Under this model, we assume

[TABLE]

where $\boldsymbol{\Gamma}$ is a $P\times P$ precision matrix. Under this model, the vectorization of $\mathbf{Y}$ is MVN-distributed with $NP\times NP$ variance-covariance matrix $\boldsymbol{\Sigma}\otimes\left(\boldsymbol{\Upsilon}+\kappa_{0}^{-1}\mathbf{J}_{N}\right)+\boldsymbol{\Gamma}^{-1}\otimes\mathbf{I}_{N}$ where $\mathbf{I}_{N}$ is the $N\times N$ identity matrix. Unlike the case where $\mathbf{Y}_{i}=\mathbf{X}_{i}$, $\mathbf{Y}$ cannot be expressed as matrix-normal even in the data-complete case because the variance cannot be expressed as the Kronecker product of two matrices. As such, our post-order likelihood computation algorithm is useful for this extended model, even when there are no missing data points.

3.1   Inference of Residual Variance

Similar to our inference of $\boldsymbol{\Sigma}$ in the diffusion process, we place a conjugate $\text{Wishart}\left(\boldsymbol{\Lambda}_{s},\nu_{s}\right)$ prior on $\boldsymbol{\Gamma}$ using the rate parameterization. This yields the full conditional distribution

[TABLE]

Because $\mathbf{X}$ is latent in this model, each time we update $\boldsymbol{\Gamma}$ we first draw from the full conditional posterior of $\mathbf{X}$ using the algorithm described in Section 2.2.1. For cases where $\mathbf{Y}$ is not completely observed, we must perform an additional data augmentation step where we draw from $\mathbf{Y}^{\text{\tiny mis}}\left|\mathbf{Y}^{\text{\tiny obs}},\mathbf{X},\boldsymbol{\Gamma}\right.$. To do this, we decompose the sampling precision matrix into blocks such that

[TABLE]

From Equation 18, we see that

[TABLE]

As such, we can directly sample $\mathbf{Y}^{\text{\tiny mis}}_{i}$ from its full conditional above and update $\mathbf{Y}_{i}=\mathbf{C}_{i}\left[\mathbf{Y}^{\text{\tiny obs}}_{i},\mathbf{Y}^{\text{\tiny mis}}_{i}\right]^{t}$ for $i=1,\ldots,N$. This process also has computational complexity $\mathcal{O}\left(NP^{3}\right)$.

Note that we can draw from the joint full conditional of $\boldsymbol{\Sigma}$ and $\boldsymbol{\Gamma}$ by performing a single pre-order data augmentation where we draw from $p(\mathbf{X},\mathbf{Y}^{\text{\tiny mis}}\,|\,\boldsymbol{\Sigma},\boldsymbol{\Gamma})$ and subsequently draw from $p(\boldsymbol{\Sigma},\boldsymbol{\Gamma}\,|\,\mathbf{X},\mathbf{Y}^{\text{\tiny mis}})=p(\boldsymbol{\Sigma}\,|\,\mathbf{X})p(\boldsymbol{\Gamma}\,|\,\mathbf{X},\mathbf{Y})$. These distributions are conditionally independent due to the fact that $\mathbf{X}$ and $\mathbf{X}-\mathbf{Y}$ are independent by construction. This procedure effectively halves the computation time as we only need to perform a single post-order likelihood computation/pre-order data augmentation step to sample both $\boldsymbol{\Sigma}$ and $\boldsymbol{\Gamma}$, rather than each time we sample one.

3.2   Heritability Statistic

The residual variance extension enables us to estimate phenotypic heritability over evolutionary time. We use a definition analogous to the broad-sense heritability in statistical genetics (see Visscher et al. 2008). Namely, we seek to quantify the proportion of variance in a trait attributable to the Brownian diffusion process on the phylogeny (as opposed to the residual variance). Note that we are primarily interested in heritability in the HIV example below, for which we use data from a recent paper by Blanquart et al. (2017). As such, we use a multivariate generalization of the heritability statistic from that paper. Specifically, we estimate phylogenetic heritability by taking the expectation of the empirical sample variance under our extended model. We define the $P\times P$ empirical covariance matrix as

[TABLE]

where $\bar{\mathbf{y}}=\frac{1}{N}\sum_{i=1}^{N}\mathbf{Y}_{i}=\frac{1}{N}\mathbf{Y}^{t}\mathbf{1}_{N}$ and $\bar{\mathbf{Y}}=\mathbf{1}_{N}\bar{\mathbf{y}}^{t}=\frac{1}{N}\mathbf{J}_{N}\mathbf{Y}$. The expectation of this quantity reduces to the following expression (see SI Section 2 for details):

[TABLE]

Because $\mathbb{E}\left[\mathbf{S}^{2}\left({\mathbf{Y}}\right)\right]$ is a linear combination of $\boldsymbol{\Sigma}$ and $\boldsymbol{\Gamma}^{-1}$, we propose the $P\times P$ heritability matrix $\mathbf{H}=\{h_{kl}\}$ with entries

[TABLE]

where $c_{\sigma}=\frac{1}{N}\mbox{tr}\left[\boldsymbol{\Upsilon}\right]-\frac{1}{{N}^{2}}\mathbf{1}_{N}^{t}\boldsymbol{\Upsilon}\mathbf{1}_{N}$ and $c_{\gamma}=\frac{N-1}{N}$. The diagonal entries $h_{kk}=h^{2}_{k}$ represent the marginal phylogenetic heritability of each trait, and the off-diagonal entries represent pair-wise co-heritability (Falconer 1960, chap. 19) between traits. Note that naive computation of $c_{\sigma}$ relies on constructing $\boldsymbol{\Upsilon}$ that has complexity $\mathcal{O}\left(N^{2}\right)$. SI Section 3 describes a novel post-order algorithm that computes $c_{\sigma}$ in $\mathcal{O}\left(N\right)$ time without explicitly constructing $\boldsymbol{\Upsilon}$.

While the breadth of research in heritability is extensive across both statistical genetics and phylogenetics (see in particular the recent paper by Mitov and Stadler 2018), we choose the same heritability statistic as used by Blanquart et al. (2017) for direct comparison with their analysis. That being said, our methods could be readily adapted to approximate the posterior distribution of several of the alternative heritability statistics presented in Mitov and Stadler (2018). Additionally, our pre-order data augmentation procedure allows us to generate samples directly from the posterior of the latent trip traits $\mathbf{X}$, from which we can directly compute the genetic covariance $\mathbf{S}^{2}\left({\mathbf{X}}\right)$ rather than relying on expectations.

.   esearch Materials

We have implemented these methods in the development version of BEAST (Drummond et al. 2012). The XML files for all analyses and instructions for running them are available at https://github.com/suchard-group/missing_traits_paper.

.   omputational Efficiency

Our method dramatically increases computational efficiency over the current best-practice method. This latter procedure, developed by Cybis et al. (2015), treats the missing and latent values of $\mathbf{X}$ as unknown parameters and numerically integrates them out by placing a Gibbs sampler on each tip $\mathbf{X}_{i}$ that draws from its full conditional distribution $p\left({\mathbf{X}_{i}}\left|{\mathbf{Y}_{i},\mathbf{X}_{\lceil{i}\rceil}}\right.\right)$ for $i=1,\ldots,N$ where $\mathbf{X}_{\lceil{i}\rceil}=\mathbf{X}\backslash\mathbf{X}_{i}$. Because the full conditional distribution of $\mathbf{X}_{i}$ relies on the other missing and latent values in $\mathbf{X}$, we sample each tip individually. The advantage of this is that the likelihood calculation, the Gibbs sampler of the diffusion variance $\boldsymbol{\Sigma}$, and the data augmentation procedure for each tip all have complexity $\mathcal{O}\left(NP^{2}\right)$ rather than our $\mathcal{O}\left(NP^{3}\right)$. As such, this numerical integration procedure has overall complexity $\mathcal{O}\left(MNP^{2}\right)$ where $M$ is the number of tips with missing or latent values. For any extended model where $p\left(\mathbf{Y}_{i}\left|\mathbf{X}_{i}\right.\right)$ is not degenerate at $\mathbf{X}_{i}$, all values of $\mathbf{X}$ are latent and $M=N$.

We formalize our comparison by computing the median and minimum effective sample size (ESS) per hour for all parameters of interest under both our analytical integration method and the sampling method discussed above. We also compute the ESS per sample and samples per hour to understand how our improved method influences both the autocorrelation between MCMC samples and the amount of computation work required to generate a single draw from the posterior. We define the number of samples as the number of states in which the MCMC simulation updates the parameters of interest (as opposed to missing trait values). Table 1 presents the results of our efficiency comparisons.

.   pplications

6.1   Mammalian Life History

A major task for life history theory is to understand the ecological and evolutionary significance of correlation between life history traits such as age at sexual maturity, the number of offspring per reproductive event, and reproductive lifespan (Roff 2002). Establishing patterns of such correlation grants insight into whether life history variation between individuals, populations or species is consistent with pace-of-life theory (Reynolds 2003; Réale et al. 2010). This theory predicts that ‘fast’ traits such as early maturity, large broods, small offspring, frequent reproduction and a short lifespan are positively associated with each other as a consequence of organisms pursuing strategies that prioritize either current or future reproduction. Existing approaches using comparative life history data to investigate fast-slow trait covariation patterns (e.g. mammals: Bielby et al. 2007; hymenoptera: Blackburn 1991; lizards: Clobert et al. 1998; birds: Sæther and Bakke 2000; plants: Salguro-Gómez 2017; fish: Wiedmann et al. 2014) generally support the fast-slow hypothesis; however, results are rarely consistent across taxa. This may reflect important taxonomic differences in life history evolution, but there is concern that differences are an artifact of different methodologies (Jeschke and Kokko 2009).

One key limitation is that previous methods have required complete data for each species. As complete measurements across a rich suite of varied life history traits are not yet available for most species, this means that researchers must choose to either reduce the number of traits or reduce the number of species included in analyses. By integrating out missing traits, we resolve this issue and analyze the life history dataset used in Capellini et al. (2015), which is based largely on the final PanTHERIA dataset (Jones et al. 2009), supplemented with measurements from Ernest (2003) and additional sources. Our analysis includes all the variables analyzed by Bielby et al. (2007) (gestation length, weaning age, neonatal body mass, litter size, litter frequency, and age at first birth) plus reproductive lifespan (maximum lifespan minus age at first birth). We include female body mass as a trait rather than analyze size-corrected residuals and log-transform and standardize all traits prior to analysis. The analysis assumes the phylogeny of Fritz et al. (2009) that remains the most complete phylogeny for mammals. In total, 3690 species in the phylogeny have measurement of at least one trait and are included. Table 2 reports the number of species with measurements for each trait. Only 518 species have complete data on all 8 traits; thus the ability to include species with partially missing traits enables inclusion of 208% more measurements.

To estimate the correlation between these traits throughout mammalian evolution, we jointly model them with an MBD process on the tree with residual variance. In this analysis, we are primarily interested in the correlation between traits during the MBD process on the tree and estimate trait correlations from the marginal posterior of $\boldsymbol{\Sigma}$. Figure 2 summarizes these findings.

Our results are clearly consistent with the fast-slow trait covariation patterns that pace-of-life theory predicts. The ‘slow’ life history traits (longer gestation, later weaning, larger neonatal body mass, later age at first birth, and longer reproductive lifespan) are all positively correlated with each other and negatively correlated with the two ‘fast’ life history traits (greater litter size and more frequent litters). All correlations are significant (determined by $<5\%$ posterior tail probability) with the notable exception of that between litter size and litter frequency. This apparent lack of correlation may be due to the opposing effects of their joint positive correlation with body mass combined with a trade-off between these two traits that life history theory predicts (Stearns 1989). Nevertheless, our results demonstrate that larger animals tend to have slower life history traits, confirming known patterns and reflecting the central role of body size in life history evolution.

We compare the computational efficiency of our method against that of the sampling method using the MBD model both with and without residual variance. Table 1 shows an increase in overall computational efficiency of two orders-of-magnitude as indicated by the change in ESS per hour. Additionally, we see that our method succeeds at reducing both the amount of computational work per MCMC sample (as indicated by the increase in samples per hour) and autocorrelation (as indicated by the increase in ESS per sample).

6.2   Prokaryote evolution

Comparative genomics has greatly assisted in the formulation of prokaryote evolutionary theories. Several such theories have been inspired by and tested through measuring correlation among different phenotypic and genomic traits. For example, the thermal adaptation hypothesis posits that higher GC content is involved in adaptation to high temperatures because it may offer thermostability to genetic material (Bernardi and Bernardi 1986). The genome streamlining hypothesis attempts to explain the compactness of prokaryotic genomes through natural selection favoring small genomes (Doolittle and Sapienza 1980; Orgel and Crick 1980; Giovannoni et al. 2014). Sabath et al. (2013) argue that lower cell volume is an adaptive response to high temperature. The field is well-aware of the need to account for phylogenetic relationships when measuring correlation, but statistical analyses generally rely on fixed, poorly resolved trees and simple models of trait evolution.

Here, we estimate correlation among a set of genotypic and phenotypic traits while simultaneously accounting for phylogenetic uncertainty and accommodating complexity in the trait evolutionary process. We construct our data set from a study by Goberna and Verdú (2016), who collated cell diameter, cell length, optimum temperature and pH measurements for a large set of prokaryotes. Prior experience in resolving large, unknown trees suggests that we limit our analysis to less than $\sim$750 taxa. As such, we include all taxa with three or more measurements and a selection of the taxa with only two measurements in our analysis. For our selection of 705 taxa, we obtain data on genome length, the number of coding sequences, and GC content from the prokaryotes table in NCBI Genome. Table 2 presents the number of measurements for each trait. We log-transform and standardize all traits (except for GC content which we logit-transform and standardize). To infer the phylogeny, we obtain matching 16S sequences via the ARB software package (Ludwig et al. 2004) that we then align using the SINA Alignment Service (Pruesse et al. 2012) and manually edit.

Through MCMC simulation, we simultaneously infer the sequence and trait evolutionary process. We model the sequence evolutionary process using a general time-reversible model (Tavaré 1986) with gamma-distributed rate variation among sites (Yang 1994). We use an uncorrelated lognormal relaxed clock to model rate variation among branches (Drummond et al. 2006) and specify a Yule birth prior process on the unknown tree (Gernhard 2008). For the trait evolutionary process, we assume an MBD model with residual variance.

Figure 3 displays our estimated maximum clade credibility phylogeny with associated trait measurements, and Figure 4 presents the phylogenetic correlation between those traits. One notable result is the positive correlation between optimal temperature and GC content ($0.22$ posterior mean, $[0.08,0.37]$ 95% highest posterior density interval) that the thermal adaptation hypothesis predicts (Bernardi and Bernardi 1986). Researchers have discussed this hypothesis for years with mixed support (Hurst and Merchant 2001; Musto et al. 2004; Wang et al. 2006; Wu et al. 2012; Sabath et al. 2013; Aptekmann and Nadra 2018). Our analysis, however, includes 435 taxa with measurements for both GC content and optimal growth temperature, making it the largest study we are aware of that accounts for phylogenetic relationships. Interestingly, while cell diameter and cell length are not significantly correlated, they are both positively correlated with genome length. Smaller cells have been associated with smaller genomes in both prokaryotes and unicellular eukaryotes (Shuter et al. 1983; Lynch 2007), but the reasons for this are not fully understood (Dill et al. 2011). We also estimate a relatively strong negative correlation between genome length and optimal temperature ($-0.52$ $[-0.67,-0.37]$), supporting the genomic streamlining hypothesis during thermal adaptation. Note that we do not compare computation times here, as simultaneous inference of the phylogenetic tree makes the sampling method prohibitively slow.

6.3   HIV-1 virulence

Recent years have witnessed a strong interest in using phylogenetic comparative methods to study the heritability of HIV-1 virulence. Initially, Alizon et al. (2010) employed Pagel’s $\lambda$ (Pagel 1999) to measure the extent to which HIV-1 set-point viral load reflects viral shared evolutionary history in the Swiss HIV Cohort Study (Swiss HIV Cohort Study et al. 2009) patients. A relatively high heritability estimate of set-point viral load, a predictive measure of clinical outcome, motivated others to examine to what extent the viral genotype can control for this trait (e.g. Hodcroft et al. 2014; Vrancken et al. 2015). These efforts have resulted in widely varying estimates, from 6% to 59%, prompting a discussion on the methods used to estimate the heritability of virulence (see Mitov and Stadler 2018; Bertels et al. 2018). Here, we revisit the most comprehensive data set recently analyzed (Blanquart et al. 2017) to determine the extent to which variability in HIV-1 virulence is attributable to viral genetic variation. We focus on the dataset of subtype B viruses from Blanquart et al. (2017) that encompasses 1581 taxa with associated measures of set-point viral load and CD4 cell count decline. We rely on the maximum likelihood phylogeny inferred for this data set, but convert it to a time-measured tree with dated tips using a heuristic procedure (To et al. 2016). A prior examination of the correlation between sampling time and root-to-tip divergence using TempEst (Rambaut et al. 2016) indicated the presence of outliers, most of which can be attributed to a basal lineage in the phylogeny. As the subtyping of the taxa in this basal lineage also was ambiguous (Blanquart, personal communication), we remove this lineage (36 taxa) together with 9 other outlier taxa. We note that this resulted in a time to the most recent common ancestor (TMRCA) estimate of about 1960 that is much more in line with a recent subtype B TMRCA estimate (1967, 95% Bayesian credibility interval of 1963–1970; Worobey et al. 2016) than the estimate including the basal lineage ($\sim$1930).

Two measures of set-point viral load are available for all remaining taxa: (i) one based on a standardized choice of assay on a single sample taken between 6 and 24 months after infection and before the initiation of antiretroviral therapy (“gold standard viral load”, GSVL) and (ii) a more classical measure of set-point viral load (SPVL) based on the mean of all log viral loads measured between 6 and 24 months after infection. Figure 5 presents the phylogeny and associated trait values. To estimate heritability of both set-point viral load measures and CD4 slope, we model all three measurements as a multivariate trait in our MBD model with residual variance and approximate the posterior distribution of the heritability statistic $\mathbf{H}$ via MCMC. Our estimated heritabilities are $0.21$ $[0.11,0.3]$ for GSVL, $0.18$ $[0.1,0.26]$ for SPVL, and $0.16$ $[0.07,0.25]$ for CD4 cell decline. These estimates are consistent with similar estimates reported by Blanquart et al. (2017)

We further asses model fit by assessing predictive performance of GSVL on SPVL. We omit CD4 slope from our analysis as it is measured concurrently with SPVL. We randomly remove 5% of the SPVL measurements from the data set and consider four different models. We consider both a bivariate case where we assume a multivariate process and a univariate case where we analyze SPVL alone. For both the bivariate and univariate cases, we use the MBD model both with and without the residual variance extension. For each removed SPVL measurement, we compute the mean squared error (MSE) between the predicted and true values. We repeat each analysis 50 times and report the cumulative results in Figure 6, from which two results emerge. First, the MSE of prediction in the bivariate cases are lower than those in the univariate cases. This is unsurprising given the strong correlation between SPVL and GSVL. Second, addition of residual variance to the model results in modestly better prediction of SPVL in both the bivariate and univariate cases. This emphasizes the importance of including model extensions like residual variance in these analyses.

We again demonstrate improvements in computational efficiency (see Table 1). While less dramatic than the mammals example, we still see an order-of-magnitude increase in effective sample size per hour in the MBD model without residual variance. This attenuation is to be expected, as there are far fewer missing measurements in the HIV data set than the mammal data set. Nevertheless, our method still outperforms the sampling method in the simple MBD model even when only 9.4% of data points are missing. For the model with residual variance, our method outperforms the sampling method by two orders-of-magnitude.

.   iscussion

Oftentimes comparative biologists are interested in phylogenetically adjusted methods for assessing relationships between traits of organisms. However, frequently when the number of taxa grows large the level of missing data increases, making inference challenging. Here, we have developed a method for evaluating the likelihood of observed traits given a tree while integrating out missing values analytically that dramatically outperforms current best-practice methods. In the mammalian data set, with $N=3690$ and $56.8\%$ missing data, we achieve a minimum effective sample size per hour $410\times$ greater than previous methods. This increase in speed brings computation times down from more than a week to less than an hour. Even in the more tractable HIV data set, with $N=1536$ and $9.4\%$ missing data, we increase the minimum ESS per hour by a factor of $65$. Both increases in speed are due to an overall decrease in both autocorrelation between MCMC samples and the amount of computational work required per sample. Importantly, this increase in computational efficiency allows for previously intractable analyses on large trees. Specifically, we incorporate residual variance into the model and (in the prokaryotes example) simultaneously infer $\boldsymbol{\Sigma}$, $\boldsymbol{\Gamma}$, and an unknown phylogeny $\mathscr{F}$. Further, the residual variance extension is only one of several potential extensions. Other possible extensions could incorporate data sets with repeated measurements at the tips of the tree and factor analyses (Tolkoff et al. 2017).

An important limitation of our and previous methods is that they assume an ignorable missing data mechanism (i.e. that the data are missing at random and that the prior on any model parameters is independent of the missing data mechanism) (Little and Rubin 1987). While these conditions may hold in some comparative biology examples, possible violations abound (e.g. any data set where values below some detection limit are omitted). One potential solution explicitly includes these thresholds in the model and renormalizes the observed-data likelihood appropriately.

Finally, and perhaps most importantly, we propose our method as a special case solution to the long-standing statistical problem involving multivariate normal distributions with missing data. Specifically, our method applies to any MVN distribution with a three-point structured covariance matrix (see Ho and Ané 2014). Intuitively, this condition arises in covariance matrices generated from processes that are additive on an acyclic graph (see Figure 7). This restriction, however, is not overly limiting and applies to a broad range of normal models including multilevel hierarchical models and matrix-normal distributions such as the one we use here. Additionally, our pre-order data augmentation procedure enables $\mathcal{O}\left(N\right)$ imputation in these highly structured models. While Allen and Tibshirani (2010) and Glanz and Carvalho (2018) have utilized the EM algorithm (Dempster et al. 1977) to efficiently perform maximum likelihood imputation in similar problems, our method could serve as an alternative for approaches that base inference on the observed-data likelihood.

upplementary Material

Detailed Calculations:

Includes SI Sections 1 through 3 (pdf file)

unding

This work was partially supported by the NIH under training programs T32-GM008185 and T32-HG002536 and Grants R01 AI107034 and U19 AI135995; NSERC Discovery under Grant RGPIN-2018-05447 and Launch Supplement DGECR-2018-00181; the Artic Network via the Wellcome Trust under project 206298/Z/17/Z; the KU Leuven Special Research Fund (‘Bijzonder Onderzoeksfonds’) under Grant OT/14/115; the Research Foundation – Flanders (‘Fonds voor Wetenschappelijk Onderzoek – Vlaanderen’) under Grants G066215N, G0D5117N and G0B9317N; the NSF under Grant DMS 1264153; the European Research Council via the European Union’s Horizon 2020 Research and Innovation Programme under Grant no. 725422-ReservoirDOCS; and startup funds from Dalhousie University and the Canada Research Chairs Program.

cknowledgements

The authors thank François Blanquart and Christophe Fraser for their assistance in compiling the HIV-1 data set, and Marta Goberna for advice in compiling the prokaryote data set.

SI 1.   atrix inversion computations

To evaluate the observed data likelihood, we must compute branch-deflated precisions $\mathbf{P}^{\star}_{i}=\left(\mathbf{P}_{i}^{-}+t_{i}\boldsymbol{\delta}_{i}\boldsymbol{\Sigma}\boldsymbol{\delta}_{i}\right)^{-}$ for $i=1,\ldots,2N-2$. We demonstrate below that this matrix exists and is well-defined under the definition of our pseudo-inverse. Using the permutation matrix $\mathbf{C}_{i}$ from Section 2.1.2, we decompose the diffusion variance $\boldsymbol{\Sigma}$ and node precision $\mathbf{P}_{i}$ such that

[TABLE]

for $i=1,\ldots,2N-2$. We use this decomposition to identify that:

[TABLE]

where

[TABLE]

The matrix $\mathbf{T}$ is the sum of a positive-definite matrix and positive-semidefinite matrix and is therefore invertible.

SI 2.   eritability Statistic

We compute the expectation of the empirical variance $\mathbb{E}\left[\mathbf{S}^{2}\left({\mathbf{Y}}\right)\right]$ under the MBD model with residual variance as follows:

[TABLE]

The multivariate normal distribution of $\mbox{vec}\left[\mathbf{Y}\right]$ implies $\text{Cov}\left(Y_{ik},Y_{jl}\right)=\Sigma_{kl}\Upsilon_{ij}+\Gamma^{-1}_{kl}1_{\left\{i=j\right\}}$ where $1_{\left\{i=j\right\}}$ is an indicator function. Using this information in SI Equation 3,

[TABLE]

Note that $\mathbb{E}\left[\mathbf{Y}_{i}\right]=\mathbf{Y}_{2N-1}$ for $i=1\ldots N$, which implies

[TABLE]

As such, our expression for the expected empirical variance reduces to the following:

[TABLE]

SI 3.   lgorithm for Efficiently Computing $\mathbf{H}$

Note that naive computation of $c_{\sigma}=\frac{1}{N}\mbox{tr}\left[\boldsymbol{\Upsilon}\right]+\frac{1}{{N}^{2}}\mathbf{1}_{N}^{t}\boldsymbol{\Upsilon}\mathbf{1}_{N}$ in Equation 24 would require constructing the $N\times N$ matrix $\boldsymbol{\Upsilon}$ and summing over all its elements, which has computation complexity of at least $\mathcal{O}\left(N^{2}\right)$. For cases where $\mathscr{F}$ is random and changes throughout the MCMC simulation, this quantity must be re-computed each time we compute the statistic. To avoid this issue, we implement an algorithm that avoids constructing $\boldsymbol{\Upsilon}$ in its entirety and simply calculates both $\mbox{tr}\left[\boldsymbol{\Upsilon}\right]$ and $\mathbf{1}_{N}^{t}\boldsymbol{\Upsilon}\mathbf{1}_{N}$ in $\mathcal{O}\left(N\right)$ time. The algorithm performs a post-order traversal of the tree where at each internal node $\nu_{\ell}$ we compute $N_{\lfloor\ell\rfloor}$ (the number of tips below $\nu_{\ell}$), $s_{\lfloor\ell\rfloor}$ (the sum of all elements in $\boldsymbol{\Upsilon}_{\lfloor\ell\rfloor}$), and $d_{\lfloor\ell\rfloor}$ (the sum of the diagonal elements in $\boldsymbol{\Upsilon}_{\lfloor\ell\rfloor}$). We define $\boldsymbol{\Upsilon}_{\lfloor\ell\rfloor}$ as the tree variance-covariance matrix constructed from the sub-tree $\mathscr{F}_{\lfloor\ell\rfloor}$ that is simply the tree that contains only the nodes below $\nu_{\ell}$ with node $\nu_{\ell}$ as its root. For internal nodes $\nu_{\ell}$ with child nodes $\nu_{j}$ and $\nu_{k}$, we accumulate

[TABLE]

At the tips, we initialize with $s_{\lfloor i\rfloor}=d_{\lfloor i\rfloor}=0$ and $N_{\lfloor i\rfloor}=1$. At the root, $s_{\lfloor 2N-1\rfloor}=\mathbf{1}_{N}^{t}\boldsymbol{\Upsilon}\mathbf{1}_{N}$ and $d_{\lfloor 2N-1\rfloor}=\mbox{tr}\left[\boldsymbol{\Upsilon}\right]$. This algorithm visits each node in $\mathscr{F}$ exactly once and has run time $\mathcal{O}\left(N\right)$.