Efficient recycled algorithms for quantitative trait models on phylogenies

Phylogenetic models for continuous trait evolution, like those for discrete traits, are specified by the density of trait values at the root and the transition densities along the branches. We use $f(x_{r}|\theta_{r})$ to denote the density for the trait value at the root, where $\theta_{r}$ is a set of relevant model parameters. We use $f(x_{i}|x_{j},\theta_{i})$ to denote the transitional density for the value at node $i$, conditional on the trait value at its parent node $j$. Here, $\theta_{i}$ represents a bundle of parameters related to node $i$ such as branch length, population size, and mutation rate. All of these parameters could vary throughout the tree.

To see how the model works, consider how continuous traits might be simulated. A state $X_{r}$ is sampled from the root density $f(X_{r}|\theta_{r})$. We now proceed through the phylogeny from the root to the tips, each time visiting a node only after its parent has already been visited. For each node $i$, we generate the value at that node from the density $f(X_{i}|x_{j},\theta_{v})$, where $x_{j}$ is the simulated trait value at node $j$, the parent of node $i$. In this way, we will eventually generate trait values for the tips.

We use $X_{1},\ldots,X_{n}$ to denote the random trait values at the tips and $X_{n+1},\ldots,X_{2n-1}$ to denote the random trait values at the internal nodes, ordered so that children come before parents. Hence $X_{2n-1}$ is the state assigned to the root. Let

$$\mathcal{E}(T)=\{(i,j):\mbox{node $i$ is a child node $j$}\}$$

(2)

denote the set of branches in the tree. The joint density for all trait values, observed and ancestral, is given by multiplying the root density with all of the transition densities

$$f(x_{1},\ldots,x_{n},x_{n+1},\ldots,x_{2n-1}|\theta)=f(x_{2n-1}|\theta)\!\!\!\prod_{(i,j)\in\mathcal{E}(T)}\!\!\!f(x_{i}|x_{j},\theta_{i}).$$

(3)

The probability of the observed trait values $x_{1},\ldots,x_{n}$ is now determined by integrating out all of the ancestral trait values:

$$\mathcal{L}(T)=f(x_{1},\ldots,x_{n}|\theta)=\int\int\cdots\int f(x_{2n-1}|\theta_{r})\!\!\!\prod_{(i,j)\in\mathcal{E}(T)}\!\!\!f(x_{i}|x_{j},\theta_{i})\,\,\mathrm{d}x_{n+1},\ldots,\,\mathrm{d}x_{2n-1}.$$

(4)

In these integrals, the bounds of integration will vary according to the model.

The oldest, and most widely used, continuous trait models assume that traits (or transformed gene frequencies) evolve like Brownian motion (Cavalli-Sforza and Edwards, 1967, Felsenstein, 1973). For these models, the root density $f(x_{r}|\theta)$ is Gaussian (normal) with mean $0$ and unknown variance $\sigma_{r}^{2}$. The transition densities $f(x_{i}|x_{j},\theta_{v})$ are also Gaussian, with mean $x_{j}$ (the trait value of the parent) and standard deviation proportional to branch length. Note that there are identifiability issues which arise with the inference of the root position under this model, necessitating a few tweaks in practice.

It can be shown that when the root density and transitional densities are all Gaussian, the joint density (4) is multivariate Gaussian. Furthermore, the covariance matrix for this density has a special structure which methods such as the pruning technique of Felsenstein (1973) exploit. This technique continues to work when Brownian motion is replaced by an OU process (Felsenstein, 1988, Hansen, 1997, Lande, 1976).

Landis et al. (2013) discuss a class of continuous trait models which are based on Lévy processes and include jumps. At particular times, as governed by a Poisson process, the trait value jumps to a value drawn from a given density. Examples include a compound Poisson process with Gaussian jumps and a variance Gamma model with Gamma distributed jumps. Both of these processes have analytical transition probabilities in some special cases.

Lepage et al. (2006) use the Cox-Ingersoll-Ross (CIR) process to model rate variation across a phylogeny. Like the OU process (but unlike Brownian motion), the CIR process is ergodic. It has a stationary Gamma density which can be used for the root density. The transition density is a particular non-central Chi-squared density and the process only assumes positive values.

