Hi Rob.
Your function seems to work fine. The reason ltt() returns a vectors
(typically) 2 elements longer is because, by design, my vectors are
flanked by time=0 and lineages=0 at the start of the tree; and
time=tree.length and lineages=num.taxa at the end of the tree. If you
add these to either end, both functions are identical (except that yours
is much faster, of course).
I also noticed, though, that with your function, Rob, if you have a
speciation event exactly at the present, it will not be recorded (in
which case the number of elements returned by ltt() versus lttplot()
differ by 3 instead of 2). I'm sure that this can be fixed easily.
Thanks for sharing this.
Best, Liam
--
Liam J. Revell
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://phytools.blogspot.com
On 8/11/2011 7:52 PM, Rob Lanfear wrote:
Hi Liam, Emmanuel, Nick, et al
Thanks for all the help.
Liam, I can confirm that your new function works fine on all the trees
I've tested (about 20). Thanks!
Some of my trees are pretty large (a few thousand taxa), so for
interest's sake I coded up Emmanuel's suggestion too.
The code is below (liberally borrowed from suggestion on this list, and
from Liam's code). For reasons I don't entirely understand, my code
gives my an ltt vector which is 2 shorter than Liam's. Part of the
difference is that I start with 2 taxa at the root node, so I count one
less event than Liam at the start of the vector. But I'm not sure where
the other difference lies - I've probably made a mistake somewhere.
Any clues gratefully received.
Rob
lttplot<- function(phy, plot=TRUE, log.lineages=FALSE){
tol<-1e-6
x<- dist.nodes(phy)
n<- Ntip(phy)
ROOT <- n+1
#distances from root to all nodes and tips
events<- x[ROOT,]
#tips are extinctions (we remove the extant tips later)
names(events)[1:n] <- -1
#nodes are speciations
names(events)[-(1:n)] <- +1
#but extant tips don't count
present<- max(events)
events<- events[-which(events>present-tol)]
#sort w.r.t. time
events<- sort(events)
#the LTT plot is the cumsum of the names
ltt<- cumsum(names(events))
#add 1 to LTT, because the root node actually represents the creation of
2 taxa
#and we don't know the timing of the creation of the first taxon (root
branch)
ltt <- ltt+1
if(plot==TRUE){
if(log.lineages==TRUE)
plot(events,log(ltt),"s",xlab="time",ylab="log(lineages)")
else
plot(events,ltt,"s",xlab="time",ylab="lineages")
}
return(list(ltt=ltt, times=events))
}
On 12 August 2011 07:44, Liam J. Revell <liam.rev...@umb.edu
<mailto:liam.rev...@umb.edu>> wrote:
Hi Rob.
Thanks for identifying the bug. This should be fixed, I hope, in
the newest version of "phytools" (v0.0-6; available:
http://anolis.oeb.harvard.edu/__~liam/R-phylogenetics/
<http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/>).
The way this function works is by first computing the heights above
the root of all the nodes (including tip nodes) in the tree, and
then at each event (that is, a speciation or extinction), it counts
the number of edges that include that time point. This is slow. It
seems likely that Emmanuel's suggestion (in a previous email) to use
dist.nodes() would run much faster than this.
- Liam
--
Liam J. Revell
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.__revell/
<http://faculty.umb.edu/liam.revell/>
email: liam.rev...@umb.edu <mailto:liam.rev...@umb.edu>
blog: http://phytools.blogspot.com
On 8/11/2011 12:42 AM, Rob Lanfear wrote:
Hi Liam,
Thanks for the help. Extremely useful.
And thanks for clearing up my lack of understanding about the
differences in maximum ages!
Cheers,
Rob
On 11 August 2011 14:05, Liam J. Revell <liam.rev...@umb.edu
<mailto:liam.rev...@umb.edu>
<mailto:liam.rev...@umb.edu <mailto:liam.rev...@umb.edu>>> wrote:
Hi Rob.
I can reproduce your error, but I haven't figured out the
problem yet.
You can try an earlier version of this function, which seems
to work:
source("http://anolis.oeb.__ha__rvard.edu/~liam/R-____phylogenetics/ltt/v0.3/ltt.R
<http://harvard.edu/~liam/R-__phylogenetics/ltt/v0.3/ltt.R>
<http://anolis.oeb.harvard.__edu/~liam/R-phylogenetics/ltt/__v0.3/ltt.R
<http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/ltt/v0.3/ltt.R>>")
p2 <- ltt(t1, log.lineages=FALSE, drop.extinct=FALSE)
Sorry about this.
Also:
max(p1$times)==max(p2$times)
can be FALSE because if drop.extinct is set to TRUE, then
the crown
age of the pruned tree can be smaller than in the full tree
if some
lineages arising at the root of the tree do not leave any extant
descendants.
- Liam
--
Liam J. Revell
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.____revell/
<http://faculty.umb.edu/liam.__revell/>
<http://faculty.umb.edu/liam.__revell/
<http://faculty.umb.edu/liam.revell/>>
email: liam.rev...@umb.edu <mailto:liam.rev...@umb.edu>
<mailto:liam.rev...@umb.edu <mailto:liam.rev...@umb.edu>>
blog: http://phytools.blogspot.com
On 8/10/2011 9:50 PM, Rob Lanfear wrote:
library(TreeSim)
library(phytools)
library(geiger)
#simulate tree of 100 taxa with initial diversification
followed
by a
period of B=D
t1 <- sim.rateshift.taxa(100, 1, c(0.2, 0.2), c(0.2, 0.05),
c(1,1), c(0,
20))[[1]]
t1
#confirm number of extant taxa is 100
n.extant<- length(prune.extinct.taxa(t1)$____tip.label)
n.extant
#do ltt plot without extinct taxa (works fine)
p1<- ltt(t1, log.lineages=FALSE, drop.extinct=TRUE)
#do ltt plot with extinct taxa (looks odd)
p2<- ltt(t1, log.lineages=FALSE, drop.extinct=FALSE)
#max times don't seem to correspond between the two plots.
max(p1$times)==max(p2$times)
--
Rob Lanfear
Postdoc,
Centre for Macroevolution and Macroecology,
Research School of Biology,
Australian National University
Tel: +61 2 6125 7270 <tel:%2B61%202%206125%207270>
www.robertlanfear.com <http://www.robertlanfear.com>
<http://www.robertlanfear.com>
--
Rob Lanfear
Postdoc,
Centre for Macroevolution and Macroecology,
Research School of Biology,
Australian National University
Tel: +61 2 6125 7270
www.robertlanfear.com <http://www.robertlanfear.com>
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
