Update to rerootingMethod for ancestral state reconstruction to permit polytomies

A phytools user recently reported that the ancestral character estimation method for discrete traits in phytools, rerootingMethod, does not work for trees with polytomies. This is in fact correct - although I had long forgotten that this was the case.

The function rerootingMethod uses the approach of Yang et al. (1995) to compute the marginal ancestral states by re-rooting the tree at all internal nodes. For instances in which some states at the tips are uncertain, it can also re-root at the tip to compute empirical Bayesian posterior probabilities for the leaves as well. The function is less relevant now that ace in the ape package also can compute the marginal ancestral states. (Although, as far as I know, ace can also not handle polytomies.) rerootingMethod stays somewhat relevant only by virtue of being able to compute posterior probabilities at tips, and for handling uncertainty at tip states (two sides of the same coin, I guess).

It was fairly straightforward to update rerootingMethod to handle polytomies, and the code is posted here. The way that this current version works is pretty simple. It just first takes the input tree, checks for polytomies using is.binary.tree in the ape package, resolves polytomies randomly with branches of zero length using multi2di, estimates ancestral states for all internal nodes (& tips, if tips=TRUE), and the uses the phytools function matchNodes to back-translate the reconstructed nodes in the resolved tree to the original, input tree. Note that although it is doing something here that ace does not do - the code internally uses (modified code from) ace.

Here - we can check it out by simulating a tree with some polytomies, generating trait data for tips & nodes, and then reconstructing ancestral states on the tree:

library(phytools)
dt<-rtree(n=26)
dt$tip.label<-LETTERS
## set some of the internal branches
ii<-which(dt$edge[,2]>Ntip(dt))
dt$edge.length[sample(ii,4)]<-0
mt<-di2multi(dt)
plotTree(mt)

Next let's simulate a character up the tree:

Q<-matrix(c(-1,0.5,0.5,
    0.5,-1,0.5,0.5,0.5,-1),3,3,
    dimnames=list(letters[1:3],letters[1:3]))
st<-sim.history(mt,Q)

## Done simulation(s).

x<-st$states
y<-getStates(st,"nodes")

Now, we can reconstruct & plot ancestral states using ace on the bifurcating tree; and rerootingMethod on the bi- and on the multifurcating trees. If the methods are working properly then the reconstructions will be the same - although in the bifurcating node there will be additional nodes with states identical to other nodes from which they are separated by branches of zero length.

## first ace
fit1<-ace(x,mt,type="discrete",model="ER") ## doesn't work

## Error in ace(x, mt, type = "discrete", model = "ER"): "phy" is not rooted AND fully dichotomous.

fit2<-ace(x,dt,type="discrete",model="ER") ## doesn't work

## Error in ace(x, dt, type = "discrete", model = "ER"): some branches have length zero or negative

## set zero-length branches to be 1/1000000 total tree length
dst<-dt
dst$edge.length[dst$edge.length==0]<-max(nodeHeights(dt))*1e-6
fit3<-ace(x,dst,type="discrete",model="ER")
fit3

## 
##     Ancestral Character Estimation
## 
## Call: ace(x = x, phy = dst, type = "discrete", model = "ER")
## 
##     Log-likelihood: -25.3342 
## 
## Rate index matrix:
##   a b c
## a . 1 1
## b 1 . 1
## c 1 1 .
## 
## Parameter estimates:
##  rate index estimate std-err
##           1   0.3278  0.1042
## 
## Scaled likelihoods at the root (type '...$lik.anc' to get them for all nodes):
##          a          b          c 
## 0.82837257 0.09781827 0.07380916

Now with the new version of rerootingMethod:

packageVersion("phytools")

## [1] '0.4.57'

fit4<-rerootingMethod(mt,x,model="ER")
fit4

## $loglik
## [1] -25.33421
## 
## $Q
##            a          b          c
## a -0.6556595  0.3278298  0.3278298
## b  0.3278298 -0.6556595  0.3278298
## c  0.3278298  0.3278298 -0.6556595
## 
## $marginal.anc
##               a          b            c
## 27 0.8283833959 0.09781145 0.0738051507
## 28 0.2121485528 0.69489845 0.0929529991
## 29 0.0532939602 0.91028552 0.0364205221
## 30 0.0627113681 0.88952034 0.0477682947
## 31 0.0132238730 0.97569370 0.0110824277
## 32 0.3870898035 0.30541903 0.3074911711
## 33 0.4035347684 0.15548939 0.4409758369
## 34 0.5865138808 0.07266794 0.3408181834
## 35 0.3530969369 0.06523786 0.5816652004
## 36 0.7045263050 0.13922658 0.1562471186
## 37 0.1288634178 0.06087131 0.8102652730
## 38 0.3401006710 0.49952777 0.1603715570
## 39 0.4129576403 0.32707964 0.2599627226
## 40 0.1266883464 0.35550347 0.5178081798
## 41 0.0553700550 0.21371716 0.7309127832
## 42 0.8088376946 0.10680265 0.0843596595
## 43 0.8570510600 0.08945954 0.0534894017
## 44 0.6274100032 0.30199393 0.0705960682
## 45 0.0836944308 0.86315981 0.0531457619
## 46 0.0119373448 0.97893201 0.0091306462
## 47 0.0006664523 0.99874487 0.0005886784

