## 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.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.000035