Relaxed random walks at scale

TL;DR

This paper introduces a scalable Bayesian method using Hamiltonian Monte Carlo for fitting relaxed random walk models of trait evolution, enabling efficient analysis of large phylogenies with many taxa.

Contribution

The paper presents a novel HMC-based approach that efficiently infers branch-specific rate variation in RRW models, scaling linearly with the number of taxa.

Findings

Achieves over 300-fold speed-up compared to previous methods.

Successfully applied to a phylogeny with 3650 tips.

Enables analysis of trait evolution in large-scale phylogenetic data.

Abstract

Relaxed random walk (RRW) models of trait evolution introduce branch-specific rate multipliers to modulate the variance of a standard Brownian diffusion process along a phylogeny and more accurately model overdispersed biological data. Increased taxonomic sampling challenges inference under RRWs as the number of unknown parameters grows with the number of taxa. To solve this problem, we present a scalable method to efficiently fit RRWs and infer this branch-specific variation in a Bayesian framework. We develop a Hamiltonian Monte Carlo (HMC) sampler to approximate the high-dimensional, correlated posterior that exploits a closed-form evaluation of the gradient of the trait data log-likelihood with respect to all branch-rate multipliers simultaneously. Our gradient calculation achieves computational complexity that scales only linearly with the number of taxa under study. We compare the efficiency of our HMC sampler to the previously standard univariable Metropolis-Hastings approach while studying the spatial emergence of the West Nile virus in North America in the early 2000s. Our method achieves an over 300-fold speed-increase over the univariable approach. Additionally, we demonstrate the scalability of our method by applying the RRW to study the correlation between five mammalian life history traits in a phylogenetic tree with 3650 tips.

Equations

$${\mathbf{Y}}_i\sim \text{MVN}\left({\mathbf{Y}}_{\text{pa}(i)},t_i\mathbf{V}\left({\varphi }_i\right)\right),$$

$$\mathbf{\mathbf{V}}\left(\phi_{i}\right)=\left\{\begin{array}[]{rl}\phi_{i}\boldsymbol{\Sigma}&\text{rate-scalar parameterization},\phi_{i}>0,\\ \frac{1}{\phi_{i}}\boldsymbol{\Sigma}&\text{scale-mixture-of-normals parameterization},\phi_{i}>0,\\ e^{\phi_{i}}\boldsymbol{\Sigma}&\text{unconstrained parameterization},\phi_{i}\in\mathbb{R},\ \text{and}\\ \boldsymbol{\Sigma}&\text{standard Brownian diffusion}.\end{array}\right.$$

$$p({\mathbf{Y}}_{2N-1})=\text{MVN}\left({\mathbf{\nu }}_0,{\kappa }_0^{-1}\mathbf{\Sigma }\right)$$

$$p(\mathbf{\varphi },\mathbf{\Sigma },s,\mathcal{F},\mathbf{\theta }|\mathbf{Y},\mathbf{S})\propto \underset{\text{likelihood}}{\underbrace{p(\mathbf{Y}|\mathbf{\varphi },\mathbf{\Sigma },\mathcal{F})p(\mathbf{S}|\mathcal{F},\mathbf{\theta })}}\underset{\text{priors}}{\underbrace{p(\mathbf{\varphi }|s)p(s)p(\mathbf{\Sigma })p(\mathcal{F},\mathbf{\theta })}},$$

$$p(\mathbf{\varphi }|\mathbf{\Sigma },s,\mathcal{F},\mathbf{\theta },\mathbf{Y},\mathbf{S})$$

$$p(\mathbf{\Sigma }|\mathbf{\varphi },s,\mathcal{F},\mathbf{\theta },\mathbf{Y},\mathbf{S})$$

$$p(s|\mathbf{\varphi },\mathbf{\Sigma },\mathcal{F},\mathbf{\theta },\mathbf{Y},\mathbf{S})$$

$$p(\mathcal{F},\mathbf{\theta }|\mathbf{\varphi },\mathbf{\Sigma },s,\mathbf{Y},\mathbf{S})$$

$$H(\mathbf{\varphi },\mathbf{\rho })=\underset{\text{potential energy}}{\underbrace{-\log p(\mathbf{\varphi }|\mathbf{\Sigma },s,\mathcal{F},\mathbf{\theta },\mathbf{Y},\mathbf{S})}}+\underset{\text{kinetic energy}}{\underbrace{\frac{1}{2}{\mathbf{\rho }}^t\mathbf{M}\mathbf{\rho }}},$$

$$p(\mathbf{\varphi },\mathbf{\rho }|\mathbf{\Sigma },s,\mathcal{F},\mathbf{\theta },\mathbf{Y},\mathbf{S})\propto e^{-H(\mathbf{\varphi },\mathbf{\rho })}.$$

$$\frac{d{\varphi }_i}{dt}$$

$$\frac{d{\rho }_i}{dt}$$

$$\frac{\partial }{\partial {\varphi }_i}\log p(\mathbf{\varphi }|\mathbf{\Sigma },s,\mathcal{F},\mathbf{Y},\mathbf{S})=\frac{\partial }{\partial {\varphi }_i}\log p(\mathbf{Y}|\mathbf{\varphi },\mathbf{\Sigma },\mathcal{F})+\frac{\partial }{\partial {\varphi }_i}\log p(\mathbf{\varphi }|s).$$

$$\frac{\partial }{\partial {\varphi }_i}\left[\log p(\mathbf{Y})\right]=$$

$$=$$

$$\frac{\partial }{\partial {\varphi }_i}\left[p({\mathbf{Y}}_i|{\mathbf{Y}}_{\lceil i\rceil })\right]$$

$$-\text{mbox}tr\left[{\mathbf{Q}}_i{t}_i\frac{\partial }{\partial {\varphi }_i}\left[\mathbf{V}\left({\varphi }_i\right)\right]\right]\}p({\mathbf{Y}}_i|{\mathbf{Y}}_{\lceil i\rceil }),$$

