Fast ML estimation of ancestral states for a continuously valued trait

I just wrote a new function for ancestral character estimation that takes advantage of the fact that the ancestral value for the root node computed during the calculation of independent contrasts is also the MLE of the root. Re-rooting at all internal nodes and the recomputing the PIC root state, cumbersome as it sounds, is actually much faster than finding ancestral states via numerical optimization of the likelihood. The main part of the code for this function is therefore as follows:

M<-btree$Nnode
N<-length(btree$tip)
anc<-vector()
for(i in 1:M+N){
  anc[i-N]<-ace(x,multi2di(root(btree,node=i)), method="pic")$ace[1]
  names(anc)[i-N]<-i
}

Obviously, although this function is designed as a faster version of ace(...,type="continuous",method="ML"), the function actually runs by using many calls to ace(...,method="pic")!!

The object name btree in this code is used to denote binary tree. This highlights the fact that, obviously, contrasts can only be computed for bifurcating trees. We would still like our function to run, though, if the tree contains multifurcating. The remainder of the function code is dedicated to, first, changing a multifurcating tree to a bifurcating one; and then, after computing ancestral states using the code given above, lining up the nodes on the binary tree to the nodes of the original tree containing multifurcations.

The way a did this is a little ad hoc, but it seems to work. I basically went through all the nodes on both the binary and multifurcating tree and, for each node, pulled out a list of all the descendant tips. Then, I used these lists of descendant tips to match nodes between trees.

The code for this function, fastAnc, is here. Let's try it and compare to both ape::ace and phytools::anc.ML.

> library(phytools)
> # first, simulate bifurcating tree & data
> tree<-rtree(200)
> x<-fastBM(tree)
> # now load the source & estimate states with each function
> source("fastAnc.R")
> system.time(res1<-ace(x,tree,method="ML")$ace)
  user  system elapsed
 19.69    0.05   20.38
> system.time(res2<-anc.ML(tree,x,maxit=10000)$ace)
  user  system elapsed
236.67    0.34  255.00
> system.time(res3<-fastAnc(tree,x))
  user  system elapsed
  2.60    0.00    2.98
> # plot to compare
> par(mfrow=c(2,1),mai=c(1.02,0.82,0.1,0.1))
> plot(res1,res3,xlab="ace",ylab="fastAnc")
> plot(res2,res3,xlab="anc.ML",ylab="fastAnc")

Obviously, we get the same estimates from each function - but at much greater computational cost (particularly for anc.ML, which is pretty terrible).

That was for a fully bifurcating tree. We can also do a smaller example to make sure we are back-translating our node IDs correctly. Let's try it:

> tree<-rtree(12)
> # collapse two shortest branches into multifurcations
> tree<-di2multi(tree,tol= sort(tree$edge.length[tree$edge[,2]>length(tree$tip)])[3])
> plot(tree,no.margin=TRUE); nodelabels()

> # simulate data
> x<-fastBM(tree)
> # estimate ancestral states using the three methods
> res1<-ace(x,tree,method="ML")$ace
Error in ace(x, tree, method = "ML") :
 "phy" is not rooted AND fully dichotomous.
> res2<-anc.ML(tree,x)$ace
> res3<-fastAnc(tree,x)
> par(mai=c(1.02,0.82,0.2,0.2))
> plot(res2,res3,xlab="anc.ML",ylab="fastAnc")
> res2
        13          14          15          16          17
0.28663746  0.05138613 -0.25089177  0.50494223  0.13529718
        18          19          20          21
0.03920500  0.49088849  1.24192050  1.63727735
> res3
        13          14          15          16          17
0.28665585  0.05139666 -0.25089395  0.50494620  0.13528936
        18          19          20          21
0.03920408  0.49088277  1.24191236  1.63727663

First off - ace doesn't work at all if the tree is not bifurcating (not sure why this is). anc.ML works fine, but will be very slow for large trees (as we've discovered). fastAnc gives the same estimates and node names - even though it had to first convert to a bifurcating tree, and then back-translate the node numbers. Great!

New version of phytools on CRAN

I just submitted a new version of phytools (version 0.2-0) to CRAN. This is basically the same as the latest "non-CRAN" release (0.1-98) with a couple of minor changes and corrections in the manual pages and to the function fancyTree. However, there have been lots of new additions and changes to phytools since the last CRAN update, version 0.1-9. Here's a (probably non-comprehensive) sample:

