## Tuesday, October 11, 2011

### Estimating a single rate on the tree with brownie.lite()

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:

> require(phytools)
> 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():

> brownie.lite(tree,x,maxit=10000)
\$sig2.single
 5.441812
...
\$logL1
 -347.2953
...

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:

> mean(pic(x,tree)^2)*(length(tree\$tip)-1)/length(tree\$tip)
 5.441812

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.)