Sunday, November 27, 2011

Simulating an irreversible discrete trait using sim.history()

A user asked how one should go about simulating a character that evolves from state "0" (say) to state "1", but not the reverse. This corresponds to the following transition matrix:

> Q
   0 1
0 -1 0
1  1 0


in which we can substitute any instantaneous transition rate for the 1.0 and -1.0. Unfortunately, for various reasons sim.history doesn't work for transition rates of exactly zero - however it works fine if we substitute a value slightly greater than zero. To see this work, try the following example code that uses "geiger" and "phytools":

require(phytools)
require(geiger)
tree<-drop.tip(birthdeath.tree(b=1,d=0,taxa.stop=101),"101")
tree<-rescaleTree(tree,1)
Q<-matrix(c(-1,1,1e-10,-1e-10),2,2)
rownames(Q)<-colnames(Q)<-c("0","1")
mtree<-sim.history(tree,Q,"0")
cols<-c("blue","red"); names(cols)<-c(0,1)
plotSimmap(mtree,cols,ftype="off")


which should produce something that looks like this:



Don't forget, of course, that the tip states are stored in mtree$states, i.e.:

> mtree$states
 53  98  99  24  25  26  18  13  30  31  56  80 
"0" "0" "0" "0" "0" "0" "0" "0" "1" "1" "1" "1" 
 81  77  ...
"1" "1"  ...

No comments:

Post a Comment