Kutsukake and Innan (2013) examine a family of compound Poisson models, focusing particularly on a model where the trait values make exponentially distributed jumps upwards or downwards. In the case that the rates of upward and downward jumps are the same, the model reduces to the variance Gamma model of Landis et al. (2013) and has an analytical probability density function.

Sirén et al. (2011) propose a simple and elegant model for gene frequencies whereby the root value is drawn from a Beta distribution and each transitional density is Beta with appropriately chosen parameters.

Trait values at the tips are not always observed directly. A simple, but important, example of this is the threshold model of Wright (1934), explored by Felsenstein (2005). Under this model, the trait value itself is censored and we only observe whether or not the value is positive or negative. A similar complication arises when dealing with gene frequency data as we typically do not observe the actual gene frequency but instead a binomially distributed sample based on that frequency (Sirén et al., 2011).

If the trait values at the tip are not directly observed we integrate over these values as well. Let $\pi(z_{i}|x_{i},\theta_{i})$ denote the probability of observing $z_{i}$ given the trait value $x_{i}$. The marginalised likelihood is then

$$\mathcal{L}(T|z_{1},\ldots,z_{n})=\int\int\cdots\int f(x_{r}|\theta)\!\!\!\prod_{(i,j)\in\mathcal{E}(T)}\!\!\!f(x_{i}|x_{j},\theta_{v})\!\prod_{i=1}^{n}\pi(z_{i}|x_{i},\theta_{i})\,\,\mathrm{d}x_{1},\ldots,\,\mathrm{d}x_{2n-1}.$$

(5)

Analytical integration can be difficult or impossible. For the most part, it is unusual for an integral to have an analytical solution, and there is no general method for finding it when it does exist. In contrast, numerical integration techniques (also known as numerical quadrature) are remarkably effective and are often easy to implement. A numerical integration method computes an approximation of the integral from function values at a finite number of points. Hence we can obtain approximate integrals of functions even when we don’t have an equation for the function itself. See Cheney and Kincaid (2012) for an introduction to numerical integration, and Dahlquist and Björck (2008) and Davis and Rabinowitz (2007) for more comprehensive technical surveys.

The idea behind most numerical integration techniques is to approximate the target function using a function which is easy to integrate. In this paper we will restrict our attention to Simpson’s method which approximates the original function using piecewise quadratic functions. To approximate an integral $\int_{a}^{b}f(x)\,\mathrm{d}x$ we first determine $N+1$ equally spaced points

$$x_{0}=a,\,\,\,\,x_{1}=a+\frac{b-a}{N},\,\,\,\,x_{2}=a+2\frac{b-a}{N},\ldots,x_{k}=a+k\frac{b-a}{N},\ldots,\,\,x_{N}=b.$$

(6)

We now divide the integration into $N/2$ intervals

$$\int_{a}^{b}f(x)\,\,\mathrm{d}x=\sum_{\ell=1}^{N/2}\,\,\,\,\int\displaylimits_{x_{2\ell-2}}^{x_{2\ell}}f(x)\,\,\mathrm{d}x.$$

(7)

Within each interval $[x_{2\ell-2},x_{2\ell}]$, there is a unique quadratic function which equals $f(x)$ at each the three points $x=x_{2\ell-2}$, $x=x_{2\ell-1}$ and $x=x_{2\ell}$. The integral of this quadratic on the interval $[x_{2\ell-2},x_{2\ell}]$ is $\frac{(b-a)}{3N}\left(f(x_{2\ell-2})+4f(x_{2\ell-1})+f(x_{2\ell})\right)$. Summing over $\ell$, we obtain the approximation

$$\int_{a}^{b}f(x)\,\,\mathrm{d}x\approx\sum_{\ell=1}^{N/2}\frac{(b-a)}{3N}\left(f(x_{2\ell-2})+4f(x_{2\ell-1})+f(x_{2\ell})\right).$$

(8)

With a little rearrangement, the approximation can be written in the form

$$\int_{a}^{b}f(x)\,\,\mathrm{d}x\approx\frac{(b-a)}{N}\sum_{k=0}^{N}w_{k}f(x_{k})$$

(9)

where $w_{k}=4/3$ when $k$ is odd and $w_{k}=2/3$ when $k$ is even, with the exception of $w_{0}$ and $w_{N}$ which both equal $1/3$. Simpson’s method is embarrassingly easy to implement and has a convergence rate of $O(N^{-4})$. Increasing the number of intervals by a factor of $10$ decreases the error by a factor of $10^{-4}$. See Dahlquist and Björck (2008) and Davis and Rabinowitz (2007) for further details.

