A user just contacted me to inquire as to how he could obtain the MLE of the single evolutionary rate for a character on a phylogenetic tree. This can easily be done with brownie.lite(), but it involves 1 simple extra step. To describe this, let me illustrate it with an example. First, let's load "phytools" and simulate a random tree:
> tree<-rtree(150) # from "ape"
Now, generate random data on the tree using fastBM():
>x<-fastBM(tree,sig2=5) # evolutionary rate is sig^2=5
OK, now the trick. brownie.lite() is designed to fit multiple evolutionary rates to the data & tree, and it reads the rate regimes from the stochastic map style trees created by read.simmap() & make.simmap(). The problem here is that we only have one rate regime. We can solve this easily by appending the matrix, $mapped.edge to our "phylo" object tree. We do this as follows:
> tree$mapped.edge<-matrix(tree$edge.length,nrow(tree$edge), 1,dimnames=list(NULL,"0"))
Now we can run brownie.lite():
Note that if all we want is the MLE of the rate, and don't care at all about the likelihood, we can skip a lot of these steps & just do:
which gives the same rate! (Note that we multiply by n-1/n to get the ML rate, but the mean square of the contrasts is actually unbiased. See O'Meara et al. 2006.)
Post a Comment
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.