Addendum to adding a single tip to the tree

Yesterday I posted an extremely simple function for adding an extra tip to the tree. In response to a user request, here is an addendum in which the function also checks if the tree is ultrametric and (so long as a custom branch length is not supplied) computes a branch length for the added tip such that the tree remains ultrametric after the new tip is added. Since this wraps around the 'ape' function bind.tree, I have also added the argument position, which is the distance below the specified node where the new tip should be added. First, here's the function:

bind.tip<-function(tree,tip.label,edge.length=NULL, where=NULL,position=0){
  if(is.null(where)) where<-length(tree$tip)+1
  if(is.null(edge.length)&&is.ultrametric(tree)){
    H<-nodeHeights(tree)
    if(where==(length(tree$tip)+1))
     edge.length<-max(H)
    else
     edge.length<-max(H)-H[tree$edge[,2]==where,2]+position
  }
  tip<-list(edge=matrix(c(2,1),1,2),
    tip.label=tip.label,
    edge.length=edge.length,
    Nnode=1)
    class(tip)<-"phylo"
  obj<-bind.tree(tree,tip,where=where,position=position)
  return(obj)
}

And here's a demo:

> tree<-pbtree(n=10)
> plotTree(tree,node.numbers=T)

> # now add 1/2 branch length below node 16
> tree2<-bind.tip(tree,"t11",where=16, position=0.5*tree$edge.length[which(tree$edge[,2]==16)])
> plotTree(tree2,node.numbers=T)

Cool - it works.

Note that if we want to add a new terminal edge along a branch leading to a tip - we have to specify the tip by node number not tip label. So, for instance, to add another tip to our new tree halfway along the branch leading to tip t11, we would do:

> tree3<-bind.tip(tree2,"t12",
 where=which(tree2$tip.label=="t11"),
 position=0.5*tree2$edge.length[which(tree2$edge[,2]==
 which(tree2$tip.label=="t11"))])
> plotTree(tree3,node.numbers=T)