It should be remembered, however, that the convergence rate is still only an asymptotic bound, and gives no guarantees on how well the method performs for a specific function and choice of $N$. Simpson’s method, for example, can perform quite poorly when the function being integrated has rapid changes or soft peaks. We observed this behaviour when implementing threshold models, as described below. Our response was to better tailor the integration method for the functions appearing. We noted that the numerical integrations we carried out all had the form

$$\int_{a}^{b}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}f(x)\,\,\mathrm{d}x$$

(10)

where $\mu$ and $\sigma$ varied. Using the same general approach as Simpson’s rule, we approximated $f(x)$, rather than the whole function $e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}f(x)$, by a piecewise quadratic function $p(x)$. We could then use standard techniques and tools to evaluate $\int_{a}^{b}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}p(x)\,\,\mathrm{d}x$ numerically. The resulting integration formula, which we call the Gaussian kernel method, gives a significant improvement in numerical accuracy.

A further complication is that, in models of continuous traits, the trait value often ranges over the whole real line, or at least over the set of positive reals. Hence, we need to approximate integrals of the form

$$\int_{-\infty}^{\infty}f(x)\,\,\mathrm{d}x\mbox{ or }\int_{0}^{\infty}f(x)\,\,\mathrm{d}x$$

(11)

though the methods discussed above only apply to integrals on finite intervals. We truncate these integrals, determining values $U$ and $L$ such that the difference

$$\int_{-\infty}^{\infty}f(x)\,\,\mathrm{d}x-\int_{L}^{U}f(x)\,\,\mathrm{d}x$$

(12)

between the full integral $\int_{-\infty}^{\infty}f(x)\,\,\mathrm{d}x$ and the truncated integral $\int_{L}^{U}f(x)\,\,\mathrm{d}x$ can be bounded analytically. Other strategies are possible; see Dahlquist and Björck (2008) for a comprehensive review.

Felsenstein has developed pruning algorithms for both continuous and discrete characters (Felsenstein, 1981a, b). His algorithm for continuous characters works only for Gaussian processes. Our approach is to take his algorithm for discrete characters and adapt it to continuous characters.

The (discrete character) pruning algorithm is an application of dynamic programming. For each node $i$, and each state $x$, we compute the probability of observing the states for all tips which are descendants of node $i$, conditional on node $i$ having ancestral state $x$. This probability is called the partial likelihood at node $i$ given state $x$. Our algorithm follows the same scheme, with one major difference. Since traits are continuous, we cannot store all possible partial likelihoods. Instead, we store likelihoods for a finite set of values and plug these values into a numerical integration routine.

Let $i$ be the index of a node in the tree not equal to the root, let node $j$ be its parent node. We define the partial likelihood, $\mathcal{F}_{i}(x)$ to be the likelihood for the observed trait values at the tips which are descendants of node $i$, conditional on the parent node $j$ having trait value $x$. If node $i$ is a tip with observed trait value $x_{i}$ we have

$$\mathcal{F}_{i}(x)=f(x_{i}|x,\theta_{i})$$

(13)

recalling that $f(x_{i}|x,\theta_{i})$ is the density for the value of the trait at node $i$ conditional on the value of the trait for its parent. More generally, we may only observe some value $z_{i}$ for which we have the conditional probability $\pi(z_{i}|x_{i},\theta_{i})$ conditional on the trait value $x_{i}$. In this case the partial likelihood is given by

$$\mathcal{F}_{i}(x)=\int f(\tilde{x}|x,\theta_{i})\pi(z_{i}|\tilde{x})\,\,\mathrm{d}\tilde{x}.$$

(14)

Suppose node $i$ is not the root and that it has two children $u,v$. Since trait evolution is conditionally independent on disjoint subtrees, we obtain the recursive formula

$$\mathcal{F}_{i}(x)=\int f(\tilde{x}|x,\theta_{i})\mathcal{F}_{u}(\tilde{x})\mathcal{F}_{v}(\tilde{x})\,d\tilde{x}.$$

(15)

Finally, suppose that node $i$ is the root and has two children $u,v$. We evaluate the complete tree likelihood using the density of the trait value at the root,

$$\mathcal{L}(T)=\int f(x|\theta_{r})\mathcal{F}_{u}(x)\mathcal{F}_{v}(x)\,\,\mathrm{d}x.$$

