New message in force.ultrametric: this function is not a formal method to ultrametricize your tree!

I just updated the phytools function force.ultrametric to print the following message (unless explicitly turned off by the user by setting the optional argument message=FALSE):

## ***************************************************************
## *                          Note:                              *
## *    force.ultrametric does not include a formal method to    *
## *  ultrametricize a tree & should only be used to coerce to   *
## *   a phylogeny that fails is.ultramtric due to rounding --   *
## *   not as a substitute for formal rate-smoothing methods.    *
## ***************************************************************

The reason that I decided that this warning was necessary because I've seen force.ultrametric used in some papers in lieu of a formal rate-smoothing method for generating an ultrametric tree.

This is not appropriate.

What is force.ultrametric for then? It's essentially a coercion method that forces a tree that should be ultrametric, but that fails is.ultrametric (perhaps, for example, due to numerical precision or rounding of the edge lengths when it was written to file), to be so.

Here's what I mean.

Let's load the data object anoletree that comes loaded with phytools.

library(phytools)
data(anoletree)
anoletree<-as.phylo(anoletree)
plotTree(anoletree,ftype="i",fsize=0.7,lwd=1)

Looks pretty ultrametric, right? It should be.

Think again, though – it's not!

is.ultrametric(anoletree)

## [1] FALSE

What's wrong? We can see if we pull out the total height of all the tips of the tree as follows:

h<-setNames(sapply(1:Ntip(anoletree),
    nodeheight,tree=anoletree),
    anoletree$tip.label)
h

##            ahli         allogus     rubribarbus           imias          sagrei         bremeri quadriocellifer      ophiolepis 
##               6               6               6               6               6               6               6               6 
##         mestrei           jubar      homolechis        confusus           guafe         garmani        opalinus         grahami 
##               6               6               6               6               6               6               6               6 
##     valencienni      lineatopus       evermanni       stratulus           krugi      pulchellus       gundlachi       poncensis 
##               6               6               6               6               6               6               6               6 
##           cooki    cristatellus    brevirostris        caudalis          marron        websteri       distichus         alumina 
##               6               6               6               6               6               6               6               6 
##    semilineatus         olssoni       insolitus       whitemani       haetianus        breslini         armouri         cybotes 
##               6               6               6               6               6               6               6               6 
##         shrevei   longitibialis         strahmi        marcanoi        baleatus       barahonae        ricordii         cuvieri 
##               6               6               6               6               6               6               6               6 
##   altitudinalis        oporinus        isolepis        allisoni        porcatus        loysiana         guazuma        placidus 
##               6               6               6               6               6               6               6               6 
##        sheplani         alayoni     angusticeps        paternus       alutaceus    inexpectatus       clivicola    cupeyalensis 
##               6               6               6               6               6               6               6               6 
##    cyanopleurus         alfaroi      macilentus       vanidicus        baracoae          noblei      smallwoodi    luteogularis 
##               6               6               6               6               6               6               6               6 
##       equestris   bahorucoensis dolichocephalus      hendersoni     darlingtoni        aliniger      singularis    chlorocyanus 
##               6               6               6               6               6               6               6               6 
##     coelestinus        occultus 
##               6               6

Ok…. All the tree heights seem the same. But wait a sec….

options(digits=10)
h

