Stochastic character mapping on the tree

I'm just now returning from the 'Evolution' meeting (joint meeting of SSE, ASN, and SSB) in Norman, Oklahoma. I saw many good and exciting talks at the meeting. Hopefully I will be able to channel some of this excitement into new functions and methods for my nascent {phytools} package in the coming weeks.

Rich Fitzjohn, a Ph.D student with Sally Otto at UBC, gave an excellent talk on fitting and using state-dependent diversification models; however, he also began his talk with a very clear and simple explanation of the algorithm for stochastic character mapping, originally developed by Nielsen (2002) and Huelsenbeck et al. (2003). A couple of my functions (in particular, brownie.lite(), which implements the method of O'Meara et al. (2006); and evol.vcv(), which implements Revell & Collar (2009)) take stochastic map trees created by my function read.simmap() as input.

As it turns out, Rich has already programmed stochastic character mapping in his impressive package {diversitree} that I am just beginning to explore; however after Rich's explanation the exercise seemed to be so straightforward that I thought I would try it myself. This way it will be easy for users to seamless go from discrete and continuous character vectors to the rate variation analyses of brownie.lite() and evol.vcv().

To conduct stochastic character mapping we first need the conditional likelihoods of each of our states for our mapped discrete character at each node in the tree. Making the task even easier, it is my understanding that these are returned as the object $lik.anc by the {ape} function ace(...,type="discrete"). [Note, this is not well documented - so I am still trying to confirm that these likelihoods are indeed the "right" ones.] Once we have these, we just traverse the tree once, from the root to the tips. At the root, we just assign the root state with probability determined by the relative conditional likelihoods at each root. Then, for each node, we compute the element-wise product of the vector of conditional likelihoods for that node and the vector P_i·, where P for a branch length t is given by: e^Q*t; and then we stochastically assign internal nodes using the relative probabilities we have computed. (Here, e^X denotes the matrix exponential of X.)

After we have obtained a set of probable states for all the nodes of the tree, we then generate plausible character histories along each edge. We do this by generating waiting times for character shifts from an exponential distribution with rate -Q_ii for starting state i. We do this along each edge of the tree, rejecting branches that don't start and finish with the states simulated in the previous step of our analysis.

I have just posted a beta version of this function to my R phylogenetics page (direct link here). Assuming that there doesn't prove to be any major bugs in the function, I will also incorporate into the next version of {phytools}.

To use the function, simply load the source, and then type:

> mtree<-make.simmap(tree,x)

for "phylo" object tree and discrete character vector x then cross your fingers.* [*This is not a required step, but it can never hurt.] Right now the function fits only a general time-reversible model (i.e., symmetric Q, all rates different); but in future I will add additional options.