(16)

The bounds of integration in (14)—(16) will vary according to the model.

We use numerical integration techniques to approximate (14)—(16) and dynamic programming to avoid an exponential explosion in the computation time. Let $N$ denote the number of function evaluations for each node. In practice, this might vary over the tree, but for simplicity we assume that it is constant. For each node $i$, we select $N+1$ trait values

$$X_{i}[0]<X_{i}[1]<\cdots<X_{i}[N].$$

(17)

How we do this will depend on the trait model and the numerical integration technique. If, for example, the trait values vary between $a$ and $b$ and we are applying Simpson’s method with $N$ intervals we would use $X_{i}[k]=a+\frac{b-a}{N}k$ for $k=0,1,2,\ldots,N$.

We traverse the tree starting at the tips and working towards the root. For each non-root node $i$ and $k=0,1,\ldots,N$ we compute and store an approximation $F_{i}[k]$ of $\mathcal{F}_{i}(X_{j}[k])$, where node $j$ is the parent of node $i$. Note that this is an approximation of $\mathcal{F}_{i}(X_{j}[k])$ rather than of $\mathcal{F}_{i}(X_{i}[k])$ since $\mathcal{F}_{i}(x)$ is the partial likelihood conditional on the trait value for the parent of node $i$. The value approximation $F_{v}[i]$ is computed by applying the numerical integration method to the appropriate integral (14)—(16), where we replace function evaluations with approximations previously computed. See below for a worked example of this general approach.

The numerical integration methods we use run in time linear in the number of points being evaluated. Hence if $n$ is the number of tips in the tree, the algorithm will run in time $O(nN^{2})$. For the integration techniques described above, the convergence rate (in $N$) for the likelihood on the entire tree had the same order as the convergence rate for the individual one-dimensional integrations (see below for a formal proof of a specific model). We have therefore avoided the computational blow-out typically associated with such high-dimensional integrations, and achieve this without sacrificing accuracy.

The algorithms we have described compute the joint density of the states at the tips, given the tree, the branch lengths, and other parameters. As with discrete traits, the algorithms can be modified to infer ancestral states for internal nodes in the tree. Here we show how to carry out reconstruction of the marginal posterior density of a state at a particular node. The differences between marginal and joint reconstructions are reviewed in (Yang, 2006, pg 121).

First consider marginal reconstruction of ancestral states at the root. Let $u$ and $v$ be the children of the root. The product $\mathcal{F}_{u}(x)\mathcal{F}_{v}(x)$ equals the probability of the observed character conditional on the tree, branch lengths, parameters and a state of $x$ at the root. The marginal probability of $x$, ignoring the data, is given by the root density $f(x|\theta_{r})$. Integrating the product of $\mathcal{F}_{u}(x)\mathcal{F}_{v}(x)$ and $f(x|\theta_{r})$ gives the likelihood $\mathcal{L}(T)$, as in (16). Plugging these into Bayes’ rule, we obtain the posterior density of the state at the root:

$$f(x_{r}|z_{1},\ldots,z_{n})=\frac{\mathcal{F}_{u}(x_{r})\mathcal{F}_{v}(x_{r})f(x_{r}|\theta_{r})}{\mathcal{L}(T)}.$$

(18)

With general time reversible models used in phylogenetics, the posterior distributions at other nodes can be found by changing the root of the tree. Unfortunately the same trick does not work for many qualitative trait models, including the threshold model we study here. Furthermore, recomputing likelihoods for each possible root entails a large amount of unneccessary computation.

Instead we derive a second recursion, this one starting at the root and working towards the tips. A similar trick is used to compute derivatives of the likelihood function in Felsenstein and Churchill (1996). For a node $i$ and state $x$ we let $\mathcal{G}_{i}(x)$ denote the likelihood for the trait values at tips which are not descendants of node $i$, conditional on node $i$ having trait value $x$. If node $i$ is the root $r$, then $\mathcal{G}_{r}(x)$ is $1$ for all $x$.

Let node $i$ be any node apart from the root, let node $j$ be its parent and let node $u$ be the other child of $j$ (that is, the sibling of node $i$). We let $\tilde{x}$ denote the trait value at node $j$. Then $\mathcal{G}_{i}(x)$ can be written

