ML estimation of ancestral states for large trees using ace

A phytools user recently contacted me to report differences between several different functions in R to estimate ancestral states at internal nodes for continuous characters using likelihood.

Since I had no idea to expect this, I ran a quick comparison of ace(method="ML"), fastAnc, and anc.ML. Ostensibly, all three methods are performing ML estimation assuming a Brownian model of evolutionary change. The only difference is the method of optimization. fastAnc, for example, using a 're-rooting method' in which the tree is re-rooted at every internal node and the contrasts algorithm is used to obtain the ML state for that node. anc.ML, also in phytools, uses numerical optimization - but initiates the search using the values from fastAnc.

It occurred to me immediately that the problem may be related not to the implementation of the model, but rather to the simple issue of numerical optimization. Here is a demo suggesting that this could indeed be the case:

library(phytools)

## tiny example:
tree<-pbtree(n=26,tip.label=LETTERS)
x<-fastBM(tree)
fit1<-ace(x,tree,type="continuous",method="ML")
fit2<-fastAnc(tree,x,vars=TRUE,CI=TRUE)
fit3<-anc.ML(tree,x)
obj<-cbind(fit1$ace,fit2$ace,fit3$ace)
colnames(obj)<-c("ace(method=\"ML\")","fastAnc","anc.ML")
pairs(obj,pch=21,bg="grey",cex=1.5)

## small example
tree<-pbtree(n=100)
x<-fastBM(tree)
fit1<-ace(x,tree,type="continuous",method="ML")
fit2<-fastAnc(tree,x,vars=TRUE,CI=TRUE)
fit3<-anc.ML(tree,x)
obj<-cbind(fit1$ace,fit2$ace,fit3$ace)
colnames(obj)<-c("ace(method=\"ML\")","fastAnc","anc.ML")
pairs(obj,pch=21,bg="grey",cex=1.5)

## medium example
tree<-pbtree(n=200)
x<-fastBM(tree)
fit1<-ace(x,tree,type="continuous",method="ML")
fit2<-fastAnc(tree,x,vars=TRUE,CI=TRUE)
fit3<-anc.ML(tree,x)
obj<-cbind(fit1$ace,fit2$ace,fit3$ace)
colnames(obj)<-c("ace(method=\"ML\")","fastAnc","anc.ML")
pairs(obj,pch=21,bg="grey",cex=1.5)

## large example:
tree<-pbtree(n=1000)
x<-fastBM(tree)
fit1<-ace(x,tree,type="continuous",method="ML")

## Warning in sqrt(diag(solve(h))): NaNs produced

fit2<-fastAnc(tree,x,vars=TRUE,CI=TRUE)
fit3<-anc.ML(tree,x)
obj<-cbind(fit1$ace,fit2$ace,fit3$ace)
colnames(obj)<-c("ace(method=\"ML\")","fastAnc","anc.ML")
pairs(obj,pch=21,bg="grey",cex=1.5)

What we should see is that the results are identical for 'tiny' and relatively modest-sized trees, but go awry for large & possibly medium-sized trees where the number of parameters (the states at all the nodes) really explodes. This might suggest that the problem is with numerical optimization on large trees, rather than an issue with how the model is implemented or how the likelihood computed.

This should serve as a reminder that any method depending on numerical optimization is only as good as its optimization routine.

Bug fix in fastBM for non "clade wise" ordered trees

I just pushed a small fix to fastBM for cases in which the input trees are not in "cladewise" order. Since re-ordering is very fast, the fix first checks to see if the current order is "cladewise" (or if no order attribute is specified), and if not, reorders the tree's edges in a clade-wise fashion.

This bug was identified by Josselin Cornuault. Thanks!

Function to set custom tip-spacing for subtrees of a plotted tree now in phytools (along with custom plot and print methods)

Earlier today I posted about a general method to expand or contract the tip-spacing of all the taxa descended from a particular node or set of nodes.

Now, this function has been added to the phytools package, along with a pair of S3 methods to print & plot the computed object.

The plot method I hadn't thought of earlier. All it does is use do.call internally to send any arguments we give the method to plotTree or plotSimmap internally, depending on the tree object class. This means that any argument of plotTree or plotSimmap can also be passed to the S3 plot method for the object class "expand.clade".

For instance:

