So, in addition to tracking new development activity with this page, I thought I would also post a few words about some of the existing functions that I am distributing from my R phylogenetics page.Today, I will discuss the function evol.rate.mcmc(). This function takes a continuous character vector and a phylogeny with branch lengths, and then uses Bayesian MCMC to identify the phylogenetic position of a shift in the evolutionary rate of phenotypic evolution through time. Evolutionary rates and shift-points between rates are sampled from their joint posterior distribution. The underlying method for this function is described in a manuscript presently in review:
Revell, L. J., D. L. Mahler, P. R. Peres-Neto, and B. D. Redelings. Submitted. A new method for identifying exceptional evolutionary diversification.
To give a quick example of how this works, take the data and tree presented at right. These data were generated with a high rate (σ2=20) on the red branches of the tree; and a low rate (σ2=1.0) on the blue branches. The data are in a labeled vector, x, and the phylogenetic tree is given in an {ape} "phylo" object tree.
To load the source and required dependencies, we will first enter:
> source("evol.rate.mcmc.R")
> require(ape); require(coda);
The package {coda} will be used to analyze and compute various diagnostics for our MCMC chain (next post) - in this case our posterior sample - but we may as well load it now.
To run the function, we simply give R (some variant of) the following command:
> res<-evol.rate.mcmc(tree,x,ngen=20000,control=list(sd1=1.0, sd2=1.0,sda=0.5,sdlnr=2.0,rand.shift=0.05,print=1000))
This tells our MCMC to run under the following conditions:
ngen: the number of generations for our MCMC run is 20,000;
sd1: standard deviation of the proposal distribution for rate 1 is 1.0;
sd2: standard deviation of the proposal distribution for rate 2 is 1.0;
sda: standard deviation on the proposal distribution for the ancestral state at the root node of the tree is 0.5.
rand.shift: probability of a random shift to a different position in the tree is 0.05 [otherwise the MCMC will proceed as a random walk with a default variance];
print: print frequency for the chain (to screen - all generations are stored in res$mcmc).
We also need to assign a prior probability distribution to our parameters. I chose to specify a uniform prior on the absolute values of rates 1 and 2, the root state, and the position of the rate shift; however, I have used a log-normal prior distribution for the ratio of rate 1 / rate 2, centered on 0.0, and with a standard deviation given by sdlnr - in this case we set this to 2.0.
The run returns a list with two components:
$mcmc: is our matrix of parameters from our MCMC run, including the sample of shift points between rates on the tree; and
$tips: is a list of vectors containing all the tips in the clade defined as rate 1.
In my next post, I will discuss analysis of the posterior sample from this run.
Liam, very nice post example. I'm unclear how the shift is being modeled -- there is constant probability rand.shift that at any moment in time a shift in rate occurs, with transition from sd1->sd2 occurring at the same rate as sd2->sd1?
ReplyDeleteThis assumes the trait diversification rate (sd1 or 2) is determined by a discrete character with constant known transition rate rand.shift? Would it be possible to MCMC over possible values of rand.shift as well, given some prior for it?
Exciting developments. I was working over the difference between this and Harmon/Eastman's RJMCMC of diversification rates -- I guess they can relax the assumption of two rates and don't infer a transition rate between the faster and slower diversification, but rather just infer transition points?
Hi Carl. Sorry this was unclear. "rand.shift" is just for the MCMC chain. [I.e., the chain samples shift points between rates as a random walk, but then jumps to a random position in the tree with probability "rand.shift." We did this mainly to help improve mixing.]
ReplyDeleteOur model is one in which the Brownian rate changes once at some point somewhere in the tree. It is meant as a starting point to model shifts in the OU process and multiple shifts in a similar way. It differs from related methods, such as O'Meara et al. (2006), in that we bring in no extrinsic information about where the rate changes in the tree. This is inferred only from our data for the continuous character.
I saw an early version of the Eastman et al. manuscript (which, I will note, seemed very impressive) shortly after submitting our paper. Since their article is unpublished, I won't discuss its details except to comment that indeed the two methods are related.
Is this method valid to apply to datasets containing non-coeval tips? (A stereotypical question for a paleontologist, I know.)
ReplyDeleteChris Venditti presented an abstract on a MCMC method for identifying shifts in the rate of BM evolution on a tree, at the International Paleontological Congress in London last summer. Meade and Pagel were his co-authors. He applied the method to body-size evolution on the dinosaur supertree. I guess it's just an idea whose time had come?
-Dave
Hi Dave.
ReplyDeleteYes, the method is completely appropriate to a phylogeny with non-contemporaneous tips.
I don't know of the Venditti method and can't find it on his website. I would be happy to learn more.
You're probably right that the time has come for this kind of approach. I starting working on this just as soon as I knew how. . . .
By the way, if you try to run this with your data, please let me know if it works and if I can be of any help in tuning the control parameters of the function.
ReplyDeleteThe result matrix, $mcmc, can be printed to a tab delimited text file and then read into Tracer for typical MCMC diagnostics. [These will be more informative for the likelihoods and quantitative parameters than for the phylogenetic shift-position.] If you want to stick with R, most of the same diagnostics can be computed using functions from the package "coda." I will cover some of this in a subsequent post.
Hey Liam,
ReplyDeleteCool stuff. One possible extension would be to have, say, two continuous traits and to look for shifts in trait correlation (i.e. to investigate the degree of correlation itself as the trait of interest). That could be a great tool to study phenotypic integration/dissociation.
Does that make sense to you?
Hi Marcio.
ReplyDeleteI think this is a great idea. I already have a likelihood method (described in Revell & Collar 2008; Evol.) and R function [evol.vcv()] for fitting multiple evolutionary correlations, conditioned on a discrete character map. A natural next step using MCMC would be to find the shift-points based on only the data for our continuously valued traits.
Liam,
ReplyDeleteHow would you apply this method to datasets with multiple valid phylogenetic hypotheses? Average the parameter estimates across MCMC analyses on each tree?
-Dave
Dave,
ReplyDeleteThis is a really good question - and not so easy to answer. Part of the posterior sample in this case is a position on an edge of the tree - the location of our Brownian rate shift.
If you have multiple trees (say, from a Bayesian posterior sample), then you can imagine having posterior density for the rate shift on some edges of one tree that are absent from a second. I think in this case you might want to compute the posterior probability of the change on each edge as the average probability across trees - setting this probability to zero if the edge is absent from the tree.
This is just an idea, though.