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> 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
> 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>> 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>
>>
>> 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
>> 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
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
