Thursday, January 26, 2012

Getting tip states from sim.history

A user just emailed me to ask how you get the tip states from a character history simulated using sim.history. This is pretty easy because sim.history stores them in the vector $states.

So, for one tree:

> require(phytools)
> tree<-pbtree(n=25,scale=1)
> Q<-matrix(c(-1,1,1,-1),2,2)
> mtree<-sim.history(tree,Q)
> mtree$states
t1  t2  t3  t4  t5  t6  t7  t8  t9 t10 t11 t12 t13 t14 t15 t16
"2" "1" "2" "2" "2" "2" "2" "1" "2" "2" "2" "2" "2" "2" "1" "2"
t17 t18 t19 t20 t21 t22 t23 t24 t25
"2" "2" "2" "1" "2" "2" "2" "2" "2"


Or, for multiple trees:

> mtrees<-sim.history(tree,Q,nsim=10)
> mtrees[[1]]$states # first simulated history
t1  t2  t3  t4  t5  t6  t7  t8  t9 t10 t11 t12 t13 t14 t15 t16
"1" "2" "1" "2" "2" "2" "2" "2" "1" "1" "2" "1" "1" "1" "1" "1"
t17 t18 t19 t20 t21 t22 t23 t24 t25
"1" "1" "1" "1" "2" "1" "1" "1" "1"
> mtrees[[2]]$states # second history
t1  t2  t3  t4  t5  t6  t7  t8  t9 t10 t11 t12 t13 t14 t15 t16
"2" "1" "1" "2" "2" "2" "1" "1" "1" "2" "1" "1" "1" "2" "2" "2"
t17 t18 t19 t20 t21 t22 t23 t24 t25
"1" "1" "1" "2" "1" "1" "1" "1" "1"
> # etc, or to get all states from all simulated histories (in columns)
> X<-sapply(mtrees,function(x) x$states) # all histories
> X
    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
t1  "1"  "2"  "2"  "1"  "2"  "1"  "2"  "1"  "1"  "2"  
t2  "2"  "1"  "1"  "2"  "2"  "2"  "1"  "1"  "2"  "2"  
t3  "1"  "1"  "2"  "2"  "1"  "1"  "2"  "2"  "2"  "1"
...


That's it!

No comments:

Post a Comment