If we plot both of these on our trees in turn, we should see that the marginal reconstructions (excepting the small deviation that was required to give our tree non-zero branch lengths throughout) should be equal:

## first ace
plotTree(dt,offset=0.5)
nodelabels(pie=fit3$lik.anc)
tiplabels(pie=to.matrix(x,seq=letters[1:3]),cex=0.5)

## now phytools
plotTree(mt,offset=0.5)
nodelabels(pie=fit4$marginal.anc)
tiplabels(pie=to.matrix(x,seq=letters[1:3]),cex=0.5)

Finally, for fun let's overlay the true, known states to see how close (or off) we were:

plotSimmap(st,colors=setNames(c("red","green","blue"),letters[1:3]),
    offset=0.5)
nodelabels(pie=fit4$marginal.anc)
tiplabels(pie=to.matrix(x,seq=letters[1:3]),cex=0.5)

That's it.

Calculating the times spent with each of number of lineages on the tree

Yesterday a phytools user emailed me the following inquiry:

“Is there an easy way in phytools to infer the times while there are exactly k ancestral lineages in the tree? I would need to calculate these based on a newick tree input. However, I havent found a function in phytools for it. Somthing like k=2 from x to y k=3 from z to v”

This can be done pretty easily using the phytools function ltt. Of course, one might allow for the possilibity that - if lineages both increase & decrease through time, there could be more than one period of time during which there is 2 lineages, 3 lineages etc.

The following gives a quick demo - first using a tree in which the number of lineages both increases & descreases in time since the root:

library(phytools)
tree<-rtree(26,tip.label=LETTERS)
plotTree(tree)

## create an object with the number of lineages through time:
obj<-ltt(tree,log=FALSE)

Now let's go through all the number of possible lineages and tabulate all the time intervals in which the tree has each number of lineages:

foo<-function(n,obj){
    ii<-which(obj$ltt==n)
    ii<-ii[ii!=length(obj$ltt)]
    sapply(ii,function(i,x) c(x[i],x[i+1]),x=obj$times)
}
Nlineages<-sapply(1:max(obj$ltt),foo,obj=obj)
if(is.matrix(Nlineages)){
    colnames(Nlineages)<-1:max(obj$ltt)
    rownames(Nlineages)<-c("start","end")
} else if(is.list(Nlineages)){
    names(Nlineages)<-1:max(obj$ltt)
    for(i in 1:length(Nlineages)) 
        rownames(Nlineages[[i]])<-c("start","end")
}

Our object is a list, with a matrix for each number of lineages showing the time intervals spent with that number of lineages.

Nlineages

## $`1`
##               16
## start 0 4.128545
## end   0 4.198241
## 
## $`2`
##               27       15
## start 0.00000000 4.049935
## end   0.07445996 4.128545
## 
## $`3`
##               28       13       11
## start 0.07445996 1.501722 3.962453
## end   0.15032225 1.711749 4.049935
## 
## $`4`
##              30         1         4        2       41       22
## start 0.1503222 0.3991084 0.5425447 1.499231 1.711749 3.655242
## end   0.2915629 0.4658409 0.8433361 1.501722 2.151097 3.962453
## 
## $`5`
##              31        50        29       24       35       21
## start 0.2915629 0.4658409 0.8433361 1.464400 2.151097 3.651133
## end   0.3991084 0.5425447 0.8714423 1.499231 2.194213 3.655242
## 
## $`6`
##              51        3       37        5        9       18
## start 0.8714423 1.458014 2.194213 2.331990 3.612215 3.642776
## end   0.9738692 1.464400 2.274517 2.353893 3.618789 3.651133
## 
## $`7`
##              32       26       25       45       42       14       12
## start 0.9738692 1.096267 1.325357 2.274517 2.353893 2.607869 3.509315
## end   1.0088520 1.263415 1.458014 2.331990 2.573532 2.785814 3.612215
##             44
## start 3.618789
## end   3.642776
## 
## $`8`
##             33       34       46       38        7
## start 1.008852 1.263415 2.573532 2.785814 3.484318
## end   1.096267 1.325357 2.607869 2.810022 3.509315
## 
## $`9`
##             47       20        6
## start 2.810022 3.318845 3.382049
## end   2.953015 3.380557 3.484318
## 
## $`10`
##             39       23       19        8       40
## start 2.953015 3.103372 3.199463 3.277445 3.380557
## end   2.988247 3.114297 3.252606 3.318845 3.382049
## 
## $`11`
##             43       36       17       49
## start 2.988247 3.114297 3.179681 3.252606
## end   3.103372 3.149976 3.199463 3.277445
## 
## $`12`
##             48
## start 3.149976
## end   3.179681

