Monday, August 15, 2011

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.

No comments:

Post a Comment