$$\mathcal{G}_{i}(x)=\int f(\tilde{x}|x,\theta_{i})\mathcal{G}_{j}(\tilde{x})\mathcal{F}_{u}(\tilde{x})\,d\tilde{x}.$$

(19)

This integral can be evaluated using the same numerical integrators used when computing likelihoods. Note that $f(\tilde{x}|x,\theta_{i})$ is the conditional density of the parent state given the child state, which is the reverse of the transition densities used to formulate the model. How this is computed will depend on the model and its properties; see below for an implementation of this calculation in the threshold model.

Once $\mathcal{G}_{i}(x)$ has been computed for all nodes, the actual (marginal) posterior densities are computed from Bayes’ rule. Letting $u,v$ be the children of node $i$,

$$f(x_{i}|z_{1},\ldots,z_{n})=\frac{\mathcal{G}_{i}(x_{i})\mathcal{F}_{u}(x_{i})\mathcal{F}_{v}(x_{i})f(x_{i})}{\mathcal{L}(T)}.$$

(20)

In this section we show how the general framework can be applied to the threshold model of Wright (1934) and Felsenstein (2005, 2012). Each trait is modelled by a continuously varying liability which evolves along branches according to a Brownian motion process. While the underlying liability is continuous, the observed data is discrete: at each tip we observe only whether the liability is above or below some threshold.

We will use standard notation for Gaussian densities. Let $\phi(x|\mu,\sigma^{2})$ denote the density of a Gaussian random variable $x$ with mean $\mu$ and variance $\sigma^{2}$; let

$$\Phi(y|\mu,\sigma^{2})=\int_{-\infty}^{y}\phi(x|\mu,\sigma^{2})$$

(21)

denote its cumulative density function, with inverse $\Phi^{-1}(\alpha|\mu,\sigma^{2})$.

Let $X_{1},\ldots,X_{2n-1}$ denote the (unobserved) liability values at the $n$ tips and $n-1$ internal nodes. As above we assume that the $i<j$ whenever node $i$ is a child of node $j$, so that the root has index $2n-1$.

The liability value at the root has a Gaussian density with mean $\mu_{r}$ and variance $\sigma_{r}^{2}$:

$$f(x_{2n-1}|\theta_{r})=\phi(x_{2n-1}|\mu_{r},\sigma_{r}^{2}).$$

(22)

Consider any non-root node $i$ and let $j$ be the index of its parent. Let $t_{i}$ denote the length of the branch connecting nodes $i$ and $j$. Then $X_{i}$ has a Gaussian density with mean $x_{j}$ and variance $\sigma^{2}t_{v}$:

$$f(x_{i}|x_{j},\theta_{i})=\phi(x_{i}|x_{j},\sigma^{2}t_{i}).$$

(23)

Following Felsenstein (2005), we assume thresholds for the tips are all set at zero. We observe $1$ if the liability is positive, $0$ if the liability is negative, and $?$ if data is missing. We can include the threshold step into our earlier framework by defining

$$\pi(z_{i}|x_{i})=\begin{cases}1&\mbox{ if $z_{i}=1$ and $x_{i}>0$, or $z_{i}=0$ and $x_{i}\leq 0$, or $z_{i}=?$}\\ 0&\mbox{ otherwise.}\end{cases}$$

(24)

The likelihood function for observed discrete values $z_{1},\ldots,z_{n}$ is then given by integrating over liability values for all nodes on the tree:

$$\mathcal{L}(T|z_{1},\ldots,z_{n})=\int\displaylimits_{-\infty}^{\infty}\!\cdots\!\int\displaylimits_{-\infty}^{\infty}\phi(x_{2n-1}|\mu_{r},\sigma_{r}^{2})\prod_{(i,j)}\!\!\phi(x_{i}|x_{j},\sigma^{2}t_{i})\prod_{i=1}^{n}\pi(z_{i}|x_{i})\,\,\mathrm{d}x_{1}\ldots\,\mathrm{d}x_{2n-1}.$$

(25)

The first step towards computing $\mathcal{L}(T|z_{1},\ldots,z_{n})$ is to bound the domain of integration so that we can apply Simpson’s method. Ideally, we would like these bounds to be as tight as possible, for improved efficiency. For the moment we will just outline a general procedure which can be adapted to a wide range of evolutionary models.

