Wednesday, November 8, 2017

(New) print method for (old) Bayesian MCMC function to identify the phylogenetic position of a rate-shift

A number of years ago some colleagues & I described a new Bayesian method (Revell et al. 2012) for identifying shifts in the rate of phenotypic evolution on a phylogeny.

Believe it or not, at the time we thought that this was cutting-edge - and in a sense it was. To that point, there had been very limited use of Bayesian MCMC in phylogenetic comparative methods, particularly for continuous traits. Furthermore, though techniques to test a priori hypotheses about changes in the evolutionary rate or process for continuous traits were already well-developed, no one (at that point) had published a technique for identifying shifts absent a prior hypothesis. This, of course, changed rapidly as members of the Harmon Lab at the University of Idaho published various rjMCMC (reversible-jump) Bayesian methods designed to sample rates and the location of rate-shifts from their joint posterior distribution.

Nonetheless, our method is still valid. The function implementing the method, evol.rate.mcmc, has long been part of the phytools package. Today, I decided to create a custom object class & add a print method to the function - just to make it a little easier to use.

Here's a quick demo:

## Loading required package: ape
## Loading required package: maps

plot of chunk unnamed-chunk-1

## Starting MCMC run....
## Done MCMC run.
## Object of class "evol.rate.mcmc" containing the results from a
## the Bayesian MCMC analysis of a Brownian-motion rate-shift model.
## MCMC was conducted for 1e+05 generations sampling every 100 
## generations.
## The most commonly sampled rate shift(s) occurred on the edge(s)
## to node(s) 121.
## Use the functions posterior.evolrate and minSplit for more detailed
## analysis of the posterior sample from this analysis.

In fact, the method has easily identified the precise edge on which I simulated a rate-shift. Here is a plot showing the simulated shift & then I will calculate the rate-shift location with the minimum distance to the rest in the posterior sample & overlay it.

text(xx[2],yy[2],"location of the minimum\ndistance split",pos=2,

plot of chunk unnamed-chunk-2

Neat - we're almost right on.

I simulated the tree & data as follows (although I ran it a few times to ensure that the shift was not on a terminal edge or for a very small clade):

nchanges<-function(tree) sum(sapply(tree$maps,length))-nrow(tree$edge)
while(nchanges(mtree)!=1) nchanges(mtree<-sim.history(tree<-foo(),Q))

No comments:

Post a Comment