##            ahli         allogus     rubribarbus           imias          sagrei         bremeri quadriocellifer      ophiolepis 
##      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994 
##         mestrei           jubar      homolechis        confusus           guafe         garmani        opalinus         grahami 
##      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994 
##     valencienni      lineatopus       evermanni       stratulus           krugi      pulchellus       gundlachi       poncensis 
##      5.99999994      5.99999994      5.99999994      5.99999994      6.00000000      6.00000000      5.99999994      5.99999994 
##           cooki    cristatellus    brevirostris        caudalis          marron        websteri       distichus         alumina 
##      5.99999994      5.99999994      5.99999994      5.99999988      5.99999988      5.99999994      5.99999994      5.99999994 
##    semilineatus         olssoni       insolitus       whitemani       haetianus        breslini         armouri         cybotes 
##      5.99999994      5.99999994      5.99999994      6.00000000      5.99999994      5.99999994      5.99999994      5.99999994 
##         shrevei   longitibialis         strahmi        marcanoi        baleatus       barahonae        ricordii         cuvieri 
##      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994 
##   altitudinalis        oporinus        isolepis        allisoni        porcatus        loysiana         guazuma        placidus 
##      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994      6.00000000 
##        sheplani         alayoni     angusticeps        paternus       alutaceus    inexpectatus       clivicola    cupeyalensis 
##      6.00000000      6.00000000      6.00000000      6.00000000      5.99999994      5.99999994      5.99999994      5.99999994 
##    cyanopleurus         alfaroi      macilentus       vanidicus        baracoae          noblei      smallwoodi    luteogularis 
##      5.99999994      6.00000000      6.00000000      5.99999994      6.00000000      6.00000000      6.00000000      6.00000000 
##       equestris   bahorucoensis dolichocephalus      hendersoni     darlingtoni        aliniger      singularis    chlorocyanus 
##      5.99999994      6.00000000      5.99999994      5.99999994      5.99999994      5.99999994      5.99999994      6.00000000 
##     coelestinus        occultus 
##      6.00000000      5.99999994

Now we can see that some of the tips have total height of 6.00000000 while others have a total height of 5.99999994. This is probably because one of the internal edges of the tree was rounded downward by about 0.00000006 at some point – perhaps when the phylogeny was written to or read from file.

We can fix this using force.ultrametric as follows:

new.tree<-force.ultrametric(anoletree)

## ***************************************************************
## *                          Note:                              *
## *    force.ultrametric does not include a formal method to    *
## *  ultrametricize a tree & should only be used to coerce to   *
## *   a phylogeny that fails is.ultramtric due to rounding --   *
## *   not as a substitute for formal rate-smoothing methods.    *
## ***************************************************************

Our new tree has different total height as follows:

h<-setNames(sapply(1:Ntip(new.tree),
    nodeheight,tree=new.tree),
    new.tree$tip.label)
h

##            ahli         allogus     rubribarbus           imias          sagrei         bremeri quadriocellifer      ophiolepis 
##     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959 
##         mestrei           jubar      homolechis        confusus           guafe         garmani        opalinus         grahami 
##     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959 
##     valencienni      lineatopus       evermanni       stratulus           krugi      pulchellus       gundlachi       poncensis 
##     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959 
##           cooki    cristatellus    brevirostris        caudalis          marron        websteri       distichus         alumina 
##     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959 
##    semilineatus         olssoni       insolitus       whitemani       haetianus        breslini         armouri         cybotes 
##     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959 
##         shrevei   longitibialis         strahmi        marcanoi        baleatus       barahonae        ricordii         cuvieri 
##     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959 
##   altitudinalis        oporinus        isolepis        allisoni        porcatus        loysiana         guazuma        placidus 
##     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959 
##        sheplani         alayoni     angusticeps        paternus       alutaceus    inexpectatus       clivicola    cupeyalensis 
##     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959 
##    cyanopleurus         alfaroi      macilentus       vanidicus        baracoae          noblei      smallwoodi    luteogularis 
##     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959 
##       equestris   bahorucoensis dolichocephalus      hendersoni     darlingtoni        aliniger      singularis    chlorocyanus 
##     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959     5.999999959 
##     coelestinus        occultus 
##     5.999999959     5.999999959

It also now passes is.ultrametric:

is.ultrametric(new.tree)

## [1] TRUE

In this case, we have a new total height of the tree that is the (weighted) average of our two different previous heights. This is because the default method of force.ultrametric is to use non-negative least squares to minimize the sum of squared difference between the edge lengths of our tree and the edge lengths of an ultrametric tree with the same topology.

The implied patristic distances on these two trees should also be (virtually) identical:

D1<-cophenetic(anoletree)
D2<-cophenetic(new.tree)[rownames(D1),
    colnames(D1)]
plot(D1,D2,bty="n",pch=21,cex=1.2,bg="grey",
    xlab="patristic distances on original tree",
    ylab="patristic distances on NNLS tree")
lines(c(0,12),c(0,12),lty="dotted")
legend("topleft","1:1 line",lty="dotted",bty="n")

If our tree comes from a method (like ML estimation using RAxML, for instance) that does not result in an ultrametric tree, my recommendation would be to use a formal method such as penalized likelihood (Sanderson 2002). In R, this is implemented in the ape function chronos.