Sampling edge lengths under a Yule process

There has been a little bit of discussion today on R-sig-phylo listserve about transforming branch lengths.

One thing that wasn't mentioned was the possibility of sampling branch lengths under a model. I thought it would be straightforward to sample branch lengths under a Yule model (that is, a pure-birth speciation model).

The following is code that does this. (Set plot=TRUE to see a cool animation of the tree being 'grown' from left to right.)

pb_edgelength<-function(tree,b=1,plot=TRUE,...){
    ll<-rexp(n=Ntip(tree)-1,rate=2:Ntip(tree)*b)
    tree$edge.length<-rep(0,nrow(tree$edge))
    live.nodes<-Descendants(tree,Ntip(tree)+1,"children")
    tips<-vector()
    for(i in 1:length(ll)){
        tips<-c(tips,live.nodes[live.nodes<=Ntip(tree)])
        live.nodes<-setdiff(live.nodes,tips)
        ii<-which(tree$edge[,2]%in%c(live.nodes,tips))
        tree$edge.length[ii]<-tree$edge.length[ii]+ll[i]
        node<-if(length(live.nodes)<=1) live.nodes else 
            sample(live.nodes,1) ## choose one node
        live.nodes<-c(setdiff(live.nodes,node),
            Descendants(tree,node,"children"))
        if(plot) plotTree(tree,...)
    }
    tree
}

OK, now let's try it out with a tree obtained using rtree from the ape package:

library(phytools)
library(phangorn)
tree<-rtree(n=100,br=NULL) ## no branch lengths
t.pb<-pb_edgelength(tree,plot=FALSE)
plotTree(t.pb,ftype="off")

par(mar=c(5.1,4.1,2.1,2.1))
obj<-ltt(t.pb)

obj$gamma

## [1] -0.1807099

Compare this to Grafen's edge lengths from compute.brlen:

t.grafen<-compute.brlen(tree)
plotTree(t.grafen,ftype="off")

par(mar=c(5.1,4.1,2.1,2.1))
obj<-ltt(t.grafen)

obj$gamma

## [1] 6.732706

Obviously, the branching times from Grafen's branch length transformation are very different from those obtained under a Yule process!

That's it, really.