library(phytools)
packageVersion("phytools")

## [1] '0.6.15'

obj

## An object of class "expand.clade" consisting of:
## (1) A phylogenetic tree (x$tree) with 100 tips and
##     99 internal nodes.
## (2) A vector (x$tips) containing the desired tip-spacing.

obj$tree

## 
## Phylogenetic tree with 100 tips and 99 internal nodes.
## 
## Tip labels:
##  t4, t30, t47, t86, t87, t97, ...
## 
## The tree includes a mapped, 2-state discrete character with states:
##  a, b
## 
## Rooted; includes branch lengths.

colors

##      a      b 
## "blue"  "red"

plot(obj,colors=colors,ftype="off",direction="downwards")

That's all for now.

General function for expanding (or contracting) the tip-spacing of subtrees during visualization

A couple of days ago a phytools user asked me about using custom tip-spacing in a fan-style tree plot in R. I posted a fairly straightforward solution; however how it works is not super obvious, and it involves a lot of scripting, so I thought it might be worth adding this as a function in phytools.

The idea with the function is not to plot a tree with custom tip-spacing - we can already do that using plotTree or plotSimmap. Rather, the function reorders the tree (if it is not in “cladewise” order), computes the desired spacing of each specified subtree, and then returns an object consisting of the tree & a vector that we can use in our tree plot.

Here's the function:

expand.clade<-function(tree,node,factor=5){
    cw<-reorder(tree)
    tips<-setNames(rep(1,Ntip(tree)),cw$tip.label)
    get.tips<-function(node,tree){
            dd<-getDescendants(tree,node)
            tree$tip.label[dd[dd<=Ntip(tree)]]
    }
    desc<-unlist(lapply(node,get.tips,tree=cw))
    for(i in 2:Ntip(cw)){
        tips[i]<-tips[i-1]+
            if(names(tips)[i]%in%desc){
                1 
            } else if(names(tips)[i-1]%in%desc){
                1
            } else 1/factor
    }
    obj<-list(tree=tree,tips=tips)
    class(obj)<-"expand.clade"
    obj
}

## S3 method for the object class
print.expand.clade<-function(x,...){
    cat("An object of class \"expand.clade\" consisting of:\n")
    cat(paste("(1) A phylogenetic tree (x$tree) with",Ntip(x$tree),
        "tips and\n   ",obj$tree$Nnode,"internal nodes.\n"))
    cat("(2) A vector (x$tips) containing the desired tip-spacing.\n\n")
}

Here's an example:

tree

## 
## Phylogenetic tree with 100 tips and 99 internal nodes.
## 
## Tip labels:
##  t4, t30, t47, t86, t87, t97, ...
## 
## The tree includes a mapped, 2-state discrete character with states:
##  a, b
## 
## Rooted; includes branch lengths.

nodes<-c(108,140,168)
for(i in 1:length(nodes)) tree<-paintSubTree(tree,nodes[i],"b","a")
obj<-expand.clade(tree,nodes,factor=2.5)
obj

## An object of class "expand.clade" consisting of:
## (1) A phylogenetic tree (x$tree) with 100 tips and
##     99 internal nodes.
## (2) A vector (x$tips) containing the desired tip-spacing.

plot(obj$tree,colors=setNames(c("blue","red"),c("a","b")),
    tips=obj$tips,type="fan",fsize=0.7,lwd=3)

It also works for other tree styles. For instance:

obj<-expand.clade(tree,nodes,factor=3)
obj

## An object of class "expand.clade" consisting of:
## (1) A phylogenetic tree (x$tree) with 100 tips and
##     99 internal nodes.
## (2) A vector (x$tips) containing the desired tip-spacing.

plot(obj$tree,colors=setNames(c("blue","red"),c("a","b")),
    tips=obj$tips,fsize=0.4,lwd=1,ftype="off",xlim=c(0,6.5))
foo<-function(node,tree) extract.clade(tree,node)$tip.label
tips<-unlist(lapply(nodes,foo,tree=tree))
nn<-sapply(tips,function(tip,tree) which(tree$tip.label==tip),tree=tree)
linklabels(tips,nn,link.type="curved",cex=0.6)

To get the final example to work, I had to push a fix to the function linklabels. The fix is already on GitHub & the functions of this post will be added soon.

That's it.