Thursday, May 9, 2013

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.

4 comments:

  1. Hi Liam,
    This is exciting, I've just started playing with densityMap and it's a great visualization tool. The links to the updated function and build don't appear to be active (yet?).
    One suggestion for a future iteration - for print publication figures, people may want a greyscale option instead of the red-blue colors. I've tried it out by tweaking the line of the function that defines the colors to:

    cols <- grey(seq(0,1,length.out=1001))

    and while it probably doesn't have quite the same ability to pick up small differences than the red-blue, it looks pretty good in combination with outline=TRUE.

    ReplyDelete
    Replies
    1. Hi Travis.

      Cool. Sorry about the bum links. Should be fixed.

      What you suggest is very easy now, just do:

      maps<-densityMap(trees)
      maps$cols[]<-grey(seq(0,1,length.out=1001))
      plot(maps,outline=TRUE)

      with any other options that you want.

      I will post this to the blog as I think it is a great idea!

      - Liam

      Delete
    2. Longer post on this here. Thanks for the suggestion Travis. Awesome. Liam

      Delete
    3. Works great, thanks!

      Delete