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.

Simple new S3 plot method for the "fastAnc" object class

I just committed a new update to phytools that involved adding a new generic plot method for the "fastAnc" object class. This was pretty simple, with the main trick to it being that I modified the "fastAnc" object just by adding the tree and input data from fastAnc as object attributes.

Here’s a simple demo.

## load package
library(phytools)

## check version number
packageVersion("phytools")

## [1] '2.5.4'

## load tree and data
data(mammal.tree)
data(mammal.data)

## extract log body mass
lnBodyMass<-setNames(log(mammal.data$bodyMass),
  rownames(mammal.data))

## run ASR using phytools::fastAnc
mammal.anc<-fastAnc(mammal.tree,lnBodyMass)
mammal.anc

## Ancestral character estimates using fastAnc:
##       50       51       52       53       54       55       56       57       58       59 
## 4.616864 3.928453 3.570720 3.372081 5.008271 5.416959 2.503895 1.783068 2.074242 2.092048 
##       60       61       62       63       64       65       66       67       68       69 
## 2.102201 2.427742 2.879777 2.935790 4.000602 3.660331 4.314389 4.607878 4.760248 4.846334 
##       70       71       72       73       74       75       76       77       78       79 
## 5.299692 5.545168 7.076476 5.516430 5.552013 5.123064 5.352999 5.190895 5.148764 5.146091 
##       80       81       82       83       84       85       86       87       88       89 
## 5.221873 6.105385 4.596915 3.642387 3.618018 3.586540 4.613914 4.591181 4.847749 4.777965 
##       90       91       92       93       94       95       96       97 
## 4.873585 4.882390 5.109305 4.960400 4.878856 4.928119 4.832228 4.246744

plot(mammal.anc,ftype="i",fsize=0.6,offset=0.5,
  title="log(body mass)")

Cool.

Let’s simulate a large tree and then try it in a different style.

## simulate tree & data
sim_tree<-pbtree(n=500,scale=100)
sim_data<-fastBM(sim_tree)
## do ASR
sim_asr<-fastAnc(sim_tree,sim_data,CI=TRUE)
## print ASR
print(sim_asr,printlen=4)

## Ancestral character estimates using fastAnc:
##        501       502       503       504     
##  -1.968242 -2.989587 -4.870655 -9.828655 ....
## 
## Lower & upper 95% CIs:
##          lower     upper
## 501   -8.59862  4.662136
## 502  -8.817602  2.838427
## 503 -10.380942  0.639632
## 504 -15.847845 -3.809464
##           ....      ....

## plot results
plot(sim_asr,type="fan",ftype="off",lwd=1,color="grey",
  node.cex=0.8,tip.cex=0.5,ylim=c(-120,100),legend="bottomleft")

That’s pretty much the idea.

User-supplied bug fix for fastAnc

A phytools user, David Labonte from the University of Cambridge, recently reported the following bug with the phytools function fastAnc:

“I am using your fastAnc function to estimate ancestral states of a continuous variable, and it runs smoothly, however only if CI=FALSE. I also noted that it returns a longer variance than ancestral state estimation vector. Notably, this problem only arises for non-dichotomous trees.”

David also was kind enough to provide the following example which reproduces the error very nicely (modified slightly):

library(phytools)
set.seed(1)
tree<-pbtree(n=10,scale=1)
plotTree(tree,node.numbers=TRUE)

x<-fastBM(tree)
fastAnc(tree,x,vars=TRUE,CI=TRUE) ## runs smoothly without error

## $ace
##          11          12          13          14          15          16 
## -0.35538758 -0.39739912 -0.32825378 -0.51263348 -0.32987276 -0.13938453 
##          17          18          19 
## -0.14780235 -0.07371257 -0.09712277 
## 
## $var
##          11          12          13          14          15          16 
## 0.081863455 0.063948796 0.045742058 0.019817948 0.006889495 0.027865560 
##          17          18          19 
## 0.003468311 0.022557094 0.009799787 
## 
## $CI95
##          [,1]        [,2]
## 11 -0.9161787  0.20540351
## 12 -0.8930459  0.09824762
## 13 -0.7474467  0.09093912
## 14 -0.7885549 -0.23671206
## 15 -0.4925586 -0.16718692
## 16 -0.4665670  0.18779789
## 17 -0.2632314 -0.03237331
## 18 -0.3680854  0.22066021
## 19 -0.2911508  0.09690522

tree<-collapse.to.star(tree,fastMRCA(tree,"t6","t9"))
plotTree(tree,node.numbers=TRUE)

fastAnc(tree,x) ## works no problem

##         11         12         13         14         15         16 
## -0.3498724 -0.3911156 -0.3188186 -0.5111098 -0.3295891 -0.1216447

