Plotting densityMap using grayscale

Travis Ingram commented that it would be nice to be able to plot density maps using the function densityMap in grayscale. He gave a line of code that could be modified internally to do this; however (as of today!) this is not necessary as we can instead modify our object of class "densityMap" (now returned invisibly by the function) and replot using plot.densityMap. This is what that would look like:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.60’
> # this is just to simulate some data
> tree<-pbtree(n=80,scale=1)
> Q<-matrix(c(-1,1,1,-1),2,2)
> rownames(Q)<-colnames(Q)<-c(0,1)
> tree<-sim.history(tree,Q)
> x<-tree$states
> mtrees<-make.simmap(tree,x,nsim=100)
make.simmap is sampling character histories conditioned on the transition matrix
Q =
        0        1
0 -1.25418  1.25418
1  1.25418 -1.25418
(estimated using likelihood);
and (mean) root node prior probabilities
pi =
  0  1
0.5 0.5
Done.
> # we don't care about this plot
> maps<-densityMap(mtrees,res=500)
sorry - this might take a while; please be patient
> # now let's change our colormap using Travis's code
> maps$cols[]<-grey(seq(1,0,length.out=length(maps$cols)))
> plot(maps,fsize=c(0.6,1),outline=TRUE,lwd=5)

(Click for higher res version.) Cool.

Thanks to Travis for the great suggestion & code.

New version of densityMap

The full title of this post should read New much faster version of densityMap that returns the plotted map invisibly combined with generic plotting method for special mapping object class, or something like that - but that seemed like a tongue twister. densityMap is a function to visualize the posterior sample of maps from a stochastic mapping analysis and is described in an article that I have in press at Methods in Ecology & Evolution.

I just posted a new version of densityMap (code here) that simultaneously addresses a couple of significant issues with this (otherwise pretty cool, in my opinion) function.

Firstly, prior versions of the function are way too slow. The function is still slow - just no longer way too slow. This speed-up was achieved entirely using the trick of removing the class attribute of our "multiPhylo" object which - for some reason that is not entirely clear to me - dramatically speeds up handling of the object both by apply family functions, and even in for loops.

Secondly, what makes densityMap especially annoying to work with - given that it is slow - is that to alter any of the plotting options, you need to re-compute the aggregate mappings. In the new version of densityMap (and new minor phytools build, phytools 0.2-60), densityMap plots the mapped tree as before, but also returns a special object of class "densityMap" invisibly. This can then be plotted using a call of the generic plot (to which I have added the phytools method plot.densityMap).

Here's a demo:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.60’
> # simulate tree & data
> tree<-pbtree(n=70,scale=1)
> Q<-matrix(c(-1,1,1,-1),2,2)
> rownames(Q)<-colnames(Q)<-c(0,1)
> tree<-sim.history(tree,Q)
> x<-tree$states
> # generate stochastic maps
> mtrees<-make.simmap(tree,x,nsim=100)
make.simmap is sampling character histories conditioned on the transition matrix
Q =
          0          1
0 -0.8857354  0.8857354
1  0.8857354 -0.8857354
(estimated using likelihood);
and (mean) root node prior probabilities
pi =
  0  1
0.5 0.5
Done.
> # generate density map
> # this would have been much slower before
> system.time(map<-densityMap(mtrees))
sorry - this might take a while; please be patient
  user  system elapsed
  7.85    0.12    7.97

