Thursday, September 22, 2011

New function: sim.history()

I just posted a new function that simulates stochastic character histories for a discretely valued character trait. This function is called sim.history(), and a directly link to the source is here. The analysis is highly similar to stochastic character mapping (e.g., here), except that the simulation is unconstrained by actual data. Unlike the "geiger" function sim.char(), which can also do discrete character simulation, sim.history() simulates not only the states at the tips of the tree, but also the states at all the internal nodes of the tree as well as the timing of all the simulated transitions between states. The format for storing this information in memory is the same as used by read.simmap() (>=v0.2) and make.simmap().

To simulate, we first need our instantaneous transition matrix, Q. Let's make Q symmetric, but with uneven rates among states. For fun, let's also call our three states "blue", "green", and "red".

> source("sim.history.R")
> Q<-matrix(c(-1,0.8,0.2,0.8,-1.2,0.4,0.2,0.4,-0.6),3,3, dimnames=list(c("red","green","blue"),c("red","green","blue")))
> Q
red green blue
red -1.0 0.8 0.2
green 0.8 -1.2 0.4
blue 0.2 0.4 -0.6

Now let's get a stochastic birth-death tree:

> require(ape)
> tree<-rbdtree(birth=1.0,death=0.0,Tmax=log(25/2))

Finally, let's simulate on the tree & plot:

> mtree<-sim.history(tree,Q)
> cols<-rownames(Q); names(cols)<-rownames(Q)
> plotSimmap(mtree,cols)

And we get the following (obviously individual results will vary):

Both the node & tip states are also stored in the modified "phylo" object: here as $node.states and $states. I include this mainly to facilitate some analyses that I have been working on.

This function should be added to the next version of "phytools".


  1. Hi, Liam,

    Thanks so much for your post!

    I have been trying to use an ultrametric tree resulting from Phylobayes (I made sure all tips have the same distance from the root) as input for make.simmap and sim.history. However, for some reason I get an error caused by a wrong Q-matrix. force.ultrametric() does not work on my tree, only chronos(), however this last function changes the branches length dramatically, and I am interested in the timing of the stochastic mapping results.

    Any advice?

    Thanks so much in advance!


  2. Just to follow up, this is the error I keep getting:
    Error in if (dt[2] < tree$edge.length[j]) new.state <- rstate(Q[, state][-match(state, :
    missing value where TRUE/FALSE needed
    In addition: Warning message:
    In rexp(n = 1, rate = -Q[state, state]) : NAs produced


Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.