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.

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

ReplyDeleteHi Liam,

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

Hi Joan.

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

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

ReplyDeleteJoan.