Monday, December 19, 2011

Faster way to compute the maximum value of Pagel's lambda

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