Bayesian inference in phylogeny

One of the most common MCMC methods used is the Metropolis–Hastings algorithm,[8] a modified version of the original Metropolis algorithm.[9] It is a widely used method to sample randomly from complicated and multi-dimensional distribution probabilities. The Metropolis algorithm is described in the following steps:[10] [11]

1. An initial tree, T_i, is randomly selected.
2. A neighbour tree, T_j, is selected from the collection of trees.
3. The ratio, R, of the probabilities (or probability density functions) of T_j and T_i is computed as follows: R = f(T_j)/f(T_i)
4. If R ≥ 1, T_j is accepted as the current tree.
5. If R < 1, T_j is accepted as the current tree with probability R, otherwise T_i is kept.
6. At this point the process is repeated from Step 2 N times.

The algorithm keeps running until it reaches an equilibrium distribution. It also assumes that the probability of proposing a new tree T_j when we are at the old tree state T_i, is the same probability of proposing T_i when we are at T_j. When this is not the case Hastings corrections are applied. The aim of Metropolis-Hastings algorithm is to produce a collection of states with a determined distribution until the Markov process reaches a stationary distribution. The algorithm has two components:

1. A potential transition from one state to another (i → j) using a transition probability function q_i,j
2. Movement of the chain to state j with probability α_i,j and remains in i with probability 1 – α_i,j.[2]

Metropolis-coupled MCMC algorithm (MC³) [12] has been proposed to solve a practical concern of the Markov chain moving across peaks when the target distribution has multiple local peaks, separated by low valleys, are known to exist in the tree space. This is the case during heuristic tree search under maximum parsimony (MP), maximum likelihood (ML), and minimum evolution (ME) criteria, and the same can be expected for stochastic tree search using MCMC. This problem will result in samples not approximating correctly to the posterior density. The (MC³) improves the mixing of Markov chains in presence of multiple local peaks in the posterior density. It runs multiple (m) chains in parallel, each for n iterations and with different stationary distributions ${\displaystyle \pi _{j}(.)\ }$, ${\displaystyle j=1,2,\ldots ,m\ }$, where the first one, ${\displaystyle \pi _{1}=\pi \ }$ is the target density, while ${\displaystyle \pi _{j}\ }$, ${\displaystyle j=2,3,\ldots ,m\ }$ are chosen to improve mixing. For example, one can choose incremental heating of the form:

${\displaystyle \pi _{j}(\theta )=\pi (\theta )^{1/[1+\lambda (j-1)]},\ \ \lambda >0,}$

so that the first chain is the cold chain with the correct target density, while chains ${\displaystyle 2,3,\ldots ,m}$ are heated chains. Note that raising the density ${\displaystyle \pi (.)}$ to the power ${\displaystyle 1/T\ }$ with ${\displaystyle T>1\ }$ has the effect of flattening out the distribution, similar to heating a metal. In such a distribution, it is easier to traverse between peaks (separated by valleys) than in the original distribution. After each iteration, a swap of states between two randomly chosen chains is proposed through a Metropolis-type step. Let ${\displaystyle \theta ^{(j)}\ }$ be the current state in chain ${\displaystyle j\ }$, ${\displaystyle j=1,2,\ldots ,m\ }$. A swap between the states of chains ${\displaystyle i\ }$ and ${\displaystyle j\ }$ is accepted with probability:

${\displaystyle \alpha ={\frac {\pi _{i}(\theta ^{(j)})\pi _{j}(\theta ^{(i)})}{\pi _{i}(\theta ^{(i)})\pi _{j}(\theta ^{(j)})}}\ }$

At the end of the run, output from only the cold chain is used, while those from the hot chains are discarded. Heuristically, the hot chains will visit the local peaks rather easily, and swapping states between chains will let the cold chain occasionally jump valleys, leading to better mixing. However, if ${\displaystyle \pi _{i}(\theta )/\pi _{j}(\theta )\ }$ is unstable, proposed swaps will seldom be accepted. This is the reason for using several chains which differ only incrementally.