fastAnc(tree,x,vars=TRUE) ## works, but vars is wrong length

## $ace
##         11         12         13         14         15         16 
## -0.3498724 -0.3911156 -0.3188186 -0.5111098 -0.3295891 -0.1216447 
## 
## $var
##          11          12          13          14          15          16 
## 0.077702044 0.060153127 0.040920299 0.018985301 0.006625218 0.015899879

fastAnc(tree,x,vars=TRUE,CI=TRUE) ## doesn't work at all

## $ace
##         11         12         13         14         15         16 
## -0.3498724 -0.3911156 -0.3188186 -0.5111098 -0.3295891 -0.1216447 
## 
## $var
##          11          12          13          14          15          16 
## 0.077702044 0.060153127 0.040920299 0.018985301 0.006625218 0.015899879 
## 
## $CI95
##          [,1]        [,2]
## 11 -0.8962241  0.19647933
## 12 -0.8718278  0.08959662
## 13 -0.7153024  0.07766520
## 14 -0.7811727 -0.24104701
## 15 -0.4891242 -0.17005403
## 16 -0.3687903  0.12550100

Even better, David solved the bug by correctly identifying the error in the code. In his own words, he says:

“I looked through the code of the function, and while I am certainly not a R coding expert, I believe the problem lies in line 28:

27    if (vars || CI) {
28        v[as.character(ancNames[, 2])]
29        names(v) <- ancNames[, 1]
30    }

analogous to the previous lines, I think this should read

27    if (vars || CI) {
28        v <- v[as.character(ancNames[, 2])]
29        names(v) <- ancNames[, 1]
30    }

This is exactly correct, and if we fix this then we find that the function now works perfectly:

source("fastAnc.R")
fastAnc(tree,x,vars=TRUE,CI=TRUE)

## $ace
##         11         12         13         14         15         16 
## -0.3498724 -0.3911156 -0.3188186 -0.5111098 -0.3295891 -0.1216447 
## 
## $var
##          11          12          13          14          15          16 
## 0.077702044 0.060153127 0.040920299 0.018985301 0.006625218 0.015899879 
## 
## $CI95
##          [,1]        [,2]
## 11 -0.8962241  0.19647933
## 12 -0.8718278  0.08959662
## 13 -0.7153024  0.07766520
## 14 -0.7811727 -0.24104701
## 15 -0.4891242 -0.17005403
## 16 -0.3687903  0.12550100

Cool! If only fixing all the bugs in phytools was this easy!

Fast ML estimation of ancestral states for a continuously valued trait

I just wrote a new function for ancestral character estimation that takes advantage of the fact that the ancestral value for the root node computed during the calculation of independent contrasts is also the MLE of the root. Re-rooting at all internal nodes and the recomputing the PIC root state, cumbersome as it sounds, is actually much faster than finding ancestral states via numerical optimization of the likelihood. The main part of the code for this function is therefore as follows:

M<-btree$Nnode
N<-length(btree$tip)
anc<-vector()
for(i in 1:M+N){
  anc[i-N]<-ace(x,multi2di(root(btree,node=i)), method="pic")$ace[1]
  names(anc)[i-N]<-i
}

Obviously, although this function is designed as a faster version of ace(...,type="continuous",method="ML"), the function actually runs by using many calls to ace(...,method="pic")!!

The object name btree in this code is used to denote binary tree. This highlights the fact that, obviously, contrasts can only be computed for bifurcating trees. We would still like our function to run, though, if the tree contains multifurcating. The remainder of the function code is dedicated to, first, changing a multifurcating tree to a bifurcating one; and then, after computing ancestral states using the code given above, lining up the nodes on the binary tree to the nodes of the original tree containing multifurcations.

The way a did this is a little ad hoc, but it seems to work. I basically went through all the nodes on both the binary and multifurcating tree and, for each node, pulled out a list of all the descendant tips. Then, I used these lists of descendant tips to match nodes between trees.

The code for this function, fastAnc, is here. Let's try it and compare to both ape::ace and phytools::anc.ML.

> library(phytools)
> # first, simulate bifurcating tree & data
> tree<-rtree(200)
> x<-fastBM(tree)
> # now load the source & estimate states with each function
> source("fastAnc.R")
> system.time(res1<-ace(x,tree,method="ML")$ace)
  user  system elapsed
 19.69    0.05   20.38
> system.time(res2<-anc.ML(tree,x,maxit=10000)$ace)
  user  system elapsed
236.67    0.34  255.00
> system.time(res3<-fastAnc(tree,x))
  user  system elapsed
  2.60    0.00    2.98
> # plot to compare
> par(mfrow=c(2,1),mai=c(1.02,0.82,0.1,0.1))
> plot(res1,res3,xlab="ace",ylab="fastAnc")
> plot(res2,res3,xlab="anc.ML",ylab="fastAnc")