Next, we can try and ultrametric trees with all lineages terminating in the present:

tree<-pbtree(n=26,tip.label=LETTERS)
plotTree(tree)

obj<-ltt(tree,log=FALSE)

foo<-function(n,obj){
    ii<-which(obj$ltt==n)
    ii<-ii[ii!=length(obj$ltt)]
    sapply(ii,function(i,x) c(x[i],x[i+1]),x=obj$times)
}
Nlineages<-sapply(1:max(obj$ltt),foo,obj=obj)
if(is.matrix(Nlineages)){
    colnames(Nlineages)<-1:max(obj$ltt)
    rownames(Nlineages)<-c("start","end")
} else if(is.list(Nlineages)){
    names(Nlineages)<-1:max(obj$ltt)
    for(i in 1:length(Nlineages)) 
        rownames(Nlineages[[i]])<-c("start","end")
}

This time we have a matrix in which each column is a number of lineages and the rows show the start and end times. When no lineages are lost through time, there will be only one time interval for each number of lineages:

Nlineages

##           1            2         3         4        5        6        7
## start 0e+00 2.000000e-11 0.4418836 0.5968292 1.356948 1.752875 1.756543
## end   2e-11 4.418836e-01 0.5968292 1.3569481 1.752875 1.756543 1.876857
##              8        9       10       11       12       13       14
## start 1.876857 1.931771 2.000830 2.022528 2.046486 2.218467 2.283606
## end   1.931771 2.000830 2.022528 2.046486 2.218467 2.283606 2.445454
##             15       16      17       18       19       20       21
## start 2.445454 2.506024 2.57934 2.586070 2.596797 2.668859 2.671522
## end   2.506024 2.579340 2.58607 2.596797 2.668859 2.671522 2.735799
##             22       23       24       25       26
## start 2.735799 2.773744 2.779681 2.823757 2.865986
## end   2.773744 2.779681 2.823757 2.865986 2.866205

Note that the tiny zeroeth timeslice with one lineage at the start of each tree is an idiosyncratic feature of phytools::ltt and could easily be ignored.

That's it.

T-shirt design for Ilhabela workshop - in R

I'm in the midst of working on a design for the t-shirts for this year's macroevolution workshop in Ilhabela, Brazil. I've created a tentative 'front' design entirely in R using the function of phytools, as follows:

library(phytools)
## don't run
# tree<-pbtree(n=100)
# x<-fastBM(tree)

obj<-contMap(tree,x,plot=FALSE)
obj<-setMap(obj,colors=c("blue","purple","red"))
layout(mat=matrix(c(1,2),2,1),heights=c(0.8,0.2))
par(bg="black")
par(fg="white")
plotSimmap(paintSubTree(tree,Ntip(tree)+1,"1"),type="fan",
    ftype="off",colors=setNames("white","1"),lwd=6,part=0.5)

## setEnv=TRUE for this type is experimental. please be patient with bugs

plotSimmap(obj$tree,type="fan",ftype="off",colors=obj$cols,lwd=4,
    add=TRUE,part=0.5)

## setEnv=TRUE for this type is experimental. please be patient with bugs

plot.new()
text(0.5,0.5,"Latin American Macroevolution Workshop\nIlhabela Brazil 2015",
    col="white",cex=2.1,font=2)

The knitr output is kind of aliased. We can create a higher quality PDF easily:

pdf(file="t-shirt.pdf",width=8,height=4.75)
layout(mat=matrix(c(1,2),2,1),heights=c(0.8,0.2))
par(bg="black")
par(fg="white")
plotSimmap(paintSubTree(tree,Ntip(tree)+1,"1"),type="fan",
    ftype="off",colors=setNames("white","1"),lwd=6,part=0.5)

## setEnv=TRUE for this type is experimental. please be patient with bugs

plotSimmap(obj$tree,type="fan",ftype="off",colors=obj$cols,lwd=4,
    add=TRUE,part=0.5)

## setEnv=TRUE for this type is experimental. please be patient with bugs

plot.new()
text(0.5,0.5,"Latin American Macroevolution Workshop\nIlhabela Brazil 2015",
    col="white",cex=2.1,font=2)
dev.off()

## windows 
##       2

Here is a link to the file. (It looks much better!)

That's it for now.