$$\text{d}{\mathbf{Q}}_i$$

$$\text{d}\left|{\mathbf{Q}}_i^{-1}\right|$$

$$\mathbf{F}\left({\mathbf{Y}}_i\right)=\frac{1}{2}\left\{{\left({\mathbf{Y}}_i-{\mathbf{n}}_i\right)}^t{\mathbf{{\rm Y}}}_i\left({\mathbf{Y}}_i-{\mathbf{n}}_i\right)-\text{mbox}tr\left[{\mathbf{\chi }}_i\right]\right\},$$

$$\frac{\partial }{\partial {\varphi }_i}\left[\log p(\mathbf{Y})\right]=$$

$$=$$

$${\mathbf{\mu }}_i=\left[\begin{array}{c}{\mathbf{n}}_i^{\text{u}}-{\left({\mathbf{Q}}_i^{\text{u}}\right)}^{-1}{\mathbf{Q}}_i^{\text{uo}}\left({\mathbf{Y}}_i^{\text{o}}-{\mathbf{n}}_i^{\text{o}}\right)\\ {\mathbf{Y}}_i^{\text{o}}\end{array}\right]$$

$${\mathbf{Z}}_i=\left[\begin{array}{cc}{\left({\mathbf{Q}}_i^{\text{u}}\right)}^{-1}& 0\\ 0& 0\end{array}\right].$$

$$\mathbb{E}\left[\mathbf{F}\left({\mathbf{Y}}_i\right)\mid \mathbf{Y}\right]$$

$$p({\mathbf{Y}}_{\lfloor i\rfloor }|{\mathbf{Y}}_i)\propto \text{MVN}\left({\mathbf{Y}}_i;{\mathbf{m}}_i,{\mathbf{P}}_i\right),$$