Obviously, we get the same estimates from each function - but at much greater computational cost (particularly for anc.ML, which is pretty terrible).

That was for a fully bifurcating tree. We can also do a smaller example to make sure we are back-translating our node IDs correctly. Let's try it:

> tree<-rtree(12)
> # collapse two shortest branches into multifurcations
> tree<-di2multi(tree,tol= sort(tree$edge.length[tree$edge[,2]>length(tree$tip)])[3])
> plot(tree,no.margin=TRUE); nodelabels()

> # simulate data
> x<-fastBM(tree)
> # estimate ancestral states using the three methods
> res1<-ace(x,tree,method="ML")$ace
Error in ace(x, tree, method = "ML") :
 "phy" is not rooted AND fully dichotomous.
> res2<-anc.ML(tree,x)$ace
> res3<-fastAnc(tree,x)
> par(mai=c(1.02,0.82,0.2,0.2))
> plot(res2,res3,xlab="anc.ML",ylab="fastAnc")
> res2
        13          14          15          16          17
0.28663746  0.05138613 -0.25089177  0.50494223  0.13529718
        18          19          20          21
0.03920500  0.49088849  1.24192050  1.63727735
> res3
        13          14          15          16          17
0.28665585  0.05139666 -0.25089395  0.50494620  0.13528936
        18          19          20          21
0.03920408  0.49088277  1.24191236  1.63727663

First off - ace doesn't work at all if the tree is not bifurcating (not sure why this is). anc.ML works fine, but will be very slow for large trees (as we've discovered). fastAnc gives the same estimates and node names - even though it had to first convert to a bifurcating tree, and then back-translate the node numbers. Great!

New version of fastAnc; new build of phytools

I just posted a new version of the function fastAnc for (relatively) fast ancestral character estimation. The function is previously described here. The main addition to this new version is that now the function (optionally) computes the variance on the ancestral state estimates based on equation (6) of Rohlf (2001), as well as (optionally) 95% confidence intervals on the states. The updated code is here; however, I have also posted a new build of phytools - which can be downloaded here and installed from source.

Note that equation (6) of Rohlf (2001) only gives the relative variance on the ancestral state estimate at the root node. To scale that estimate to our data, we need to multiply by the phylogenetic variance for our continuous trait. This can be computed as the mean square of the contrasts. Once we have the variances, we can compute our 95% CIs on the estimates as the estimates +/- 1.96 × the square root of the variances.

I didn't realize this when I was writing the function, but it turns out to be the case that this update to fastAnc depends on ape >= 3.0-7 (i.e., the newest version of ape as of the date of writing). This is because the options in the ape function for independent contrasts, pic, were expanded in the latest release to include the option of returning the tree with branches scaled to expected variance - which we can conveniently exploit to do the calculation of equation (6) in Rohlf.

I should also point out that the 95% CIs obtained by this function differ in a substantial way from the 95% CIs computed in ace. Specifically, the 95% CIs computed in ace would seem to be too small. We can show this relatively easily by simulation, as follows:

> onCI.ace<-onCI.fastAnc<-vector()
> N<-100
> for(i in 1:1000){
+ tree<-pbtree(n=N)
+ x<-fastBM(tree,internal=TRUE)
+ a<-fastAnc(tree,x[1:N],vars=TRUE,CI=TRUE)
+ onCI.fastAnc[i]<-sum((x[1:tree$Nnode+N]>a$CI95[,1])*(x[1:tree$Nnode+N]< a$CI95[,2]))/tree$Nnode
+ b<-ace(x[1:N],tree,CI=TRUE)
+ onCI.ace[i]<-sum((x[1:tree$Nnode+N]>b$CI95[,1])*(x[1:tree$Nnode+N]< b$CI95[,2]))/tree$Nnode
+ }
There were 24 warnings (use warnings() to see them)
> warnings()
Warning messages:
1: In sqrt(diag(solve(h))) : NaNs produced
2: In sqrt(diag(solve(h))) : NaNs produced
3: ...
> # this should be 0.95
> mean(onCI.fastAnc)
[1] 0.9483737
> mean(onCI.ace,na.rm=TRUE)
[1] 0.6738759

This simulation shows that although our 95% CIs computed in fastAnc include the generating values almost exactly 95% of the time (94.8% across 1000 simulations with 99 estimated ancestral states per simulation), ace CIs only include the generating value about 67% of the time.

I'm not exactly sure why this is the case, but my best guess is based on the warnings which tell us that the Hessian is being used to compute the standard errors of the parameter estimates and thus the CIs. This is an asymptotic property of the likelihood surface, and for finite sample this approximation can be quite bad (as we see above).