Estimation of missing tip states under Brownian motion using likelihood

As I have noted previously, it is possible (and relatively straightforward) to reconstruct missing tip values under Brownian motion by jointly maximizing the likelihood of internal and (missing) terminal states using likelihood.

This is implemented in the phytools function anc.ML.

One attribute of these reconstructed values is that they will precisely match the reconstructed ancestral node state of their immediate parent. This is totally unsurprising because under Brownian motion because the expected change over any time period is zero (with variance equal to the product of the Brownian rate, σ², and the edge length).

We can visualize this phenomenon using both phenogram and contMap. Here's what I mean.

First, simulate a tree & some data:

library(phytools)
tree<-pbtree(n=26,tip.label=LETTERS,scale=1)
x<-fastBM(tree)

Subsample the data so there is missing information for some tips:

xp<-sample(x,20)

Fit our ancestral states & missing terminals using ML:

fit<-anc.ML(tree,xp)
fit

## Ancestral character estimates using anc.ML under a BM model:
##        27        28        29        30        31        32        33 
## -0.851047 -0.824691 -0.884692 -1.184654 -1.128078 -1.402134 -1.580146 
##        34        35        36        37        38        39        40 
## -1.642514 -1.551363 -1.477286 -1.374941 -0.878525 -1.194390  0.306645 
##        41        42        43        44        45        46        47 
##  0.450288  0.539937  0.568322  0.527984 -0.806273 -0.753233 -0.514332 
##        48        49        50        51 
## -0.251404 -1.112156 -0.722011 -1.412652 
## 
## Fitted model parameters & likelihood:
##      sig2 log-likelihood
##  0.468265       5.497266
## 
## R thinks it has found the ML solution.

Note that reconstructed missing tip values match their parents (to fairly high numerical precision - this is obtained via numerical optimization, remember):

tips<-sapply(names(fit$missing.x),function(x,y)
    which(y==x),y=tree$tip.label)
library(phangorn)
parents<-sapply(tips,Ancestors,x=tree,type="parent")
X<-rbind(fit$missing.x,fit$ace[as.character(parents)])
rownames(X)<-c("daughter","parent")
X

##                  H         I         K         L         O          W
## daughter -1.580114 -1.402163 -1.194513 0.3066923 0.5279829 -0.7220155
## parent   -1.580146 -1.402134 -1.194390 0.3066445 0.5279840 -0.7220109

Now let's visualize it - first with phenogram:

x.all<-c(xp,fit$missing.x,fit$ace)
tree<-paintBranches(tree,tips,"2")
phenogram(tree,x.all)
nodelabels(node=tips,pie=rep(1,length(fit$missing.x)),
    cex=0.5)

or with contMap:

obj<-contMap(tree,x.all,method="user",anc.states=fit$ace,
    plot=FALSE)
obj<-setMap(obj,colors=c("white","black"))
plot(obj)
nulo<-sapply(names(fit$missing.x),add.arrow,tree=obj,
    arrl=0.04,hedl=0.02,col="red")

Finally, let's repeat the whole thing just to show it wasn't a fluke:

tree<-pbtree(n=26,tip.label=LETTERS,scale=1)
x<-fastBM(tree)
xp<-sample(x,20)
fit<-anc.ML(tree,xp)
tips<-sapply(names(fit$missing.x),function(x,y)
    which(y==x),y=tree$tip.label)
parents<-sapply(tips,Ancestors,x=tree,type="parent")
X<-rbind(fit$missing.x,fit$ace[as.character(parents)])
rownames(X)<-c("daughter","parent")
X

##                 B         L         N        R         V         Y
## daughter 1.430460 0.9306579 0.9480089 1.449290 0.2423973 0.4115517
## parent   1.430375 0.9306499 0.9480131 1.449279 0.2425937 0.4116535

x.all<-c(xp,fit$missing.x,fit$ace)
tree<-paintBranches(tree,tips,"2")
phenogram(tree,x.all)
nodelabels(node=tips,pie=rep(1,length(fit$missing.x)),
    cex=0.5)

