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
[1] 5.441812
...
$logL1
[1] -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)
[1] 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.)