It has previously been pointed out that the maximum value of λ (such that the likelihood is still defined) for an ultrametric phylogeny can be computed as the maximum of the diagonal of vcv.phylo(tree) (or the mean, as this is a costant), divided by the maximum of the off-diagonals. In other words, we could write a function:
maxLambda<-function(tree){
C<-vcv(tree)
max(C)/max(C[lower.tri(C)])
}
Now, this function is not particularly speedy because the computation vcv.phylo(tree) is quite slow. For instance:
> tree<-drop.tip(birthdeath.tree(b=1,d=0,taxa.stop=2001),"2001")
> maxLambda<-function(tree){
+ C<-vcv(tree)
+ max(C)/max(C[lower.tri(C)])
+ }
> system.time(lambda<-maxLambda(tree))
user system elapsed
33.65 0.20 33.93
> lambda
[1] 1.000035
Well, I just realized we could do this much more efficiently using the function nodeHeights from the phytools package. Here, the function would look something like this:
maxLambda<-function(tree){
H<-nodeHeights(tree)
max(H[,2])/max(H[,1])
}
And, indeed, this is much faster:
> system.time(lambda<-maxLambda(tree))
user system elapsed
0.66 0.00 0.65
> lambda
[1] 1.000035
No comments:
Post a Comment
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.