Friday, June 24, 2011

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.

4 comments:

  1. Direct link to new version of brownie.lite() is here.

    ReplyDelete
  2. Hi Liam,

    regarding the function make.simmap, I can see in the output you get how your states are distributed across the branches of your tree (with $mapped.edge). How difficult could be to visualize exactly the same but for your nodes? I mean, how your states are distributed across the nodes of your tree and with what probability.

    I am new in R and I would appreciate any possible help with this.

    thanks a lot!

    Joan.

    ReplyDelete
  3. Hi Joan.

    I simple function that pulls out the states at the internal (& terminal nodes of the tree is as follows:

    getNodes<-function(tree){
    X<-matrix(NA,nrow(tree$edge),ncol(tree$edge))
    for(i in 1:nrow(tree$edge))
    X[i,]<-names(tree$maps[[i]])[c(1,length(tree$maps[[i]]))]
    return(X)
    }

    The function returns a matrix where the elements correspond to the nodes in matrix tree$edge. Getting, say, the frequencies of nodes in each state across a set of mapped trees is only a little bit more complicated. I will try to write a function & post for this as soon as I get a chance.

    ReplyDelete
  4. thanks a lot Liam, and congratulations for this very interesting blog!

    Joan.

    ReplyDelete