I just added a minor update to my ltt() function. The only significant difference between this function and plotLtt() in Rabosky's {laser} package, is that my function can also plot trees with non-contemporaneous (i.e., extinct) tips. The code for ltt() can be found on my R-phylogenetics page (direct link here).
Basically, in preparing for a talk that I will be giving tomorrow, I realized that my function relied on a "cladewise" order of the edges of the input "phylo" object. This is the order that is created by read.tree() (as well as my functions read.newick() and read.simmap()). However, some other function in R return trees with edges in "pruningwise" order, which is more efficient for post-order tree traversal and pruning. [For more information about this, type ?reorder.phylo at the prompt.]
Luckily, {ape} has a function that can switch from "pruningwise" to "cladewise" orders, or vice versa. So, the bug with ltt() was an easy fix. I simply added the line:
tree<-reorder.phylo(tree,order="cladewise")
to the code before beginning calculation.
[In addition to this change, I also moved "tol" - the tolerance threshold for rounding errors in node heights - to be an optional input instead of fixed. This probably won't affect most users.]
Hello Liam!
ReplyDeleteWell, it's rather slow on a big tree, like my graptolite supertree.
>system.time(ltt(grapt_mbl,plot=TRUE,log.lineages=F,tol=0.1))
user system elapsed
171.45 0.14 177.22
As it so happens, I was writing a function for the same thing just this week. Here's a simplified version of my code:
require(geiger)
ttree<-birthdeath.tree(0.1,0.05,time.stop=200,taxa.stop=100)
tblen<-0.01
ntime<-dist.nodes(ttree)[,Ntip(ttree)+1]
ntime<-max(ntime)-ntime
time<-seq(max(ntime),min(ntime),by=-tblen)
death<-ntime[1:Ntip(ttree)] birth<-ntime[(Ntip(ttree)+1):length(ntime)]
div<-sapply(time,function(x) 1+sum(birth>=x)-sum(death>=x))
plot(time,div,type="l",
xlim=c(max(time),min(time)),
xlab="Time (MY)",ylab="Lineage Richness")
This runs much faster than your function, mostly by replacing the for() and while() loops with a simple sapply() that counts the number of births and subtracts the number of deaths. Are the loops serving some other important use that I am not fully grokking?
-Dave
Hi Dave. Yes, I noticed that is is slow. I don't know why I did it with loops here. I don't think there is any particular reason. I will try your version. - Liam
ReplyDelete