$${\mathbf{P}}_k=\{\begin{array}{l}\infty \times {\mathbf{\delta }}_k\text{if k is a tip}\\ ({\mathbf{P}}_i^{\ast }+{\mathbf{P}}_j^{\ast })\text{otherwise,}\end{array}$$

$$\mathbf{P}_{i}^{*}=\Big{(}\mathbf{P}_{i}^{-}+t_{i}\boldsymbol{\delta}_{i}\mathbf{\mathbf{V}}\left(\phi_{i}\right)\boldsymbol{\delta}_{i}\Big{)}^{-}\ \ \ \text{and}\ \ \ \mathbf{P}_{j}^{*}=\Big{(}\mathbf{P}_{j}^{-}+t_{j}\boldsymbol{\delta}_{j}\mathbf{\mathbf{V}}\left(\phi_{j}\right)\boldsymbol{\delta}_{j}\Big{)}^{-},$$

$${\mathbf{P}}_k{\mathbf{m}}_k=({\mathbf{P}}_i^{\ast }{\mathbf{m}}_i+{\mathbf{P}}_j^{\ast }{\mathbf{m}}_j).$$

$$p({\mathbf{Y}}_i|{\mathbf{Y}}_{\lceil i\rceil })=\text{MVN}\left({\mathbf{Y}}_i;{\mathbf{n}}_i,{\mathbf{Q}}_i\right).$$

$$\mathbf{Q}_{i}=\begin{cases}\kappa_{0}\boldsymbol{\Sigma}^{-1}\ \ \ \ \text{if i is root}\\ \Big{(}(\mathbf{Q}_{i}^{*})^{-1}+t_{i}\mathbf{\mathbf{V}}\left(\phi_{i}\right)\Big{)}^{-1}\ \ \ \ \text{otherwise}\end{cases}$$

$${\mathbf{Q}}_i^{\ast }={\mathbf{P}}_j^{\ast }+{\mathbf{Q}}_k.$$

$$\mathbf{n}_{i}=\begin{cases}\boldsymbol{\nu}_{0}\ \ \ \ \text{if i is root}\\ (\mathbf{Q}_{i}^{*})^{-1}\big{(}\mathbf{P}_{j}^{*}\mathbf{m}_{j}+\mathbf{Q}_{k}\mathbf{n}_{k}\big{)}\ \ \ \ \text{otherwise}.\end{cases}$$

Full text

Version dated:

Relaxed Random Walks at Scale

Alexander A. Fisher1, Xiang Ji1, Philippe Lemey2, and Marc A. Suchard1,3,4

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

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

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

4Department 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

Relaxed random walk (RRW) models of trait evolution introduce branch-specific rate multipliers to modulate the variance of a standard Brownian diffusion process along a phylogeny and more accurately model overdispersed biological data. Increased taxonomic sampling challenges inference under RRWs as the number of unknown parameters grows with the number of taxa. To solve this problem, we present a scalable method to efficiently fit RRWs and infer this branch-specific variation in a Bayesian framework. We develop a Hamiltonian Monte Carlo (HMC) sampler to approximate the high-dimensional, correlated posterior that exploits a closed-form evaluation of the gradient of the trait data log-likelihood with respect to all branch-rate multipliers simultaneously. Our gradient calculation achieves computational complexity that scales only linearly with the number of taxa under study. We compare the efficiency of our HMC sampler to the previously standard univariable Metropolis-Hastings approach while studying the spatial emergence of the West Nile virus in North America in the early 2000s. Our method achieves an over 300-fold speed-increase over the univariable approach. Additionally, we demonstrate the scalability of our method by applying the RRW to study the correlation between five mammalian life history traits in a phylogenetic tree with $3650$ tips.

Introduction

Phylogenetic comparative methods are an indispensable tool to study the evolution of biological traits across taxa while controlling for their shared evolutionary history that confounds the inference of trait correlation (Felsenstein, 1985). Modern comparative methods usually entertain continuous, multivariate traits, although extensions to mixed discrete and continuous outcomes are readily available (Ives and Garland Jr, 2009; Cybis et al., 2015). Approaches typically model trait evolution as a Brownian diffusion or “random walk” process that acts conditionally independently along the branches of a known or random phylogeny. Specifically, the observed or unobserved (latent) trait value of a node in a phylogeny arises from a multivariate normal distribution centered on the latent trait value of its ancestral node with variance proportional to the units of time between nodes. A strict Brownian diffusion model, however, is unable to accommodate the overdispersion in trait data that often emerges from real biological processes (Schluter et al., 1997). One such example arises when examining the dispersal rate of measurably evolving viral pathogens (Biek et al., 2007). For example, if birds serve as the viral host, migratory patterns may induce inhomogeneous dispersal rates over time (Pybus et al., 2012). In such cases, a strict Brownian diffusion model fails to capture, and therefore can also fail to predict, the spatial dynamics of an emerging epidemic.

Lemey et al. (2010) relax the strict Brownian diffusion assumption by introducing branch-rate multipliers that scale the variance of the Brownian diffusion process along each branch of the phylogeny. This “relaxed random walk” (RRW) model requires estimating $2N-2$ correlated branch-rate multipliers, where $N$ is the number of taxa in the phylogeny. Lemey et al. (2010) take a Bayesian approach to parameter estimation where they infer the posterior distribution of the branch-rate multipliers via Markov chain Monte Carlo (MCMC) employing a simple univariable Metropolis-Hastings (UMH) proposal distribution (Hastings, 1970). Since the rates remain correlated in the posterior, a random-scan (Liu, 2008) of UMH proposals inefficiently explores branch-rate space. Specifically, univariable samplers force accepted proposals to be very close together to avoid a large number of rejection steps in the Markov chain simulation. This results in high correlation between MCMC samples from the posterior, making point estimates of the branch-rate multipliers unreliable and slow to converge. For example, in the mammalian life history study we explore here, a UMH sampler requires approximately 33 days of actual run time to achieve reasonable posterior estimates of branch-rate multipliers. Despite this present drawback, RRWs find many impactful applications, e.g., in phylodynamics and phylogeography (Bedford et al., 2014; Faria et al., 2014).

To ameliorate the difficulties that high dimensional MCMC sampling presents, we propose adopting a geometry-informed sampling approach using Hamiltonian Monte Carlo (HMC). HMC equates sampling from a probability distribution with simulating the trajectory of a puck sliding across a frictionless surface warped by the shape of the distribution (Neal, 2011). To map from this statistical problem to the physical one, we view the MCMC samples of our branch-rate multipliers as the “position” of the puck and, then, for each positional dimension we introduce an associated momentum variable. In this way, we extend a D-dimensional parameter space to 2D-dimensional phase space (Betancourt, 2017) and traverse the 2D phase space via differentiating the Hamiltonian and using a numerical integration method to offer proposal states for our MCMC chain. The major limitation to HMC is calculating the gradient of the log-posterior with respect to all position parameters simultaneously. Previous approaches for calculating gradients on phylogenies have employed “pruning”-type algorithms (Felsenstein, 1981) that scale quadratically with the number of taxa in the tree (Bryant et al., 2005). Likewise, numerical approaches also scale quadratically.

In this paper, we derive a method to calculate the gradient with computational complexity that scales only linearly with the number of taxa. We implement our method in the BEAST software package (Suchard et al., 2018), a popular tool for the study and reconstruction of rooted, time-measured phylogenies. We demonstrate the speed and accuracy of our linear-order gradient HMC versus previous best practices by examining the spread of the West Nile virus across the Americas in the early 2000s. Finally, we use our technique to apply the RRW model to study the sensitivity of correlation estimates to model misspecification between mammalian adult body mass, litter size, gestation length, weaning age and litter frequency across $3650$ mammals, thereby demonstrating the scalability of our HMC implementation to tackling a previously intractable problem.

Methods

Model and Inference

Consider a known or random phylogeny $\mathscr{F}$ with $N$ sampled tip nodes and $N-1$ internal and root nodes, each with an observed or latent continuous trait value $\mathbf{Y}_{i}\,\in\,\mathbb{R}^{P}$. To traverse the phylogeny $\mathscr{F}$, let node $\text{pa}(i)$ index the parent of node $i$ with branch length $t_{i}$ connecting the two nodes. Then under the RRW model,

[TABLE]

where the $P\times P$ matrix-valued function $\mathbf{\mathbf{V}}\left(\phi_{i}\right)$ characterizes the branch-specific multivariate normal (MVN) increment that defines the diffusion process. We parameterize this function in terms of an unknown $P\times P$ positive-definite matrix $\boldsymbol{\Sigma}$ that describes the covariation between trait dimensions after controlling for shared evolutionary history and an unknown branch-rate multiplier $\phi_{i}$. Typical choices include

[TABLE]

To complete the RRW model specification, we adopt a prior density on the unobserved trait at the parentless root node,

[TABLE]

with prior mean $\boldsymbol{\nu}_{0}$ and sample-size $\kappa_{0}$.

Letting $\boldsymbol{\phi}=(\phi_{1},\ldots,\phi_{2N-2})$ and the observed data $\mathbf{Y}=(\mathbf{Y}_{1},\ldots,\mathbf{Y}_{N})$ at the tips, we are interested in learning about the posterior

[TABLE]

where $s$ is an unknown parameter characterizing our prior on $\boldsymbol{\phi}$ and $\boldsymbol{\theta}$ represents parameters of a molecular sequence substitution model for the evolution of aligned molecular sequence data $\mathbf{S}$. Note that we follow usual convention (Cybis et al., 2015) and assume that $\mathbf{Y}$ and $\mathbf{S}$ are conditionally independent given $\mathscr{F}$. We follow the example of Lemey et al. (2010) and place a log-normal prior distribution on $\boldsymbol{\phi}$ with mean $1$ and standard deviation $s$. We further assume a relatively uninformative exponential prior on $s$ with mean $10$. Additionally, we assign a Wishart conjugate prior with scale matrix $\mathbf{I}_{P}$ and $P$ degrees of freedom to $\boldsymbol{\Sigma}^{-1}$.

We use MCMC integration to approximate this posterior using a random-scan Metropolis-within-Gibbs approach (Levine and Casella, 2006; Liu, 2008). One cycle of this scheme consists of sampling $\boldsymbol{\phi},\boldsymbol{\Sigma},s$ and then $(\mathscr{F},\boldsymbol{\theta})$ via

[TABLE]

where update (5d) is unecessary when $\mathscr{F}$ is fixed, otherwise efficient sampling from the density $p(\mathscr{F},\boldsymbol{\theta}\,|\,\boldsymbol{\phi},\boldsymbol{\Sigma},s,\mathbf{Y},\mathbf{S})$ is well described elsewhere, see for example Suchard et al. (2018). Updates (5b) and (5c) are straightforward due to the conjugate priors chosen in our model. We turn our focus to the remaining component of our scheme, namely sampling from $p(\boldsymbol{\phi}\,|\,\boldsymbol{\Sigma},s,\mathscr{F},\boldsymbol{\theta},\mathbf{Y},\mathbf{S})$.

Hamiltonian Monte Carlo

We wish to sample $\boldsymbol{\phi}$ jointly to avoid potentially high autocorrelation in the resulting MCMC chain. To this end, we propose using HMC and begin with a brief description of how HMC maps sampling from a probability distribution to simulating a physical system. In classical mechanics, the Hamiltonian is the sum of the kinetic and potential energy in a closed system. To build the connection, we introduce auxiliary momentum variable $\boldsymbol{\rho}=(\rho_{1},\ldots,\rho_{2N-2})$ and write our Hamiltonian,

[TABLE]

where the mass matrix $\mathbf{M}$ weights our momentum variables. The canonical distribution from statistical mechanics relates the joint density of state variables $\boldsymbol{\phi}$ and $\boldsymbol{\rho}$ to the energy in a system via the relationship,

[TABLE]

Substituting our Hamiltonian into (7), we observe that $\boldsymbol{\phi}$ and $\boldsymbol{\rho}$ are independent and recognize the marginal density of ${\boldsymbol{\rho}}$ to be MVN. To start the HMC algorithm, we first sample $\boldsymbol{\rho}$ from this marginal density. Then by differentiating $H(\boldsymbol{\phi},\boldsymbol{\rho})$, we generate Hamilton’s equations of motion,

[TABLE]

We can use the resulting vector field in conjunction with a variety of numerical integration techniques to propose new states of $\boldsymbol{\phi}$ for our MCMC chain. Consistent with typical construction (Neal, 2011), we use the leapfrog method for numerical integration, where we follow the trajectory of $\boldsymbol{\rho}$ for a half-step before updating $\boldsymbol{\phi}$. For a full discussion of HMC, see Neal (2011). Importantly, Hamilton’s equations elicit a need to calculate $\nabla_{\boldsymbol{\phi}}\log p(\boldsymbol{\phi}\,|\,\boldsymbol{\Sigma},s,\mathscr{F},\boldsymbol{\theta},\mathbf{Y},\mathbf{S})=\left(\frac{\partial}{\partial\phi_{1}},\dots,\frac{\partial}{\partial\phi_{2N-2}}\right)^{{}^{t}}\log p(\boldsymbol{\phi}\,|\,\boldsymbol{\Sigma},s,\mathscr{F},\boldsymbol{\theta},\mathbf{Y},\mathbf{S})$ to traverse phase space.

Gradient of trait data log-likelihood

A practical HMC sampler demands efficient calculation of $\nabla_{\boldsymbol{\phi}}\log p(\boldsymbol{\phi}\,|\,\boldsymbol{\Sigma},s,\mathscr{F},\boldsymbol{\theta},\mathbf{Y},\mathbf{S})$. Differentiating the logarithm of (5a), we obtain

[TABLE]

Our log-normal prior choice for $\boldsymbol{\phi}$ renders evaluating the second term in Equation (9) trivial. Here we develop a general recursive algorithm for calculating $\nabla_{\boldsymbol{\phi}}\log p(\mathbf{Y}\,|\,\boldsymbol{\phi},\boldsymbol{\Sigma},\mathscr{F})$. To facilitate this development, consider splitting $\mathbf{Y}$ into two disjoint sets relative to any node $i$ in $\mathscr{F}$. We define $\mathbf{Y}_{\lfloor i\rfloor}$ as the observed data descendant of node $i$ and $\mathbf{Y}_{\lceil i\rceil}$ as the observed data “above” (or not descendent of) node $i$. For clarity, see Figure (1).

In the following, we drop the dependence of the log-likelihood on $\boldsymbol{\phi}$, $\boldsymbol{\Sigma}$ and $\mathscr{F}$ for notational convenience. To begin,

[TABLE]

The last equality above follows from the fact that $\phi_{i}$ is associated only with the branch above node $i$. Therefore when we condition on $\mathbf{Y}_{i}$, $\mathbf{Y}_{\lfloor i\rfloor}$ is independent of $\phi_{i}$. Similarly, $\mathbf{Y}_{\lceil i\rceil}$ evolves independent of $\phi_{i}$. To proceed with the differential above, we use the fact that $p(\mathbf{Y}_{i}\,|\,\mathbf{Y}_{\lceil i\rceil})$ follows a MVN distribution with as of yet undetermined mean $\mathbf{n}_{i}$ and precision $\mathbf{Q}_{i}$ (see Appendix (Pre-Order Partial Likelihood) for a detailed derivation). We extract the middle term from Equation (10) and find

[TABLE]

using the differential properties

[TABLE]

found in, e.g., Petersen et al. (2008).

To simplify notation, we let function

[TABLE]

where $\boldsymbol{\Upsilon}_{i}=\mathbf{Q}_{i}\hskip 1.00006ptt_{i}\frac{\partial}{\partial\phi_{i}}\left[\mathbf{\mathbf{V}}\left(\phi_{i}\right)\right]\mathbf{Q}_{i}$ and $\boldsymbol{\chi}_{i}=\mathbf{Q}_{i}\hskip 1.00006ptt_{i}\frac{\partial}{\partial\phi_{i}}\left[\mathbf{\mathbf{V}}\left(\phi_{i}\right)\right]$. Substituting Equation (13) back into Equation (10), we observe that

[TABLE]

When $\mathbf{Y}_{i}$ is fully observed (typically $i\leq N$), this expectation collapses to the direct evaluation of $\mathbf{F}\left(\mathbf{Y}_{i}\right)$. When $i=N+1,\ldots,2N-2$ or if $\mathbf{Y}_{i}$ is partially observed for $i=1,\ldots,N$, we require $p(\mathbf{Y}_{i}\,|\,\mathbf{Y})$. From Bayes’ theorem, $p(\mathbf{Y}_{i}\,|\,\mathbf{Y})\propto p(\mathbf{Y}_{\lfloor i\rfloor}\,|\,\mathbf{Y}_{i})p(\mathbf{Y}_{i}\,|\,\mathbf{Y}_{\lceil i\rceil})$. Partial likelihood $p(\mathbf{Y}_{\lfloor i\rfloor}\,|\,\mathbf{Y}_{i})$ is proportional to a MVN density characterized by computable mean $\mathbf{m}_{i}$ and precision $\mathbf{P}_{i}$ (Pybus et al., 2012). Using this fact, $p(\mathbf{Y}_{i}\,|\,\mathbf{Y})$ becomes MVN with mean $\boldsymbol{\mu}_{i}=\mathbf{Z}_{i}\left(\mathbf{P}_{i}\mathbf{m}_{i}+\mathbf{Q}_{i}\mathbf{n}_{i}\right)$ and variance $\mathbf{Z}_{i}=\left[\mathbf{P}_{i}+\mathbf{Q}_{i}\right]^{-1}.$ When tip $i$ is partially observed, we partition $\mathbf{Y}_{i}=\left(\mathbf{Y}_{i}^{\text{u}},\ \mathbf{Y}_{i}^{\text{o}}\right)^{t}$ into its unobserved and observed entries. Using properties of the conditional MVN, $p(\mathbf{Y}_{i}\,|\,\mathbf{Y})$ becomes degenerate with mean

[TABLE]

and variance

[TABLE]

Finally, for both partially and completely unobserved cases above,

[TABLE]

Equation (17) provides a recipe to compute $\nabla_{\boldsymbol{\phi}}\log p(\mathbf{Y}\,|\,\boldsymbol{\phi},\boldsymbol{\Sigma},\mathscr{F})$ using the means and precisions that characterize partial data likelihoods $p(\mathbf{Y}_{i}\,|\,\mathbf{Y}_{\lceil i\rceil})$ and $p(\mathbf{Y}_{\lfloor i\rfloor}\,|\,\mathbf{Y}_{i})$.

Tree Traversals

We introduce post- and pre-order tree traversals to recursively calculate all partial data likelihood means and precisions in computational complexity ${\cal O}\left(N\right)$ that scales linearly with $N$. To begin, let nodes $i$ and $j$ be daughters of node $k$. Following Hassler et al. (2019), let $\boldsymbol{\delta}_{k}=\text{diag}(\delta_{k1},\ldots\delta_{kP})$ for $k=1,\ldots,N$ be a diagonal matrix with indicator elements $\delta_{kp}$ that take value $1$ if $\text{Y}_{kp}$ is observed and [math] if not. For the post-order traversal,

[TABLE]

with post-order mean $\mathbf{m}_{i}$ and precision $\mathbf{P}_{i}$. For $k=1,\ldots,2N-1$ in post-order, we build the precision via

[TABLE]

with the definition that $\infty\times 0=0$ and

[TABLE]

where the pseudo-inverse, defined and developed by Bastide et al. (2018); Hassler et al. (2019), is described in Appendix (Pseudo-inverse). At the tips, we build the mean $\mathbf{m}_{k}=\boldsymbol{\delta}_{k}\odot\mathbf{Y}_{k}$ where $\odot$ is the elementwise dot product, and for the internal nodes, $\mathbf{m}_{k}$ is a solution to

[TABLE]

For a proof of these post-order updates, see Hassler et al. (2019) (Supplemental Material).

To compute $p(\mathbf{Y}_{i}\,|\,\mathbf{Y}_{\lceil i\rceil})$, we traverse the tree in pre-order fashion according to our generalized version of the recursive algorithm proposed by Cybis et al. (2015). See Appendix (Pre-Order Partial Likelihood) for a derivation of our generalized pre-order update. For the pre-order traversal,

[TABLE]

For $i=2N-1,\ldots,1$ looking down the tree, we update our pre-order precision,

[TABLE]

at each node where

[TABLE]

We also keep track of the pre-order mean at each node via

[TABLE]

Both traversals visit each node exactly once and perform a matrix inversion as their most costly operation, providing an ${\cal O}\left(NP^{3}\right)$ algorithm. However, as we observe in Equation (2), generally $\mathbf{\mathbf{V}}\left(\phi_{i}\right)=g{}\left(\phi_{i}\right)\boldsymbol{\Sigma}$. In this case, we can further reduce the computational complexity to ${\cal O}\left(NP^{2}\right)$ by factoring out $\boldsymbol{\Sigma}$. Instead of inverting $\mathbf{\mathbf{V}}\left(\phi_{i}\right)$ at each step, we only need to invert $\boldsymbol{\Sigma}$ at most once per likelihood or gradient evaluation.

Results

West Nile Virus

West Nile virus (WNV) is responsible for more than 1,500 deaths and caused over 700,000 illnesses since first reported in North America in 1999. The virus typically spreads via mosquito bites; however, the primary host is birds. First identified in New York City, WNV spread to the Pacific coast by 2003 and reached south into Argentina by 2005 (Petersen et al., 2013). We examine whole aligned viral genomes (11,029 nt) and geographic data on $104$ cases of WNV collected between 1999 and 2007 (Pybus et al., 2012). In cases where only the year of sampling is known, we set the sampling date to the midpoint of that year. Previous authors have recorded latitude and longitude geographic sampling information by converting zip code locations using ZIPList5. For 27 of the specimens, only the U.S. or Mexican state of discovery is known and so we have augmented sampling data with the coordinates of the centroid of the state (Pybus et al., 2012).

Here we study the simultaneous evolution and dispersal of WNV as it spreads across North America, following the modeling choices of Pybus et al. (2012). We define geographic location as our trait of interest $\mathbf{Y}$ within a RRW and infer rates $\boldsymbol{\phi}$ using our new HMC method. In two separate inference scenarios, we compare the computational efficiency of our method to the random-scan UMH approach employed by Pybus et al. (2012).

To begin, we set up a RRW model with log-normal prior on rates $\boldsymbol{\phi}$ with mean $=1$ and standard deviation $s$ and use a general time-reversible (GTR) + $\Gamma$ substitution model with a log-normal relaxed molecular clock. We use the UMH transition kernel to run a 250 million state MCMC chain simulation to obtain posterior mean estimates of $\boldsymbol{\Sigma}$ and $s$ as well as a maximum clade crediblility (MCC) tree. In scenario (a) we use these fixed model parameters and topology to strictly sample $\boldsymbol{\phi}$ using both the HMC and UMH transition kernels. Under this fixed analysis we run our HMC-based chain for 1 million states and a UMH-based chain for 150 million states. We use effective sample size (ESS) of the posterior $\phi_{i}$ samples for all $i$ divided by computational runtime to evaluate the performance of each MCMC approach and report densities of ESS/second across all branches in Figure (2). ESS/second is averaged across five runs each with uniform (0-10) random initial branch-rate multipliers. The median ESS/second across $\boldsymbol{\phi}$ is $118$ and $0.377$ for the HMC and UMH transition kernels respectively. This demonstrates a $312$-fold speed increase. Additionally, the minimum ESS/second is $12.4$ with HMC and $0.183$ with UMH, exhibiting a $67.7$-fold speed-up for the “least well” explored $\phi_{i}$.

In scenario (b), we use a random starting tree and jointly estimate all parameters ($\boldsymbol{\phi},\boldsymbol{\Sigma}$, $s$, $\mathscr{F}$ and $\boldsymbol{\theta}$) of the full posterior (4). Since branch-specific $\phi_{i}$ are no longer identifiable when $\mathscr{F}$ is random, we compare square jump distance across all $\boldsymbol{\phi}$ between samples from the posterior under both MCMC regimes to compare efficiency. We run HMC chains for 22.5 million states so that we are sampling from the posterior distribution of all parameters and we save the state of BEAST. Subsequently, we run both HMC and UMH chains from the same saved states and compute lag-7 square jump distance to adjust for the relative weight of the transition kernel in the full analysis. Since the UMH sampler updates only one branch-rate multiplier at a time, we compare square jump distance between samples of our HMC chain with samples from the UMH chain that are lagged $(2N-2)\times$ (i.e. $206\times$) farther apart. We run each MCMC simulation until we obtain 5000 samples from the posterior and report the average median across five separate runs. In this comparison, we find that the average square jump distances from five separate runs is $1457$ and $128.2$ for the HMC and UMH chains respectively.

In Figure (3) we report the MCC tree, obtained from applying HMC to the RRW model as described in scenario (b), where substitution rate variation is accounted for by the molecular clock model. The branch with the highest posterior dispersal rate starts the WN02 lineage identified by Gray et al. (2010). The clade of New York isolates sampled in 1999, however, maintains a much slower dispersal rate. Sufficiently precise posterior estimates under the previous UMH approach remain unavailable in reasonable time to expose these differences.

Mammalian Life History

Life history theory aims to explain how traits such as adult body mass, litter size and lifespan evolve to optimize reproductive success (Stearns, 2000). Life history theory finds important applications in determining a species’ fecundity and predicting extinction risk in response to changing environmental stimuli (Pacifici et al., 2017; Fritz et al., 2009; de Silva and Leimgruber, 2019), but due to the sparseness of much life history data, it is essential to understand how traits covary to make meaningful predictions (Santini et al., 2016). To determine which traits covary, comparative mammalian life-history studies posit a ‘fast–slow’ continuum, claiming small mammals are typically ‘fast’, characterized by early maturation, large litters and shorter lifespans, while larger mammals are typically ‘slow’ and present contrasting characteristics (Oli, 2004; Millar and Zammuto, 1983). Under this framework, certain traits such as gestation length, weaning age and body mass are predicted to be positively correlated, but reported estimates of positive correlation from data may be artifacts of the restrictive assumptions of strict Brownian diffusion modeling. Here we re-evaluate this claim within the RRW of trait evolution made tractable through ${\cal O}\left(N\right)$ HMC sampling.

Under the RRW, we infer correlation between five life history traits from the PanTHERIA data set (Jones et al., 2009), namely body mass, litter size, gestation length, weaning age and litter frequency across $3650$ mammalian species related by the fixed supertree of Fritz et al. (2009). To obtain this subset of the supertree, we only consider taxa for which at least one of these five traits is observed. We take the intersection of this set of taxa with those in the fixed supertree of Fritz et al. (2009) and prune all other observations from the tree. We log-transform and standardize the trait measurements and subsequently estimate posterior mean correlations between each pair of traits under the RRW using an HMC-based chain for 300 thousand states. To gauge the effect of a heterogeneous diffusion process on the correlation between traits, we also make inference using the strict Brownian diffusion model where the $\boldsymbol{\phi}$ are all identically $1$. Here we perform MCMC inference on the diffusion matrix $\boldsymbol{\Sigma}$ for 100 thousand states. We report posterior mean estimates of correlation between each pair of traits under both the RRW and strict Brownian diffusion in Figure (4).

In most cases, the RRW reassuringly confirms analysis under the more limited model. However, in some cases our confidence in the sign of the correlation differs between models and in one instance the sign of the posterior mean correlation disagrees. Under the RRW, we observe positive posterior mean correlation of $0.021$ between litter frequency and litter size with posterior odds ratio $2.63$ that the correlation is positive. Under the strict Brownian diffusion model we observe a negative posterior mean correlation of $-0.018$ with posterior odds ratio $2.25$ of being negative, indicating slightly weaker belief in the correlation’s sign under the strict model.

Discussion

Previous MCMC techniques to investigate trait evolution under the RRW model scale poorly with large data sets. Specifically the UMH transition kernel is ineffective for sampling correlated, high dimensional parameter space. We provide a remedy by using an HMC transition kernel to sample all branch-rate multipliers simultaneously. To improve the speed of HMC we derive an algorithm for calculating the gradient of the trait data log-likelihood. This gradient calculation achieves ${\cal O}\left(N\right)$ computational speed, a vast improvement compared to both numerical and pruning methods for calculating the gradient that typically require ${\cal O}\left(N^{2}\right)$.

We observe over 300-fold speed-up when comparing our HMC transition kernel to UMH in the spread of the WNV across North America in the early 2000s. The resulting MCC tree reveals that the largest dispersal rate precedes the most recent common ancestor of the WN02 lineage. Subsequently, the dispersal rates slow down through the WN02 clade. This suggests that this clade developed after some rapid geographic displacement. Interestingly, the appearance of smaller branch-rate multipliers within the WN02 lineage is consistent with the slowing speed of sequence evolution as described in Snapinn et al. (2007).

As exhibited in Figure (2), ESS from posterior sampling accumulates at variable speed across the branches of the tree. To further improve the sampling of our HMC algorithm, one might use an approximation of the posterior covariance of $\boldsymbol{\phi}$ for the mass matrix $\mathbf{M}$ to appropriately weight momentum updates in the HMC algorithm (Neal, 2011). Possible approximations include the Hessian of the log-posterior (a local approximation of the curvature of branch-rate multiplier space) or the sample variance across each dimension. An important consideration in choosing an appropriate $\mathbf{M}$ is whether one is studying under a fixed or random phylogeny $\mathscr{F}$. Since varying $\mathscr{F}$ in the posterior often creates multimodal distributions of $\boldsymbol{\phi}$, local approximations such as the Hessian may be of limited assistance in such cases.

We show in our application to mammalian life history that our computationally efficient HMC algorithm imbues the RRW model with the ability to handle large trees with thousands of taxa. By applying the RRW model to this massive example we confirm that large mammals have ‘slower’ life history characteristics, exhibited by the positive correlation between body mass, gestation length and weaning age, while smaller mammals scale in the opposite manner and tend to have high litter frequency and size, see Figure (4). The posterior mean correlation between litter frequency and litter size changes sign under each model, but with low posterior probability reflecting a lack of correlation between these traits. Importantly, our method allows us to obtain results in 48 hours where the previous UMH method would take approximately 33 days.

In a time where biological data are more prolific than ever, scalable approaches to complex models of evolution such as the RRW prove increasingly useful in a variety of applications. From spatial epidemiology where determining the dispersal rate of an infectious disease is crucial, to evolutionary ecology where understanding life-history can provide insight into declining animal populations, the analysis of data is becoming a bottleneck to the scientific process and the need for computationally faster approaches stands evident. We hope that this work will serve to improve the speed of such analyses.

Acknowledgments

This work was supported by the European Research Council under the European Union’s Horizon 2020 research and innovation programme (grant agreement number 725422 - ReservoirDOCS); The Artic Network from the Wellcome Trust (grant number 206298/Z/17/Z); National Science Foundation (grant number DMS 1264153); National Institutes of Health (grants R01 AI107034, U19 AI135995, and T32 GM008185 to AAF); and Research Foundation – Flanders (‘Fonds voor Wetenschappelijk Onderzoek – Vlaanderen’, G066215N, G0D5117N and G0B9317N) to PL.

Appendices

We provide public access to all BEAST XML files used in this work at http://github.com/suchard-group/RRW_at_scale.

Pre-Order Partial Likelihood

Here we derive a generalized version of the pre-order recursive algorithm proposed by Cybis et al. (2015) to compute $p(\mathbf{Y}_{i}\,|\,\mathbf{Y}_{\lceil i\rceil})$ for all $i$ in ${\cal O}\left(N\right)$. We begin with the law of total probability,

[TABLE]

for node $i$ with parent $k$ and sibling $j$. Recalling that

[TABLE]

we identify Equation (26) as a recursive expression whose solution has the form

[TABLE]

with presently undetermined pre-order mean $\mathbf{n}_{i}$ and pre-order precision $\mathbf{Q}_{i}$.

We unravel these quantities by first identifying that $p(\mathbf{Y}_{2N-1}\,|\,\mathbf{Y}_{\lceil 2N-1\rceil})=p(\mathbf{Y}_{2N-1})$ and set $\mathbf{n}_{2N-1}=\boldsymbol{\nu}_{0}$ and $\mathbf{Q}_{2N-1}=\kappa_{0}\boldsymbol{\Sigma}^{-1}$. Then proceeding in pre-order fashion for $i=2N-2,\ldots,1$

[TABLE]

Pseudo-inverse

The pseudo-inverse used in the post-order tree traversal and defined by Bastide et al. (2018) and Hassler et al. (2019) is an operation for inverting precision and variance matrices with diagonal entries that take the value $\infty$. To invert a diagonal precision matrix, $\mathbf{P}_{i}$ with entries $\infty$ and [math], we define $\infty^{-}=0$ and $0^{-}=\infty$. To invert the variance matrix, $\left(\mathbf{P}_{i}^{-}+t_{i}\boldsymbol{\delta}_{i}\mathbf{\mathbf{V}}\left(\phi_{i}\right)\boldsymbol{\delta}_{i}\right)$ we invert the block matrix of observed trait covariation and invert the remaining diagonal elements using the convention that $\infty^{-}=0$.