The marginal (prior) density of a single liability or trait value at a single node is the density for that liability value marginalizing over all other values and data. With the threshold model, the marginal density for the liability at node $i$ is Gaussian with mean $\mu_{r}$ (like the root) and variance $v_{i}$ equal to the sum of the variance at the root and the transition variances on the path from the root to node $i$. If $P_{i}$ is the set of nodes from the root to node $i$, then

$$v_{i}=\sigma_{r}^{2}+\sigma^{2}\sum_{j\in P_{i}}t_{j}.$$

(26)

The goal is to constrain the error introduced by truncating the integrals with infinite domain. Let $\epsilon$ be the desired bound on this truncation error. Recall that the number of internal nodes in the tree is $n-1$. Define

$$L_{i}=\Phi^{-1}\left(\frac{\epsilon}{2(n-1)}\Big{|}\mu_{r},v_{i}\right)$$

(27)

and

$$U_{i}=\Phi^{-1}\left(\frac{\epsilon}{2(n-1)}\Big{|}\mu_{r},v_{i}\right)$$

(28)

so that the probability $X_{i}$ lies outside the interval $[L_{i},U_{i}]$ is at most $\epsilon/(n-1)$. By the inclusion-exclusion principle, the joint probability $X_{i}\not\in[L_{i},U_{i}]$ for any internal node $i$ is at most $\epsilon$. We use this fact to bound the contribution of the regions outside these bounds.

	$\displaystyle\int\displaylimits_{-\infty}^{\infty}$	$\displaystyle\!\cdots\!\int\displaylimits_{-\infty}^{\infty}f(x_{2n-1}|\mu_{r},\sigma_{r}^{2})\prod_{(u,v)}\!\!f(x_{v}|x_{u},\theta_{v})\prod_{i=1}^{n}\pi(z_{i}|x_{i})\,\,\mathrm{d}x_{1}\ldots\,\mathrm{d}x_{2n-1}$		(29)
		$\displaystyle\qquad-\int\displaylimits_{a_{2n-1}}^{b_{2n-1}}\!\!\cdots\!\int\displaylimits_{a_{n+1}}^{b_{n+1}}\int\displaylimits_{-\infty}^{\infty}\!\cdots\!\int\displaylimits_{-\infty}^{\infty}f(x_{2n-1}|\mu_{r},\sigma_{r}^{2})\prod_{(u,v)}\!\!f(x_{v}|x_{u},\theta_{v})\prod_{i=1}^{n}\pi(z_{i}|x_{i})\,\,\mathrm{d}x_{1}\ldots\,\mathrm{d}x_{2n-1}$		(30)
		$\displaystyle\leq\int\displaylimits_{-\infty}^{\infty}\!\cdots\!\int\displaylimits_{-\infty}^{\infty}f(x_{2n-1}|\mu_{r},\sigma_{r}^{2})\prod_{(u,v)}\!\!f(x_{v}|x_{u},\theta_{v})\,\,\mathrm{d}x_{1}\ldots\,\mathrm{d}x_{2n-1}$		(31)
		$\displaystyle\qquad-\int\displaylimits_{a_{2n-1}}^{b_{2n-1}}\!\!\cdots\!\int\displaylimits_{a_{n+1}}^{b_{n+1}}\int\displaylimits_{-\infty}^{\infty}\!\cdots\!\int\displaylimits_{-\infty}^{\infty}f(x_{2n-1}|\mu_{r},\sigma_{r}^{2})\prod_{(u,v)}\!\!f(x_{v}|x_{u},\theta_{v})\,\,\mathrm{d}x_{1}\ldots\,\mathrm{d}x_{2n-1}$		(32)
		$\displaystyle\leq P\Big{(}X_{n+1}\not\in[L_{n+1},U_{n+1}]\mbox{ or }X_{n+2}\not\in[L_{n+2},U_{n+2}]\mbox{ or }\cdots\mbox{ or }X_{2n-1}\not\in[L_{2n-1},U_{2n-1}]\Big{)}$		(33)
		$\displaystyle<\epsilon.$		(34)

We therefore compute values $L_{i},U_{i}$ for $n+1\leq i\leq 2n-1$ and use these bounds when carrying out integration at the internal nodes. We define

$$X_{i}[k]=L_{i}+\frac{U_{i}-L_{i}}{N}k$$

(35)

for $k=0,1,\ldots,N$ for each internal node $i$.