Now, in earlier versions of densityMap if we wanted to adjust the way this plot looked (say - so that the labels don't overlap!) - we would have to recompute the aggregate mapping (which, remember, took us 8s here - but much longer in earlier versions or for bigger trees). However, here we have created the object maps, which we can then pass to plot.densityMap with our revised plotting options:

> plot(map,lwd=5,fsize=c(0.7,1),legend=0.4)

for example.

Cool.

Bug fix & update to phylomorphospace

Yesterday a phytools user identified a bug in the node-coloring of phylomorphospace. I thought that this was likely introduced when I recently did a major rewrite of the function (described here) and this seems to be correct. I have fixed this, and also added a new feature to the function that (for trees with a mapped discrete character) will automatically color nodes using the color of the mapped discrete character. The code of the updated function is here; and I have also posted a new phytools build (phytools 0.2-59).

Here's a demo of the fixed node coloring:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.59’
> # first let's simulate a tree & data
> tree<-pbtree(n=30,scale=1)
> XX<-fastBM(tree,nsim=2)
> plotTree(tree,node.numbers=T)

> # now let's say we want to plot nodes
> # descended from "49" red:
> cols<-rep("black",length(tree$tip.label)+tree$Nnode)
> names(cols)<-1:length(cols)
> cols[getDescendants(tree,49)]<-"red"
> # and everything from "40" blue:
> cols[getDescendants(tree,40)]<-"blue"
> # finally, these can even be nested
> cols[getDescendants(tree,44)]<-"yellow"
> # and plot
> phylomorphospace(tree,XX,xlab="X1",ylab="X2", control=list(col.node=cols))

Cool. (This is basically the demo that I gave in an earlier post.)

Now, when I was doing this I realized that it might be cool to be able to color the nodes - as well as the edges - according to a mapped discrete character. To do this, we need to be able to compute the colors of all internal nodes on the tree from their states. Here is the code that I used to do that:

zz<-c(getStates(tree,type ="tips"), getStates(tree))
names(zz)[1:length(tree$tip.label)]<-
  sapply(names(zz)[1:length(tree$tip.label)],
  function(x,y) which(y==x),y=tree$tip.label)
con$col.node<-setNames(colors[zz],names(zz))

The first line just computes the states at all tip & internal nodes using getStates; the second is just a complicated way of translating tip labels in names(zz) into node numbers, which is what is need by phylomorphospace; finally, line three translates the node states to colors.

Here is a demo:

> # transition matrix
> Q<-matrix(c(-2,2,2,-2),2,2)
> colnames(Q)<-rownames(Q)<-letters[1:2]
> # simulate stochastic history
> tree<-sim.history(tree,Q)
> phylomorphospace(tree,XX,xlab="X1",ylab="X2", colors=setNames(c("blue","red"),letters[1:2]), node.by.map=TRUE)

That's it.

Much faster versions of countSimmap, describe.simmap, and getStates

Yesterday, Klaus Schliep pointed out that a trick to speed up lapply's handling of objects of class "multiPhylo" (i.e., just a simple list of objects of class "phylo") was just to first remove the class attribute "multiPhylo". Why this would have any effect at all is somewhat of a mystery to me. is.list(trees) evaluates TRUE regardless of whether or not the class attribute has been removed, so it doesn't seem that lapply would have to coerce our object to a list in either case. Nonetheless - not only does this work, it works tremendously! I have now included this simple trick in countSimmap, describe.simmap, and getStates, all of which use lapply or sapply if the argument tree is an object of class "multiPhylo".

Here's a demo of just how much of an improvement in speed results from this trick:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.56’
>
> # simulate data
> tree<-pbtree(n=100,scale=1)
> Q<-matrix(c(-1,1,1,-1),2,2)
> colnames(Q)<-rownames(Q)<-letters[1:2]
> tree<-sim.history(tree,Q)
> x<-tree$states
>
> # stochastic mapping
> mtrees<-make.simmap(tree,x,nsim=1000)
make.simmap is sampling character histories conditioned on the transition matrix
Q =
          a          b
a -0.8686634  0.8686634
b  0.8686634 -0.8686634
(estimated using likelihood);
and (mean) root node prior probabilities
pi =
  a  b
0.5 0.5
Done.
>
> # ok, now let's time describe.simmap for
> # various subsets of our mapped trees
> system.time(X100<-describe.simmap(mtrees[1:100], message=FALSE))
  user  system elapsed
  2.69    0.02    2.70
> system.time(X200<-describe.simmap(mtrees[1:200], message=FALSE))
  user  system elapsed
  12.71    0.02  12.75
> system.time(X400<-describe.simmap(mtrees[1:400], message=FALSE))
  user  system elapsed
  74.24    0.64  75.57

Woah. I'm not even going to try the full set of 1,000 trees.

OK, now let's compare to describe.simmap with nothing more than the trick suggested by Klaus (i.e., unclass-ing the "multiPhylo" object for every use of lapply or sapply):

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.58’
> system.time(Y100<-describe.simmap(mtrees[1:100], message=FALSE))
  user  system elapsed
  0.45    0.00    0.45
> system.time(Y200<-describe.simmap(mtrees[1:200], message=FALSE))
  user  system elapsed
  0.83    0.00    0.82
> system.time(Y400<-describe.simmap(mtrees[1:400], message=FALSE))
  user  system elapsed
  1.78    0.00    1.78

Holy cow! What a huge improvement. We can even now run it on the full set of 1,000 mapped trees:

> par(cex=0.8) # make our tip labels a little smaller
> system.time(YY<-describe.simmap(mtrees,plot=TRUE))
1000 trees with a mapped discrete character with states:
a, b

trees have 29.633 changes between states on average

changes are of the following types:
        a,b    b,a
x->y 10.692 18.941

mean total time spent in each state is:
              a          b    total
raw  12.9314025 18.3673136 31.29872
prop  0.4131608  0.5868392  1.00000

  user  system elapsed
  4.65    0.06    4.71

Wow.

Faster version of getStates, but why is it faster...?

I was tinkering with the function describe.simmap to try and speed it up, when I discovered that the very simple function getStates, which does nothing more than get the states on a mapped tree in memory for the nodes or tips, had run-time using lapply to iterate over the trees in an object of class "multiPhylo", that seems to rise non-linearly with the number of trees.

> ## simulate tree & data
> tree<-pbtree(n=200,scale=1)
> Q<-matrix(c(-1,1,1,-1),2,2)
> colnames(Q)<-rownames(Q)<-letters[1:2]
> trees<-sim.history(tree,Q,nsim=200)
> # now lets get the states for all trees at all nodes
> # using sapply
> system.time(XX1<-getStates(trees[[1]]))
  user  system elapsed
  0.02    0.00    0.02
> system.time(XX10<-sapply(trees[1:10],getStates))
  user  system elapsed
  0.05    0.00    0.04
> system.time(XX50<-sapply(trees[1:50],getStates))
  user  system elapsed
  0.28    0.00    0.29
> system.time(XX100<-sapply(trees[1:100],getStates))
  user  system elapsed
  1.16    0.00    1.16
> system.time(XX200<-sapply(trees,getStates))
  user  system elapsed
    6.6    0.0    6.6

Hmmm. What's going on here?

What I discovered (somehow - I'm not sure why I tried this) is that if I first split my list of 200 trees into, say, 20 lists of 10 trees; and then I ran sapply(...,getStates) on each of these lists; then recombined the results using cbind, this is much faster. So, for instance:

> g<-function(trees){
 ff<-as.factor(ceiling(1:length(trees)/10))
 aa<-lapply(split(trees,ff),function(x)   sapply(x,getStates))
 y<-if(length(trees)>10) aa[[1]] else aa
 for(i in 2:length(aa)) y<-cbind(y,aa[[i]])
 y
}
> # now run it
> system.time(YY200<-g(trees))
  user  system elapsed
  0.97    0.00    0.97
> # check all equal
> dim(YY200)
[1] 199 200
> dim(XX200)
[1] 199 200
> all(XX200==YY200)
[1] TRUE

Seriously - what's going on here? (I have some ideas - but I'm not sure. Feedback welcome.)

Here's a new version of getStates with this hack implemented internally:

# function to get node states from simmap style trees
# written by Liam J. Revell 2013
getStates<-function(tree,type=c("nodes","tips")){
  type<-type[1]
  if(class(tree)=="multiPhylo"){
    ff<-as.factor(ceiling(1:length(tree)/10))
    aa<-lapply(split(tree,ff),function(x)
     sapply(x,getStates))
    y<-if(length(tree)>10) aa[[1]] else aa
    for(i in 2:length(aa)) y<-cbind(y,aa[[i]])
  } else if(class(tree)=="phylo"){
    if(type=="nodes"){
      y<-setNames(sapply(tree$maps,function(x)
       names(x)[1]),tree$edge[,1])
      y<-y[as.character(length(tree$tip)+1:tree$Nnode)]
    } else if(type=="tips"){
      y<-setNames(sapply(tree$maps,function(x)
       names(x)[length(x)]),tree$edge[,2])
      y<-setNames(y[as.character(1:length(tree$tip))],
       tree$tip)
    }
  } else stop("tree should be an object of class 'phylo' or 'multiPhylo'")
  return(y)
}

Huge speed-up for rerootingMethod & estDiversity

A couple of months ago I posted an extremely simple function - basically a wrapper for the 'ape' function ace - to do marginal ancestral state reconstruction using the re-rooting method of Yang.

Well - this function is quite slow. This is partly because the tree has to be re-rooted at every internal node; but mostly this is slow for a totally unnecessary reason, and that is that by wrapping around ace (or, rather, a very lightly modified version of ace used internally by phytools), at each re-rooting the function also re-estimates the transition matrix, Q. Obviously - since only symmetric transition matrices are permitted by this method - this is totally unnecessary.

This is now fixed in the latest minor phytools build (phytools 0.2-56) and the result is an enormous speed-up in computation time. So, for instance:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.55’
>
> # simulate
> tree<-pbtree(n=100,scale=1)
> Q<-matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3)
> rownames(Q)<-colnames(Q)<-letters[1:3]
> x<-sim.history(tree,Q)$states
>
> # ok, now estimate using the old version
> system.time(XX<-rerootingMethod(tree,x))
  user  system elapsed
  28.34    0.00  28.42