obj<-contMap(tree,x.all,method="user",anc.states=fit$ace,
    plot=FALSE)
obj<-setMap(obj,colors=c("white","black"))
plot(obj)
nulo<-sapply(names(fit$missing.x),add.arrow,tree=obj,
    arrl=0.04,hedl=0.02,col="red")

Recently, it was reported to me that for some empirical dataset this has not been observed to be the case. The most likely explanation for this is non- convergence of the ML optimization by anc.ML. This can be ameliorated by increasing maxit, but this may not guarantee convergence on relatively large trees. I'm sure that a better solution exists, unfortunately I have not worked on this in a bit.

Animated optimization method for cophylo

A couple of weeks ago I pushed a few updates to the phytools function for computing a co-phylogenetic object of class "cophylo" that animates the tree-traversal & optimization process.

Note that cophylo uses a greedy optimization routine involving pre-order tree-traversals and node rotations to attempt to minimize the objective function - by default the sum of squared vertical deviations of the tip heights in the two trees.

Here is a looping .gif demo of this animation:

library(phytools)
obj<-cophylo(t1,t2,anim.cophylo=TRUE,link.type="curved",
    link.lwd=3,link.lty="solid",lwd=2,
    link.col=make.transparent("blue",0.25),fsize=1.4)

The final result is here:

plot(obj,link.type="curved",link.lwd=3,link.lty="solid",
    lwd=2,link.col=make.transparent("blue",0.25),fsize=1.2)

The updates permitting this visualization are here and can be obtained by installing phytools from CRAN.

Finally, if you too want to create an animated .gif using R, the easiest way to do it is by installing ImageMagick. Then you can do something like the following:

png(file="cophylo-%03d.png",width=600,height=600)
obj<-cophylo(t1,t2,anim.cophylo=TRUE,link.type="curved",
    link.lwd=3,link.lty="solid",lwd=2,
    link.col=make.transparent("blue",0.25),fsize=1.4)
dev.off()
system("magick convert -delay 20 -loop 0 *.png cophylo-anim.gif")
file.remove(list.files(pattern=".png"))

Evaluating all possible co-phylogenetic plots for a pair of trees

I just pushed an update to the co-phylogenetic phytools method cophylo to permit it to compute & return all possible rotated objects of class "cophylo" for a pair of trees and table of associations.

This is functionally equivalent to using allRotations to rotate both trees in all possible ways, and then cophylo to compute an object of class "cophylo" (without rotation) on each pair of trees.

This will generally result in a large number of "cophylo" objects for even a modest number of taxa in the two input trees, as follows:

X

##    N1=N2=   N(cophylo)
## 1       2            4
## 2       3           16
## 3       4           64
## 4       5          256
## 5       6         1024
## 6       7         4096
## 7       8        16384
## 8       9        65536
## 9      10       262144
## 10     11      1048576
## 11     12      4194304
## 12     13     16777216
## 13     14     67108864
## 14     15    268435456
## 15     16   1073741824
## 16     17   4294967296
## 17     18  17179869184
## 18     19  68719476736
## 19     20 274877906944

Let's try the update:

library(phytools)
t1<-rtree(n=5,tip.label=LETTERS[1:5])
t2<-rtree(n=5,tip.label=letters[1:5])
assoc<-cbind(LETTERS[1:5],letters[1:5])
obj<-cophylo(t1,t2,methods="all",assoc=assoc)

## Rotating nodes to optimize matching...
## Done.

match.score<-sapply(obj,function(x) sum(attr(x$trees[[1]],"minRotate"),
    attr(x$trees[[2]],"minRotate")))
best<-which(match.score==min(match.score))

par(mfrow=c(16,16))
par(mar=rep(0,4))

for(i in 1:length(obj)){
    if(i%in%best) par(fg="blue") else par(fg="black")
    plot(obj[[i]],
    fsize=0.4,link.lty="solid",
    link.col=if(i%in%best) "blue" else "black",
    pts=FALSE,ylim=c(-0.1,1.1))
}

Here I have plotted the best matching rotations in blue.

That's it.