The next step is to use dynamic programming and numerical integration to compute the approximate likelihood. Let node $i$ be a tip of the tree, let node $j$ be its parent and let $z_{i}$ be the binary trait value at this tip. For each $k=0,1,\ldots,N$ we use standard error functions to compute

	$\displaystyle F_{i}[k]$	$\displaystyle=$	$\displaystyle\mathcal{F}_{i}(X_{j}[k])$		(36)
		$\displaystyle=$	$\displaystyle\begin{cases}\int\displaylimits_{0}^{\infty}\phi(\tilde{x}|X_{j}[k],\sigma^{2}t_{i})\,\,\mathrm{d}\tilde{x}&\mbox{ if $z_{i}=1$}\\ \int\displaylimits_{-\infty}^{0}\phi(\tilde{x}|X_{j}[k],\sigma^{2}t_{i})\,\,\mathrm{d}\tilde{x}&\mbox{ if $z_{i}=0$}\\ 1&\mbox{ if $z_{i}=?$.}\end{cases}$		(37)

Here $\phi(x|\mu,\sigma^{2})$ is the density of a Gaussian with mean $\mu$ and variance $\sigma^{2}$.

Now suppose that node $i$ is an internal node with parent node $j$ and children $u$ and $v$. Applying Simpson’s rule to the bounds $L_{i},U_{i}$ to (15) we have for each $k=0,1,\ldots,N$:

	$\displaystyle F_{i}[k]$	$\displaystyle=$	$\displaystyle\frac{U_{i}-L_{i}}{N}\sum\displaylimits_{\ell=0}^{N}w_{\ell}\phi(X_{i}[\ell]|X_{j}[k],\sigma^{2}t_{i})F_{u}[\ell]F_{v}[\ell]$		(38)
		$\displaystyle\approx$	$\displaystyle\mathcal{F}_{i}(X_{j}[k]).$		(39)

Suppose node $i$ is the root, and $u,v$ are its children. Applying Simpson’s rule to (16) gives

	$\displaystyle L$	$\displaystyle\leftarrow$	$\displaystyle\frac{U_{2n-1}-L_{n-1}}{N}\sum_{\ell=0}^{N}w_{\ell}\phi(X_{i}[\ell]|\mu_{r},\sigma_{r}^{2})F_{u}[\ell]F_{v}[\ell]$		(40)
		$\displaystyle\approx$	$\displaystyle\mathcal{L}(T|z_{1},\ldots,z_{n},\theta).$		(41)

Pseudo-code for the algorithm appears in Algorithm 1. Regarding efficiency and convergence we have:

Simpson’s rule has $O(N^{-4})$ convergence on functions with bounded fourth derivatives (Dahlquist and Björck, 2008). The root density and each of the transition densities are Gaussians, so have individually have bounded fourth derivatives. For each node $i$, let $n_{i}$ denote the number of tips which are descendents of the node. Using induction on (15), we see that for all nodes $i$, the fourth derivate of $\mathcal{F}_{i}(x)$ is $O(n_{i})$.

Letting $\epsilon=nN^{-4}$ we have from above that the error introduced by replacing the infinite domain integrals with integrals on $[L_{i},U_{i}]$ introduces at most $nN^{-4}$ error. Using a second induction proof on (15) and (38) together with the bound on fourth derivatives we have that $|\mathcal{F}_{i}(X_{j}[k])-F_{i}[k]|$ is at most $O(n_{i}N^{-4})$ for all nodes $i$, where node $j$ is the parent of node $i$. In this way we obtain error bound of $O(n_{2n-1}N^{-4})=O(nN^{-4})$ on the approximation of $\mathcal{L}(T|z_{1},\ldots,z_{n},\theta)$. $\Box$

We can estimate posterior densities using the recursion (19) followed by equation (20). The conditional density

$$f(\tilde{x}|x,\theta_{i})=\phi\left(\tilde{x}\Big{|}\mu_{r}+\frac{\sigma_{r}^{2}+\sigma^{2}P_{j}}{\sigma_{r}^{2}+\sigma^{2}P_{i}}\left(x-\mu_{r}\right),\frac{\sigma^{2}t_{i}\left(\sigma_{r}^{2}+\sigma^{2}P_{j}\right)}{\sigma_{r}^{2}+\sigma^{2}P_{i}}\right)$$

(42)

can be obtained by plugging the transitional density

