Updates to fitMk to permit more user control of optimization

I just pushed some relatively small updates to the phytools function fitMk, which (aptly) fits the so-called Mk discrete character evolution model. The updates allow for more user control of optimization - specifically by allowing the user to specify an optimization method (opt.method), which can be "nlminb" or "optim", and for the user to set initial values for the model parameters.

Note that fitMk takes a lot of code from the ape function ace (for ancestral character estimation); however fitMk allows the user to specify a fixed value of the transition matrix, Q, a non-flat prior probability density, pi, ambiguous states at the tips of the tree in the form of a matrix, and (now) an optimization method and starting values for the optimization. The function fitDiscrete in geiger is quite powerful for fitting this model as well.

These updates are in part inspired by a recent post I wrote about an interesting failure of optimization in the function. Since we are doing numerical optimization, not exhaustive enumeration, there is always a chance we will fail to find the ML solution! Hopefully these updates will allow users to be more confident that they have in fact converged on the tree MLEs for their model.

fitMk is also used internally by make.simmap and rerootingMethod (and perhaps other functions that I'm not thinking of now); however I have not yet migrated these options to these other methods.

Here's a demo with our dataset from before. This is a real empirical dataset, but one in which the taxon names have been obfuscated:

library(phytools)
packageVersion("phytools")

## [1] '0.5.52'

tree

## 
## Phylogenetic tree with 271 tips and 270 internal nodes.
## 
## Tip labels:
##  t1, t2, t3, t4, t5, t6, ...
## 
## Rooted; includes branch lengths.

head(x,n=100); cat("....\n")

##   t1   t2   t3   t4   t5   t6   t7   t8   t9  t10  t11  t12  t13  t14  t15 
##    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0 
##  t16  t17  t18  t19  t20  t21  t22  t23  t24  t25  t26  t27  t28  t29  t30 
##    0    1    1    0    0    0    0    0    0    0    0    0    1    0    0 
##  t31  t32  t33  t34  t35  t36  t37  t38  t39  t40  t41  t42  t43  t44  t45 
##    0    0    0    0    0    0    0    0    0    0    1    0    0    0    0 
##  t46  t47  t48  t49  t50  t51  t52  t53  t54  t55  t56  t57  t58  t59  t60 
##    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0 
##  t61  t62  t63  t64  t65  t66  t67  t68  t69  t70  t71  t72  t73  t74  t75 
##    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0 
##  t76  t77  t78  t79  t80  t81  t82  t83  t84  t85  t86  t87  t88  t89  t90 
##    0    1    0    0    0    0    0    0    0    0    1    0    1    0    0 
##  t91  t92  t93  t94  t95  t96  t97  t98  t99 t100 
##    0    0    0    0    0    0    0    0    1    0

## ....

plotTree(tree,type="fan",lwd=1,ftype="off")
tiplabels(pie=to.matrix(x[tree$tip.label],c(0,1)),
    piecol=c("blue","red"),cex=0.3)
add.simmap.legend(colors=setNames(c("blue","red"),
    c(0,1)),x=-176,y=169,prompt=FALSE,shape="circle")

OK, let's fit our model:

fit.nlminb<-fitMk(tree,x,opt.method="nlminb")
fit.nlminb

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           0         1
## 0 -1.000148  1.000148
## 1  1.000148 -1.000148
## 
## Fitted (or set) value of pi:
##   0   1 
## 0.5 0.5 
## 
## Log-likelihood: -187.842905 
## 
## Optimization method used was "nlminb"

This was a rate that was about 200 × too high! Let's try the "optim" method for optimization. As you might guess, this substitutes the optim function for nlminb internally, and is a little bit slower as a consequence:

fit.optim<-fitMk(tree,x,opt.method="optim")
fit.optim

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           0         1
## 0 -0.005093  0.005093
## 1  0.005093 -0.005093
## 
## Fitted (or set) value of pi:
##   0   1 
## 0.5 0.5 
## 
## Log-likelihood: -155.900603 
## 
## Optimization method used was "optim"

This gives us the correct answer this time; however it may not in general. For instance if we initial our transition rate with a value that is too high it will fail to converge on the correct solution:

fit.optim.q0.1<-fitMk(tree,x,opt.method="optim",q.init=0.1)
fit.optim.q0.1

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           0         1
## 0 -1.389036  1.389036
## 1  1.389036 -1.389036
## 
## Fitted (or set) value of pi:
##   0   1 
## 0.5 0.5 
## 
## Log-likelihood: -187.842886 
## 
## Optimization method used was "optim"

Finally, as I showed before it is possible to visualize the likelihood surface for the transition rate. This is also possible in two dimensions for different forward & backward rates, q_0,1 & q_1,0. For instance:

## one dimension
q<-seq(0.001,0.1,by=0.001)
logL<-sapply(q,function(q,tree,x) 
    logLik(fitMk(tree,x,fixedQ=matrix(c(-q,q,q,-q),2,2))),
    tree=tree,x=x)
plot(q,logL,type="l",xlab="q[0,1] & q[1,0]",ylab="log(L)")

## two dimensions
logL<-sapply(q,function(q1,q2,tree,x) sapply(q2,
    function(q2,q1,tree,x) logLik(fitMk(tree,x,fixedQ=
    matrix(c(-q1,q1,q2,-q2),2,2))),tree=tree,x=x,q1=q1),
    tree=tree,x=x,q2=q)
contour(logL,x=q,y=q,nlevels=100,xlab="q[0,1]",ylab="q[1,0]",
    main="log(L)")
fitARD<-fitMk(tree,x,model="ARD")
points(fitARD$rates[2],fitARD$rates[1],col="red",pch=19)

Cool.

Collapsing clades of 'foo' - huh?

A R-sig-phylo subscriber recently posted a request as follows:

“[R-sig-phylo] Collapse a clade by tip labels while maintaining
phylogenetic position branchlizard . Mon, 12 Sep 2016 12:46:56 -0700

I have posted this question at Stack Overflow. I hope this doesn't violate any community rules about double posting.

I probably could have worded the title better, but I am wanting to collapse any clade within a phylogenetic tree (even if the clade has one member) which has a tip label of "foo” and then count the number of tips which were dropped from that specific clade and create a branch with a tip label displaying 35 foos.

The counting portion is easy; however, when I use

drop.tip(rooted.tree,tip=which(rooted.tree$tip.label=='foo'),subtree=TRUE)

the dropped tips do not maintain their position in the tree. Rather, they are all grouped at the end (counted properly however). Is there anyway to collapse a clade by tip labels and maintain its position"

When asked to clarify, the poster responded that what he or she wants to do is convert a tree like this:

into a tree like this:

Florian Boucher posted a solution that probably works, I haven't tried. Here is another one:

First load packages and make a tree with this property:

library(phytools)
library(phangorn)

text<-"(((((foo:0.9,foo:0.7):0.7,((foo:0.6,foo:0.9):0.9,foo:0.8):0.5):0.9,A:0.8):0.9,(B:0.1,C:0.1):0.5):0.6,((((foo:0.7,(foo:0.4,foo:0.8):0.9):0.3,(foo:0.9,(foo:0.7,foo:0.2):0.5):0.6):0.7,D:0.7):0.4,((foo:0.3,((((foo:0.3,foo:0.5):0.1,foo:0.5):0.4,foo:0.7):0.3,E:0.1):0.6):1,((F:0.2,G:0.8):0.7,((foo:0.7,foo:0.3):0.2,foo:0.7):0.1):0.8):0.9):0.8);"
tree<-read.tree(text=text)
plotTree(tree)

This is our uncollapsed tree. Now let's go about identifying all the nodes that we want to collapse and collapsing them. One complication is that each time we collapse one clade, the node numbers will, of course, change. A further wrinkle is that because tip labels repeat, it might be difficult to match tips and nodes in the tipical way. Finally, we might want to maintain the tip a certain height about the root - for instance the mean height of the 'foo's of its clade:

nodes<-1:tree$Nnode+Ntip(tree) ## all nodes
subtrees<-list()
for(i in 1:tree$Nnode) subtrees[[i]]<-extract.clade(tree,nodes[i])
names(subtrees)<-nodes ## all subtrees
## all nodes with only "foo" as descendant
all.foos<-nodes[sapply(subtrees,function(x) all(x$tip.label=="foo"))]
foo.mrcas<-all.foos[sapply(all.foos,function(x,tree,y)
    !Ancestors(tree,x,"parent")%in%y,
    tree=tree,y=all.foos)]
w.foos<-which(tree$tip.label=="foo")
## rename tips uniquely
tree$tip.label[w.foos]<-paste(tree$tip.labe[w.foos],
    replicate(length(w.foos),paste(sample(letters,6),
    collapse="")),sep="_")
collapsed<-tree
## iterate over all MRCAs of foo clades
for(i in 1:length(foo.mrcas)){
    M<-matchNodes(tree,collapsed)
    nn<-M[which(M[,1]==foo.mrcas[i]),2]
    dd<-Descendants(collapsed,nn)[[1]]
    h<-sapply(dd,nodeheight,tree=collapsed)
    collapsed$tip.label[dd[1]]<-paste(length(dd),"foo(s)")
    ind<-which(collapsed$edge[,2]==dd[1])
    collapsed$edge.length[ind]<-collapsed$edge.length[ind]+
        mean(h)-h[1]
    if(length(dd)>1)
        collapsed<-drop.tip(collapsed,
            collapsed$tip.label[dd[2:length(dd)]])
}
## finally, address 'singleton' foos, if they exist
ind<-grep("foo_",collapsed$tip.label)
if(length(ind)>0) collapsed$tip.label[ind]<-"1 foo(s)"

Now, let's check our tree:

plotTree(collapsed)

voila!

New exhaustive search for cophylo method (and new object class)

I finally figured out how to perform an exhaustive search for tip matching in phytools co-phylogenetic plotting method, cophylo.

First, I wrote the function allRotations (described here) with the idea that if some rotation of the nodes existed such that the tips matched in their left-to-right ordering, then I should just be able to exhaustively rotate all the nodes of one tree, and amongst these, there should be a rotation to match the tips of the second tree. Seems logical, right?

Well, it turns out - not so. In fact, we have to rotate the nodes of both trees in all possible ways to guarantee that we find the best matching. For some reason, it has taken me a while to get this through my skull - so I was at times very perplexed as to why the previously mentioned approach didn't seem to work.

The update with this method can be seen here. I also added a new object class for a list of "cophylo" objects called a "multiCophylo" object, with its own S3 print and plot methods, in turn. The idea is that from an exhaustive search we will invariably find at least 2 node rotations of both trees that maximize matching (one & its mirror image). Sometimes we can find various such rotations, so I thought it would make sense to return all of the corresponding rotated "cophylo", rather than arbitrarily choosing among them.

Note that the number of possible node rotations of a bifurcating tree with N tips is 2^N-1 and the number of comparisons we thus have to make for two trees containg N₁ and N₂ tips, respectively is 2^N₁-1 × 2^N₂-1. These are all quite large quantities for even modestly sized phylogenies!

Here's a demo using relatively small phylogenetic trees.

First, load packages & simulate trees that do not share phylogenetic structure, but for which a tip matching exists:

library(phytools)
packageVersion("phytools")

## [1] '0.5.51'

set.seed(99)
N<-10
t1<-rtree(n=N,tip.label=LETTERS[1:N])
t2<-rtree(n=N)
assoc<-cbind(t1$tip.label,t2$tip.label)
assoc

##       [,1] [,2] 
##  [1,] "J"  "t7" 
##  [2,] "G"  "t9" 
##  [3,] "C"  "t10"
##  [4,] "H"  "t4" 
##  [5,] "B"  "t5" 
##  [6,] "I"  "t6" 
##  [7,] "E"  "t2" 
##  [8,] "A"  "t3" 
##  [9,] "F"  "t1" 
## [10,] "D"  "t8"

plot(cophylo(t1,t2,assoc=assoc))

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

Now, let's rotate the crap out of both trees and show that with our optimization algorithm, we do not find the matching:

for(i in 1:100) t1<-rotateNodes(t1,node=sample(1:t1$Nnode+Ntip(t1),1))
for(i in 1:100) t2<-rotateNodes(t2,node=sample(1:t2$Nnode+Ntip(t2),1))
## first unrotated
plot(cophylo(t1,t2,assoc=assoc,rotate=FALSE))

## then rotated with our optimization routine
obj<-cophylo(t1,t2,assoc=assoc)

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

plot(obj)

Finally, let's try the new exhaustive search method:

obj<-cophylo(t1,t2,assoc=assoc,methods="exhaustive")

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

obj

## Object of class "multiCophylo" containg 8 objects of class "cophylo".

par(mfrow=c(4,2))
plot(obj,mar=rep(1.1,4))

That's it. Note that this is not expected to work for trees that are not perfectly bifurcating (although our optimization algorithm does, with rotate.multi=TRUE).