> # unload phytools and install the new version
> detach("package:phytools",unload=TRUE)
> install.packages("phytools_0.2-56.tar.gz",type="source", repos=NULL)
* installing *source* package 'phytools' ...
** R
...
* DONE (phytools)
> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.56’
>
> # ok, now repeate the analysis using the new version
> system.time(YY<-rerootingMethod(tree,x))
  user  system elapsed
  3.36    0.01    3.38
> plot(XX$marginal.anc,YY$marginal.anc,xlab="marginal ASRs old version",ylab="marginal ASRs new version")

Cool.

The same speed-up can also be applied to estDiversity - which estimates historical lineage diversity at all the nodes of the tree based on the approach of Mahler et al. (2010) (e.g., 1, 2). So, for instance:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.55’
>
> system.time(d.old<-estDiversity(tree,x))
Please wait. . . . Warning - this may take a while!
Completed 10 nodes
Completed 20 nodes
Completed 30 nodes
Completed 40 nodes
  user  system elapsed
 228.64    0.14  236.39
>
> detach("package:phytools",unload=TRUE)
> install.packages("phytools_0.2-56.tar.gz",type="source", repos=NULL)
...
* DONE (phytools)
> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.56’
>
> system.time(d.new<-estDiversity(tree,x))
Please wait. . . . Warning - this may take a while!
Completed 10 nodes
Completed 20 nodes
Completed 30 nodes
Completed 40 nodes
  user  system elapsed
  26.02    0.04  26.41
