I just wrote some code for my phylogeny methods class (class materials here) to simulate evolution on the phylogeny by the Ornstein-Uhlenbeck process using a constant θ and α without transforming the tree. I think that what I did is correct. I have posted the code (it is an internally called function of fastBM) here. It is also in a minor phytools version (phytools 0.3-77) which can be downloaded & installed from source.
One of the main reasons it could be interesting to simulate OU using my function - rather than by first transforming the tree using transform.phylo (formerly ouTree) in the geiger package & then simulating on the transformed tree using BM - is because it can be used to simulate θ (the position of the optimum) that is different from the x0, the state at the root node. As many of you probably know - it is not possible to fit this model to a tree with contemporaneous tips (the model is non-identifiable); however simulating under the model presents no problem theoretically.
Here's a quick demo:
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.3.77’
> tree<-pbtree(n=26,tip.label=LETTERS,scale=10)
> x<-fastBM(tree,a=0,theta=3,alpha=0.2,sig2=0.1, internal=TRUE)
> phenogram(tree,x,spread.labels=TRUE,spread.cost= c(1,0.01))
To create the plot above also requires a new version of phenogram, which handles label spreading slightly different (using the full range on the vertical axis, rather than simply the range of tip values for extant species). This is also in the latest phytools version.
No comments:
Post a Comment
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.