An obvious disadvantage of the algorithm is that ${\displaystyle m\ }$ chains are run and only one chain is used for inference. For this reason, ${\displaystyle \mathrm {MC} ^{3}\ }$ is ideally suited for implementation on parallel machines, since each chain will in general require the same amount of computation per iteration.

The LOCAL algorithms[13] offers a computational advantage over previous methods and demonstrates that a Bayesian approach is able to assess uncertainty computationally practical in larger trees. The LOCAL algorithm is an improvement of the GLOBAL algorithm presented in Mau, Newton and Larget (1999)[14] in which all branch lengths are changed in every cycle. The LOCAL algorithms modifies the tree by selecting an internal branch of the tree at random. The nodes at the ends of this branch are each connected to two other branches. One of each pair is chosen at random. Imagine taking these three selected edges and stringing them like a clothesline from left to right, where the direction (left/right) is also selected at random. The two endpoints of the first branch selected will have a sub-tree hanging like a piece of clothing strung to the line. The algorithm proceeds by multiplying the three selected branches by a common random amount, akin to stretching or shrinking the clothesline. Finally the leftmost of the two hanging sub-trees is disconnected and reattached to the clothesline at a location selected uniformly at random. This would be the candidate tree.

Suppose we began by selecting the internal branch with length ${\displaystyle t_{8}\ }$ that separates taxa ${\displaystyle A\ }$ and ${\displaystyle B\ }$ from the rest. Suppose also that we have (randomly) selected branches with lengths ${\displaystyle t_{1}\ }$ and ${\displaystyle t_{9}\ }$ from each side, and that we oriented these branches. Let ${\displaystyle m=t_{1}+t_{8}+t_{9}\ }$, be the current length of the clothesline. We select the new length to be ${\displaystyle m^{\star }=m\exp(\lambda (U_{1}-0.5))\ }$, where ${\displaystyle U_{1}\ }$ is a uniform random variable on ${\displaystyle (0,1)\ }$. Then for the LOCAL algorithm, the acceptance probability can be computed to be:

${\displaystyle {\frac {h(y)}{h(x)}}\times {\frac {{m^{\star }}^{3}}{m^{3}}}\ }$

To estimate a branch length ${\displaystyle t}$ of a 2-taxon tree under JC, in which ${\displaystyle n_{1}}$ sites are unvaried and ${\displaystyle n_{2}}$ are variable, assume exponential prior distribution with rate ${\displaystyle \lambda \ }$. The density is ${\displaystyle p(t)=\lambda e^{-\lambda t}\ }$. The probabilities of the possible site patterns are:

${\displaystyle 1/4\left(1/4+3/4e^{-4/3t}\right)\ }$

for unvaried sites, and

${\displaystyle 1/4\left(1/4-1/4e^{-4/3t}\right)\ }$

Thus the unnormalized posterior distribution is:

${\displaystyle h(t)=\left(1/4\right)^{n_{1}+n_{2}}\left(1/4+3/4{e^{-4/3t}}^{n_{1}}\right)\ }$

or, alternately,

${\displaystyle h(t)=\left(1/4-1/4{e^{-4/3t}}^{n_{2}}\right)(\lambda e^{-\lambda t})\ }$

Update branch length by choosing new value uniformly at random from a window of half-width ${\displaystyle w\ }$ centered at the current value:

${\displaystyle t^{\star }=|t+U|\ }$

where ${\displaystyle U\ }$is uniformly distributed between ${\displaystyle -w\ }$ and ${\displaystyle w\ }$. The acceptance probability is:

${\displaystyle h(t^{\star })/h(t)\ }$

Example: ${\displaystyle n_{1}=70\ }$, ${\displaystyle n_{2}=30\ }$. We will compare results for two values of ${\displaystyle w\ }$, ${\displaystyle w=0.1\ }$ and ${\displaystyle w=0.5\ }$. In each case, we will begin with an initial length of ${\displaystyle 5\ }$ and update the length ${\displaystyle 2000\ }$ times.