1) A new function, threshBayes, to analyze the threshold model from quantitative genetics (described here and here).

2) A new function, bmPlot, for visualizing Brownian motion evolution (here) and evolution under the threshold model (here).

3) A bug fix in the important phytools function nodeHeights (here).

4) An update to the canonical correlation analysis function phyl.cca, as well as some other things (described here).

5) Renaming & movement of the likelihood function for joint Pagel's λ for multiple characters to the NAMESPACE (described, along with instructions for using this function - likMlambda - to perform joint estimation of Pagel's λ, here).

6) Minor update to phylogenetic principal components function, phyl.pca (here).

7) New version of fitDiversityModel with much improved likelihood estimation (described here).

8) New tree plotting method in the function fancyTree (here) along with various improvements thereof (here).

9) A number of changes to documentation and some minor fixes and updates to internal phytools functions.

The source build of phytools 0.2-0 is already available from CRAN (although will not yet have percolated through all mirrors), and from my phytools page. Obviously, it may take a few days for Mac OS and Windows binaries to be available.

read.simmap: computation time rises more than linearly with the number of trees

Sam Price reports that computation time for phytools function read.simmap(...,version=1.5) rises more than linearly with the number of trees in the input file. That means that it takes (substantially) more than twice as long to read in a data file containing 200 trees than it does to read in a data file containing 100 trees. (Or, equivalently, it would be faster to split your trees into two files and then read them in separately.)

Well, it turns out that Sam is right, of course - but I haven't yet figured out why. If I have time, I will try and figure this out tomorrow. Here is the result of reading in tree files with various numbers of trees using read.simmap and timing the result using system.time (if it seems slow it is because I am running it on my VAIO ultrabook):

Modification to fancyTree(...,type="droptip") so that it works on non-ultrametric trees

In response to a phytools user comment I gave a trick (here) for plotting the leaves and edges to be dropped in a different color & style from the rest of the tree. I then incorporated this trick into the phytools function fancyTree (described here). What I failed to point out was that this trick will not work on non-ultrametric trees. (To be honest, it didn't even really occur to me to point this out, so jaded am I from working with simulated data and trees.) Well, I have fixed this now, so that fancyTree(...,type="droptip") no longer uses the trick I described here, although it does still, generally speaking, use the method of fancyTree(...,type="extinction").

Users can check out the updated code on my R phylogenetics page here; or just download the latest (non-CRAN) version of phytools (v0.1-98: here), and install from source.

Let's check out how it works:

> # install new version of phytools, if necessary
> install.packages("phytools_0.1-98.tar.gz",type="source", repos=NULL)
> # simulate non-ultrametric tree & drop tips
> set.seed(10)
> tree<-rtree(n=20)
> pruned<-fancyTree(tree,"drop",tip=sample(tree$tip)[1:5])

Well, that works great. The only thing I'm not satisfied with is that if pruning changes the total length of the tree, then the second panel will be automatically rescaled so that corresponding nodes in the two plotted trees no long have corresponding horizontal positions. For example:

> set.seed(100)
> tree<-rtree(n=20)
> tips<-c("t6","t8","t14","t12","t4","t20")
> pruned<-fancyTree(tree,"drop",tip=tips)

Still working on this. . . .

Matching partial strings

To allow users to specify plot type (argument type) in the phytools function fancyTree without having to write the string for the desired type out in full, I yesterday wrote a custom function to match partial strings - the internal phytools function matchType. This is what it looked like:

matchType<-function(type,types){
 for(i in 1:length(types))
  if(all(strsplit(type,split="")[[1]]==strsplit(types[i], split="")[[1]][1:length(strsplit(type,split="")[[1]])]))
   type=types[i]
 return(type)
}

The very long if statement basically asks if (when type is split into a character vector with strsplit) all it elements (up to its total length) match all the elements of the current type under evaluation. So, for instance, we can do the following:

> type="extinct"
> matchType(type,c("extinction","traitgram3d","droptip"))
[1] "extinction"
> type="drop"
> matchType(type,c("extinction","traitgram3d","droptip"))
[1] "droptip"

or any number of a virtually infinite number of variations thereof.

It did strike me, though, that there might be a pre-existing R function to do this - and sure enough, there is: pmatch in the base package. It is slightly different from matchType in that it returns the index of the (partially) matching strings, but it is otherwise pretty much identical. To use it as a substitute for matchType, we would just have to do the following:

> types=c("extinction","traitgram3d","droptip")
> type="extinct"
> types[pmatch(type,types)]
[1] "extinction"
> type="drop"
> types[pmatch(type,types)]
[1] "droptip"

Sigh - oh well!

New option in fancyTree to plot tree pruning

I just added the code snippet from yesterday as a additional option to the phytools function fancyTree (previously described here, here, and here). A link to the source code for the new function is here.

Now, fancyTree, which is the function into which I'm dropping my efforts to make it easier for R phylogenetics users to plot various types of specialized phylogenies, has three different options (specified by the option type):

1) type="extinction": this method, described here, plots branches leading only to extinct taxa, and branches before the MRCA of all extant species in the tree, using red dashed lines (whereas all the other branches are plotted in black). For example:

> set.seed(10)
> tree<-rlineage(1,0.5,4)
> fancyTree(tree,type="extinction")

2) type="traitgram3d": this method, described here, uses phytools::phylomorphospace3d to plot a three dimensional 'traitgram'. For instance:

> tree<-pbtree(n=50,scale=10)
> Y<-sim.corrs(tree,vcv=matrix(c(1,0.75,0.75,1),2,2))
> fancyTree(tree,type="traitgram3d",X=Y, control=list(spin=FALSE))

3) Finally, the newest addition, type="droptip": this method takes the tree and set of tips to drop, and then plots the branches to be pruned in a separate color. Finally, it returns the pruned tree. For instance:

> tree<-pbtree(n=30)
> tips<-sample(tree$tip.label)[1:10]
> pruned<-fancyTree(tree,type="droptip",tip=tips)

I have added this new version of fancyTree to phytools. The new package version (v0.1-97), which can be downloaded and installed from source, is here.

> install.packages("phytools_0.1-97.tar.gz",type="source", repos=NULL)
Installing package(s) into ...
...
* DONE (phytools)
> library(phytools)

Plotting with the leaves and edges to be dropped highlighted

An anonymous reader commented the following:

I am wondering, is it possible to give different colours to the taxa after pruning by using the drop.tip. I just want to show the locations of the used species in the whole tree. All the taxa in tree retains but with different colours? (e.g. used species blue or red and pruned ones grey)

Ok, I'm going to interpret this as meaning - can we paint the leaves and edges that we're going to drop from the tree a different color, plot these, then drop them? (Once we've dropped these tips and edges from the tree, they are really gone - so we definitely can't plot them any more at that point.)

The answer is, of course, yes - and I'm going to give a quick trick for doing this using functions from the packages phytools (my package, of course) and ape.

Here I go. For the purposes of illustration, I will use a simulated, 30-taxon ultrametric tree; and a random set of ten tips.

> # first simulate a tree > tree<-pbtree(n=30) > # choose a random set of tips to prune > tips<-sample(tree$tip.label)[1:10] > tips [1] "t20" "t8" "t17" "t9" "t15" "t2" "t13" "t14" "t16" [10] "t11" > # find the node numbers for those tips > nodes<-sapply(tips,function(x) which(x==tree$tip.label)) > # and the rows in our edge matrix for all terminal branches > # to be dropped > edges<-sapply(nodes,function(x) which(x==tree$edge[,2])) > # now let's shorten them by an imperceptible amount > tree$edge.length[edges]<-tree$edge.length[edges]- 0.0001*max(nodeHeights(tree)) > par(mfrow=c(2,1)) > # plot the tree using fancyTree in phytools > fancyTree(tree) > # drop tips and plot the pruned tree > tree<-drop.tip(tree,tips) > fancyTree(tree)

Pretty cool! The way this works is by identifying and then imperceptibly shortening the terminal edges leading to species that we plan to drop from the tree. We can then use the phytools function fancyTree(...,type="extinction") to automatically detect these species (as well as any internal branches leading to these species) as "extinct." It then plots all these edges with red dashed lines, instead of the standard black. (We could modify these color scheme - but as of now we'd have to change the function code).

That's it! Thanks for the suggestion.