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