Thursday, January 30, 2014

Function for midpoint rooting

Today Todd Oakley asked:

"Does anyone know of an existing midpoint rooting routine? I am displaying trees and would like to show them as midpoint rooted. I've been using the phangorn package, which does the midpoint rooting perfectly, but it has some dependencies that make it unstable on the linux machines I've been using. When I looked last, phangorn was the only package I could find with midpoint rooting."

Well, midpoint rooting is not theoretically very difficult. All one needs to do is find the longest path between any pair of tips & then locate the root midway along that path. I just posted code for this here (as well as a new phytools build, which can be downloaded & installed from source).

Here is a brief description of some of the tricks that I used:

(1) I used cophenetic.phylo from ape to get all the distances between tips, choose the longest one, and then identify the two species on either end of that path.

(2) I used reroot in phytools to re-root the tree immediately below one of these two tips. I did this so that I could use a new custom function getAncestors (which does the same thing as Ancestors in phangorn, but using no phangorn code) to find the set of all internal nodes ancestral to the other tip subtending the longest path. The new position of the root will be between two of these nodes.

(3) I used dist.nodes in ape to compute the distances between the tip of interest and all the internal nodes I had found in (2).

(4) Finally, I found the two nodes subtending the new root position & re-rooted the tree (using reroot from phytools) in the correct position between those nodes.

That's it. The function has not yet been thoroughly tested, so please give me feedback if it doesn't work as intended. Here's a quick demo:

> library(phytools)
Loading required package: ape
Loading required package: maps
Loading required package: rgl
> packageVersion("phytools")
[1] ‘0.3.89’
> tree<-rtree(n=12)
> mpt1<-midpoint.root(tree)
> plotTree(mpt1)
> require(phangorn)
Loading required package: phangorn
> mpt2<-midpoint(tree)
> plotTree(mpt2)

These trees may not look exactly the same, but they are (just with different rotations of internal nodes):

> all.equal.phylo(mpt1,mpt2)
[1] TRUE

I have not tested this function against phangorn's midpoint, but if history is any indication, midpoint is probably faster & more elegantly programmed. Hopefully for Todd's purposes, this will work.

New version of ancThresh for λ model

The current version of ancThresh (for ancestral character estimation under the threshold model; see Revell In press) permits a Brownian or Ornstein-Uhlenbeck model for the evolution of the liabilities. I just added a 3rd model, the λ model of Pagel (1999). The code for this version is here; but since it uses some functions internally that have also been updated, the best thing to do is to update phytools to the latest version.

I'm not a huge fan of the λ model in general, since it is not clear what biological process it is meant to approximate; however under some circumstances it could be a useful model for traits that evolve on the tree (and thus have phylogenetic covariance), but are also affected by contemporary factors that are not necessarily phylogenetically correlated. This is how I am using this model in the empirical study for which I have added this feature to ancThresh.

At the moment it is largely untested - so if you run into any problems, please let me know.

Saturday, January 25, 2014

Bug fix in phylANOVA

A phytools user recently reported discrepant results between phylANOVA in phytools and aov.phylo in geiger. Both functions conduct the simulation-based method of Garland et al. (2013). In theory, the only difference is that phylANOVA performs post-hoc comparison of means. It turns out, however, that phylANOVA contains the implicit assumption that y is in the order of tree$tip.label. This assumption is now only true if names(y) is NULL, in which case a warning is also issued. Updated code is here and in the latest version of phytools.