Saturday, November 28, 2015

Some bug fixes in phytools phenogram/fastAnc, rerootingMethod, & fitMk

I just spent some time this afternoon fixing a few different relatively small bugs in phytools, as follows.

(1) First, I fixed an issue identified by Emmanuel Paradis with the pre-release version of ape curently available from Emmanuel's website.

He noticed that the example:

library(phytools)
example(phenogram)

failed with the following error message (in French):

phngrmR> tree<-pbtree(n=20,scale=2) 
phngrmR> tree<-sim.history(tree,Q=matrix(c(-1,1,1,-1),2,2),anc="1")
Done simulation(s).
phngrmR> # simulate in which the rate depends on the state
phngrmR> x<-sim.rates(tree,c(1,10))
names absent from sig2: assuming same order as $mapped.edge
phngrmR> phenogram(tree,x)
Error in tree$mapped.edge[ii, ] : indice hors limites

This was because phenogram (for plotting a type of projection of the phylogeny into morphospace sometimes called a 'traitgram') calls fastAnc internally which was not working for "simmap" class objects. The reason for this is because Emmanuel has changed the code of root which fastAnc uses internally so that it no longer works properly for "simmap" class trees. The solution was pretty straightforward - and was pointed out by Emmanuel - and that is to just strip the class attribute "simmap" from the tree internally in fastAnc. This is somewhat unsatistfying, but it will have to do for now. This fix can be seen here.

(2) Next, it was reported to me that the following simply example:

tree <- read.tree(text = "((C:58.11154891,B:62.90277913):10.52003465,A:54.16296399);")
tipvals <- matrix(c(1, 0, 0, 0, 1, 1), nrow = 3, dimnames = list(c("A", "B", "C"), c("0", "1")))
rerootingMethod(tree, tipvals, model = "ER")

also generated a peculiar error as follows:

Error in fitMk(tt, yy, model = model, fixedQ = Q, output.liks = TRUE)$lik.anc[1,  : 
  incorrect number of dimensions

This turns out to be due to a peculiarity of R that when accessing a single row of a matrix you get a vector. If you try to then tree that vector as a matrix, R throws an error. This bug is thus due to sending rerootingMethod a tree with only three tips. When that tree is rerooted at its one & only internal node (besides the root) the matrix of conditional likelihoods returned by fitMk (used internally) contains only one row and is thus a vector.

My fix was merely to check if the tree has three tips, and then make sure that the object returned is a matrix not a vector. Now the function seems to work fine for the case triggerint the error, as follows:

library(phytools)
tree <- read.tree(text = "((C:58.11154891,B:62.90277913):10.52003465,A:54.16296399);")
tipvals <- matrix(c(1, 0, 0, 0, 1, 1), nrow = 3, dimnames = list(c("A", "B", "C"), c("0", "1")))
rerootingMethod(tree, tipvals, model = "ER")
## $loglik
## [1] -2.079442
## 
## $Q
##      0    1
## 0 -0.1  0.1
## 1  0.1 -0.1
## 
## $marginal.anc
##           0         1
## C 0.0000000 1.0000000
## B 0.0000000 1.0000000
## A 1.0000000 0.0000000
## 4 0.5000091 0.4999909
## 5 0.4999939 0.5000061

(3) Finally, in the processing of addressing (2), I also discovered a weird bug in fitMk in which the diagonal elements of the transition matrix Q are two times too large for model="ER". The reason for this (though simple) is kind of hard to explain. Suffice it to say that it should now be fixed. The fix can be seen here.

That's all I have to post about this, I think.

2 comments:

  1. Did you consider using [,,drop=FALSE] for (2)? I end up using that a lot to avoid the 'nasty surprises' resulting from the one-row-subsetting issue.

    ReplyDelete
    Replies
    1. I forgot about that. Yes, that is a better solution. Changed!

      Delete