Sunday, May 22, 2011

Fitting two or more evolutionary correlations to the tree and data

As of late, I am in the field and haven't had time to do much coding (although I am slowly working on implemented the phylogenetic paired t-test of my article with Patrik Lindenfors and Charlie Nunn); however, I thought I would create a quick note about a function that I have available on my R-phylogenetics page, but that I have not described in any earlier blog posts.

The function evol.vcv() fits the model of Revell & Collar (2009). In this article we fit two different BM variance-covariance matrices (we call these "evolutionary rate matrices" - because the diagonal elements are the evolutionary rates, by some account), to different parts of the tree that have been specified a priori by the user. In my implementation (direct link to code here), the user can fit an arbitrary number of traits as well as an arbitrary number of rate matrix categories. These are "painted" on the tree using SIMMAP v1.0 format Newick trees that are read into memory using the function read.simmap() (code here; prior descriptions of the development of this function in several posts).

One of the neat things about this model is that it allows the user to fit different evolutionary correlations between characters to different branches and parts of branches in the tree. For instance, in our study we fit this model for piscivorous and non-piscivorous fish species and found that the evolutionary correlation between two different aspects of buccal morphology was much higher in piscivorous fish species.

To try the function, first load the source:

> source("http://faculty.umb.edu/liam.revell/phytools/evol.vcv/v0.1/evol.vcv.R") # or download

Then, simply read in a data matrix, X, which can contain an arbitrary number of columns, e.g.:

> X<-read.csv(file="data.csv",row.names=1)

Then read the tree in simmap v1.0 format:

> tree<-read.simmap(file="treefile.tre",format="nexus") # for instance

Finally, execute the analysis:

> result<-evol.vcv(tree,X,maxit=5000)

The output is the fit and fitted parameters of both a one matrix and two evolutionary rate matrix model, as well as a comparison of the models.

No comments:

Post a Comment