$$f(x|\tilde{x},\theta_{i})=\phi(x|\tilde{x},\sigma^{2}t_{i})$$

(43)

and the two marginal densities (26)

$$f(\tilde{x})=\phi(\tilde{x},\sigma_{r}^{2}+\sigma^{2}P_{j}),\quad f(x)=\phi(x,\sigma_{r}^{2}+\sigma^{2}P_{i})$$

(44)

into the identity $f(\tilde{x}|x,\theta_{i})=f(x|\tilde{x},\theta_{i})\frac{f(\tilde{x})}{f(x)}$. We thereby obtain the recursion

$$\mathcal{G}_{i}(x)=\int\phi\left(\tilde{x}\Big{|}\mu_{r}+\frac{\sigma_{r}^{2}+\sigma^{2}P_{j}}{\sigma_{r}^{2}+\sigma^{2}P_{i}}\left(x-\mu_{r}\right),\frac{\sigma^{2}t_{i}\left(\sigma_{r}^{2}+\sigma^{2}P_{j}\right)}{\sigma_{r}^{2}+\sigma^{2}P_{i}}\right)\mathcal{G}_{j}(\tilde{x})\mathcal{F}_{u}(\tilde{x})\,d\tilde{x}$$

(45)

which we estimate using Simpson’s method. Algorithm estimates values of the posterior densities at each node, evaluated using the same set of grid points as used in Algorithm 1. An additional round of numerical integration can be used to obtain posterior means and variances.

To study the methods in practice, we reanalyse trait data published by Marazzi et al. (2012), using a fixed phylogeny. Marazzi et al. (2012) introduce and apply a new discrete state model for morphological traits which, in addition to states for presence and absence, incorporates an intermediate ‘pre-cursor’ state. Whenever the intermediate state is observed at the tips it is coded as ‘absent’. The motivation behind the model is that the intermediate state represents evoutionary pre-cursors, changes which are necessary for the evolution of a new state but which may not be directly observed. These pre-cursors could explain repeated parallel evolution of a trait in closely related traits (Marazzi et al., 2012). They compiled a data set recording presence or absence of plant extrafloral nectaries (EFNs) across a phylogeny of 839 species of Fabales, fitting their models to these data.

The threshold model also involves evolutionary pre-cursors in terms of changes in ancestral liabilities. We use these models, and our new algorithms to analyse the EFN dataset. Our analysis also makes use of the time-calibrated phylogeny inferred by Simon et al. (2009), although unlike Marazzi et al. (2012) we ignore phylogenetic uncertainty.

We conduct three separate experiments. For the first experiment, we examine the rate of convergence of the likelihood algorithm as we increase $N$. This is done for the ‘All’ EFN character (Character 1 in Marazzi et al. (2012)) for a range of estimates for the liability variance at the root, $\sigma_{r}^{2}$. The interest in $\sigma_{r}^{2}$ stems from its use in determining bounds $L_{i},U_{i}$ for each node, with the expectation that as $\sigma_{r}^{2}$ increases, the convergence of the integration algorithm will slow. The mean liability at the root, $\mu_{r}$ was determined from the data using Maximum Likelihood estimation.

We also examined convergence of the algorithm on randomly generated characters. We first evolved liabilities according to the threshold model, using the parameter settings obtained above. To examine the difference in performance for non-phylogenetic characters we also simulated binary characters by simulated coin flipping. Twenty replicates were carried out for each case.

The second experiment extends the model comparisons carried out in Marazzi et al. (2012) to include the threshold models. For this comparison we fix the transitional variance $\sigma^{2}$ at one, since changing this values corresponds to a rescaling of the Brownian process, with no change in likelihood. With only one character, the maximum likelihood estimate of the root variance $\sigma_{r}^{2}$ is zero, irrespective of the data. This leaves a single parameter to infer: the value of the liability at the root state. We computed a maximum likelihood estimate for the state at the root, then applied our algorithm with a sufficiently large value of $N$ to be sure of convergence. The Akaike Information Criterion (AIC) was determined and compared with those obtained for the model of Marazzi et al. (2012).

For the third experiment, we determine the marginal posterior densities for the liabilities at internal nodes, using Algorithm 2. These posterior probabilities are then mapped onto the phylogeny, using shading to denote the (marginal) posterior probability that a liability is larger than zero. We therefore obtain a figure analogous to Supplementary Figure 7 of Marazzi et al. (2012).