Addendum to fast Brownian simulation

Before I move on to the next full post, I wanted to add a quick addendum to the previous one.

That is, that we can use our new fastBM() function (available here) to simulate under models of evolution other than Brownian motion as well. This by using Luke Harmon et al.'s {geiger} tree transformation function deltaTree(). This function transforms the branch lengths of the tree to be proportional in length to the variances and covariances among species under various models. For instance, the function call:

> tr<-ouTree(phy,alpha)

returns a tree with branch lengths proportional to the variances and covariances among species expected under an Ornstein-Uhlenbeck process with constraint ("selection") parameter, alpha. Similarly:

> tr<-speciationalTree(phy)

returns a tree with all edge lengths set to 1.0.

To simulate under these models using the fastBM() function we came up with last week, we just send fastBM() the phylogenetic tree returned by deltaTree(). For example, to simulate under an Ornstein-Uhlenbeck model with the constraint parameter, alpha, set to alpha=3, we just compute:

> x<-fastBM(tree=ouTree(phy,alpha),alpha,sig2)

substituting our desired values for sigma-squared and the root value for sig2 and alpha (the second) in fastBM(), respectively. [Sorry for the confusing double use of "alpha" here, folks. By convention, alpha is the restraining force in an Ornstein-Uhlenbeck model e.g., Butler & King (2004); unfortunately, alpha is also sometimes given as the ancestral state at the root node of the tree. I'll try to remember to fix this in any updates of fastBM().]

Keep in mind, fastBM() has given us all the terminal and internal states, so should we want only the tip states, we need to then compute:

> x<-x[1:length(phy$tip)]

Similarly, we can simulate under the purely speciational model as follows:

> x<-fastBM(tree=speciationalTree(phy),alpha,sig2)

Here, sig2 is the variance on the speciational "jumps" of evolutionary change during punctuational evolution in our clade.

Simulating character evolution this way is also quite fast. This is because the function ouTree() is fast. (I think this is because it requires just two passes of the tree - one to get the branching times and another to rescale all the edges by the OU model.) For instance, using a 200 species stochastic tree as before, we can see that the simulation:

> system.time(x<-fastBM(tree=ouTree(phy=tree,alpha=2),
   sig2=1,alpha=0))
user system elapsed
0.34 0.00 0.34

takes only about 1/3 of a second on my netbook.