Using make.simmap() with brownie.lite()

Now that make.simmap() seems to work, I thought I would give a quick demo of how it can be used in conjunction with fastBM() and brownie.lite(), as well as a couple of functions from the {geiger} package, to simulate and then estimate under O'Meara et al.'s (2006) multi-rate Brownian motion model.

First, we load the dependencies (and source):
> require(phytools)
> require(geiger)
> source("make.simmap.R")

Now, let's simulate a tree and our discrete character:
> tree<-drop.tip(birthdeath.tree(b=1,d=0,taxa.stop=101),"101")
> x<-sim.char(tree,model.matrix=list(matrix(c(-0.2,0.2,0.2,-0.2) ,2,2)),model="discrete")[,,1]

We can create a stochastic character mapped tree for x using make.simmap():
> mtree<-make.simmap(tree,x)

Now, we take advantage of the way this map is stored to simulate on the tree. We'll use a rate of 1.0 for branches in state "1", and rate 10.0 for branches in state "2".
> simtree<-tree
> simtree$edge.length<-mtree$mapped.edge[,"1"]+ mtree$mapped.edge[,"2"]*10
> y<-fastBM(simtree)

Finally, we fit the model using brownie.lite():
> result<-brownie.lite(mtree,y)
> result
$sig2.single
[1] 3.410848
$a.single
...
$logL1
[1] -186.3361
$k1
[1] 2
$sig2.multiple
1 2
0.953248 8.705116
...
$logL.multiple
[1] -161.2660
$k2
[1] 3
$P.chisq
[1] 1.431489e-12
$convergence
[1] "Optimization has converged."

Pretty cool!

P.S., I will shortly post a version of brownie.lite(), that uses box-constraints on the optimization via optim(). This seems to be more stable. It also returns ancestral values for the root node of the tree under each model.