Update to plotSimmap()

In playing around with read.simmap(), I also identified a small issue with the stochastic map plotting function, plotSimmap(). That is, for trees with very short total length (as might be the case with a very slow evolving gene - for instance), the plotted tree is left-aligned in the plotting window and a large amount of white space can be left to right of the tip labels.

To illustrate this, I can first simulate the tree & data:

> require(phytools); require(geiger)
> set.seed(1)
> tree<-rbdtree(b=30,d=0,Tmax=0.1)
> x<-sim.char(tree,model.matrix=list(matrix(c(-10,10,10,-10),2,2)), model="discrete")[,,1]
> mtree<-make.simmap(tree,x)
> cols<-c("blue","red"); names(cols)<-c(1,2)
> plotSimmap(mtree,fsize=0.7,cols)

Which produces the following result:



(border added for illustrative effect).

I had the hardest time trying to resolve this issue; however I finally realized that if I first rescaled the sum of the total tree length and the maximum string width of the longest taxon name (which can be found using strwidth(), once a plotting window has been opened) to 1.0, then the problem is solved.

Let's load the new version of plotSimmap() and try again:

> source("http://faculty.umb.edu/liam.revell/phytools/plotSimmap/v0.6/plotSimmap.R")
> plotSimmap(mtree,fsize=0.7,cols)



Indeed, this looks much better.

Direct link to code for this function is here. I have only tried this on one system - so I would be happy to hear reports that it works on other machines (particularly Mac OS or Linux).

Discrepancy between read.simmap() and SIMMAP

Sam Price yesterday pointed out to me that my functions read.simmap() and write.simmap() assume the opposite order of the edge states along internal branches created by Jonathan Bollback's program SIMMAP. Yikes! This was news to me, but evidently the mapping ":{1,0.10:0,0.90}" in SIMMAP indicates that the branch started in state 0 and then changed to state 1. Quelle surprise! I had assumed the contrary - that the branch, in this case, started in state 1 (in which it remained for 0.10 units of time) before changing to state 2 for the remainder of the edge.

The good news is that misordering the states along the branches is irrelevant for BM based analyses, such as those implemented in brownie.lite() and evol.vcv(). In fact, in the earliest versions of read.simmap(), I discarded this information entirely. However, trees misread in this way look extremely peculiar when plotted using plotSimmap(). For instance, consider the stochastic map below which was created by a colleague using SIMMAPv1.0 and read into R using read.simmap() from the present version of "phytools." Oops!



(In case it's not immediately apparent what is wrong with this image, close examination of almost any branch with a transition shows that this transition occurs at the node, not along the edge as it should.)

Fortunately, this can easily be solved and I have added it to the latest version of read.simmap(). For backward compatibility, and because I think this ordering makes more sense anyway (no offense intended to JB), I have also retained the option of reading in trees in which the mapped states along an edge are order from root to tip in a left-to-right fashion. This alternative format (still produced by write.simmap() until I fix that too) can be read in by setting the option rev.order=FALSE (default is TRUE).

The fix was very easy to code. Basically, I apply it after reading the tree as before in the following way:

tree$maps<-lapply(tree$maps,function(x) x<-x[length(x):1])

All this does is goes through the list of mapped edges and reverses the ordering of the states, and time spent in each state, for each edge.

This can also be used to convert any tree that has already been read into memory, or even to modify a tree in memory so that it can be written with write.simmap() and then read into another program that accepts standard Newick style SIMMAP format trees.

The new version of read.simmap() is available from my webpage (direct link to code here). I will also put this in the next version of "phytools."

New function: estDiversity()

I just posted a "new" function that I actually wrote a while ago for Luke Mahler, and just wanted to return to working on.

Given a tree and a vector, x, containing the biogeographic regions of the tip species, this function computes a vector containing point estimates for the lineage density in each region at each internal node. This is accomplished by first computing the conditional probabilities of being in each region at each internal node (these should probably be the marginal probabilities, but let's just ignore that for a moment). The function then proceeds to take a slice through the tree at the height of that node and compute the marginal probabilities at each "pseudo-node" that is created by the slice bisecting an edge of the tree. The historical lineage density for the node is then estimating by summing these marginal probabilities across the pseudo-nodes and multiplying them by the conditional probabilities at the focal node. This is similar to, but not the same as, the method used in Mahler et al. (2010). Unfortunately, the present implementation requires many calls of the "ape" function ace(), which makes it quite slow. Hopefully, I can improve on this - even if just a little.

Oh yeah - the function is on my R-phylogenetic page (and will also most likely be in the next version of "phytools"). Direct link to the code is here.

Evolutionary rates paper now available online

I just discovered that my paper (with Mahler, Peres-Neto, and Redelings) on the MCMC method for identifying the location of a shift in the evolutionary rate through time is now available online as an "Accepted Article" in Evolution. I have described my work on this method in a number of prior posts (1, 2, 3 ,4, 5, 6, 7, 8, 9, 10). The link to the article is here; however Wiley did not post the supplementary appendix with the article so I will try and post it on my website soon. Check it out!

New version of "phytools" (v0.0-6)

I just posted a new alpha version of "phytools," version 0.0-6, to my website here.

Here is a list of updates or additions from the previous version:

1. New function map.overlap() that computes the fractional overlap between two stochastically mapped phylogenies (described here).

2. More robust version of brownie.lite() that avoids the problem of "false convergence" in the multi-rate model (described here).

3. New and improved versions of plotSimmap() and plotTree(). These versions have more plotting options including the option to plot without tip labels.

4. Fixed lineage-through-time function ltt(). Bug & fix described here & here, respectively.

Please check it out! Feedback welcome.

Error fixed (I hope!)

I believe I have fixed the error with ltt(), identified here. Basically, the problem seems to have been that, after computing the height above the root of all the nodes and tips in the tree, I was incorrectly counting the number of "events" (that is, increments and decrements of the number of lineages in the tree) which should be the number of nodes + the number of extinct species (that is, assuming that two species never go extinct at exactly the same time, to the extent of numerical precision)+ two (the start and the end). I have now fixed this, which I hope has resolved Rob's problem. Hopefully, the fix has not created other problems, to be discovered.

Note that this function also computes the γ statistic of Pybus & Harvey (2000), but I'm not sure of the validity of this for species with extinct taxa.

Direct link to updated version is here. I will also add this to the next version of "phytools" (which I hope to post on my website today or tomorrow).

Error with ltt()

Rob Lanfear has identified a problem with the function ltt() in "phytools" when the tree includes extinct lineages. For instance, take the following tree with 20 extant taxa:

With the present version of ltt(), we get:

which is obviously incorrect. We can substituted v0.3 of the function and try again:

> source("http://faculty.umb.edu/liam.revell/phytools/ltt/v0.3/ltt.R")
> ltt(tree,log.lineages=F)

and this seems to solve the problem. I will try and fix the most recent version (which corrected some bugs from v0.3 as well as adding additional functionality) and this will be included in the next version of "phytools." Thanks to Rob for finding this problem.

Note that the tree above was simulated using Tanja Stadler's neat package "TreeSim."