> plot(d.old,d.new,xlab="estimated historical diversity (old)",ylab="estimated historical diversity (new)")

Wow. That's an enormous difference. Cool.

Bug fix in fastAnc

A couple of days ago a phytools user reported a bug in geomorph (which reverse depends on phytools) that seemed to be due to a problem with fastAnc, which is used internally.

Well, I finally got around to looking into it. It turns out that there was a bug in fastAnc, but I'd overlooked it for good reason - it only occurs under the somewhat idiosyncratic circumstances of when is.binary.tree(tree)=FALSE and the user-specified option vars=FALSE. This is because fastAnc works by re-rooting the tree at all internal nodes of a binary tree and computing the PIC ancestral state & variance. If the input tree is not binary, then it uses multi2di, but then has to back-translate to the original tree which it does using phytools matchNodes.

The problem arose because:

if(!is.binary.tree(tree)){
  ancNames<-matchNodes(tree,btree)
  anc<-anc[as.character(ancNames[,2])]
  names(anc)<-ancNames[,1]
  if(vars) v[as.character(ancNames[,2])]
  names(v)<-ancNames[,1]
}

should have been:

if(!is.binary.tree(tree)){
  ancNames<-matchNodes(tree,btree)
  anc<-anc[as.character(ancNames[,2])]
  names(anc)<-ancNames[,1]
  if(vars||CI){
    v[as.character(ancNames[,2])]
    names(v)<-ancNames[,1]
  }
}

in the code that is executed to back-translate nodes betweeen trees. Basically, if is.binary.tree(tree)=FALSE and vars=FALSE (but under no other circumstances), fastAnc will try to assign names to a vector of zero length. Oops.

The fixed source code for this very simple function is here, but I also posted a new phytools build (phytools 0.2-55).