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
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.