Running brownie.lite() for an arbitrarily specified set of branches

A colleague asked me today if it was possible to run brownie.lite() on an arbitrarily specified set of rate regimes, rather than on the set of regimes obtained by mapping a discrete character on the phylogeny (say, using the technique of stochastic character mapping; Huelsenbeck et al. 2003). The simpler answer is "yes" - but it is a tiny bit complicated. brownie.lite() still needs a modified "phylo" object, as created by read.simmap() or make.simmap(), but we can actually build that object ourselves without too much difficulty.

To understand how to do this, we should recall the structure of the simmap "phylo" object in memory. I have blogged about this previously (e.g., here). The modified object contains two additional elements: $mapped.edge, which is a matrix containing the total time spent in each state on each edge, and $maps, which contains both the times and sequences of states on each edge. brownie.lite() uses only $mapped.edge, so we should focus on that.

Let's first consider the simplest case in which each branch, and no fractions of branches, are assigned to specific rate regimes. Say we have a list of branches in each state, where branches are indicated by the node numbers with which these end. Node numbers are the numbers given in the "phylo" object element $edge. For this example, let's say the branches in state "0" are contained in nodes0 and the branches in state "1" in nodes1. Our phylogenetic tree is tree.

First, create the element mapped.edge:

> tree$mapped.edge<-matrix(0,nrow(tree$edge),2, dimnames=list(edge=apply(tree$edge,1,function(x) paste(x,collapse=",")),state=c("0","1")))

Now fill in the matrix:

> tree$mapped.edge[match(nodes0,tree$edge[,2]),"0"]<- tree$edge.length[match(nodes0,tree$edge[,2])]
> tree$mapped.edge[match(nodes1,tree$edge[,2]),"1"]<- tree$edge.length[match(nodes1,tree$edge[,2])]

That's it. Now we have an object, tree, suitable for brownie.lite().