I just posted a new version of describe.simmap that prints an informative message and returns the results invisibly (if message is set to TRUE). This can be used on a single tree simulated with sim.history, or on a set of stochastic maps obtained using make.simmap or read in with read.simmap from the stand-alone program SIMMAP. Here's a demo:
> library(phytools)
> # simulate data
> tree<-pbtree(n=100,scale=1)
> Q<-matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3)
> rownames(Q)<-colnames(Q)<-LETTERS[1:3]
> tree<-sim.history(tree,Q,anc="A")
> # conduct stochastic mapping on simulated data
> mtrees<-make.simmap(tree,tree$states,model="ER",nsim=100)
make.simmap is sampling character histories conditioned on the transition matrix
Q =
A B C
A -2.304145 1.152073 1.152073
B 1.152073 -2.304145 1.152073
C 1.152073 1.152073 -2.304145
(estimated using likelihood);
and state frequencies (used to sample the root node following Bollback 2006)
pi =
A B C
0.3333333 0.3333333 0.3333333
> # simulated data
> describe.simmap(tree)
1 tree with a mapped discrete character with states:
A, B, C
tree has 42 changes between states
changes are of the following types:
A B C
A 0 8 6
B 3 0 11
C 9 5 0
mean total time spent in each state is:
A B C total
raw 8.2115414 5.5285739 7.9719994 21.71211
prop 0.3782009 0.2546308 0.3671683 1.00000
> # stochastic mapping results
> describe.simmap(mtrees,plot=TRUE)
100 trees with a mapped discrete character with states:
A, B, C
trees have 55.72 changes between states on average
changes are of the following types:
A,B A,C B,A B,C C,A C,B
x->y 7.43 9.67 8.48 12.57 8.79 8.78
mean total time spent in each state is:
A B C total
raw 7.7147098 6.9697895 7.0276153 21.71211
prop 0.3553182 0.3210092 0.3236725 1.00000
> # simulate data
> tree<-pbtree(n=100,scale=1)
> Q<-matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3)
> rownames(Q)<-colnames(Q)<-LETTERS[1:3]
> tree<-sim.history(tree,Q,anc="A")
> # conduct stochastic mapping on simulated data
> mtrees<-make.simmap(tree,tree$states,model="ER",nsim=100)
make.simmap is sampling character histories conditioned on the transition matrix
Q =
A B C
A -2.304145 1.152073 1.152073
B 1.152073 -2.304145 1.152073
C 1.152073 1.152073 -2.304145
(estimated using likelihood);
and state frequencies (used to sample the root node following Bollback 2006)
pi =
A B C
0.3333333 0.3333333 0.3333333
> # simulated data
> describe.simmap(tree)
1 tree with a mapped discrete character with states:
A, B, C
tree has 42 changes between states
changes are of the following types:
A B C
A 0 8 6
B 3 0 11
C 9 5 0
mean total time spent in each state is:
A B C total
raw 8.2115414 5.5285739 7.9719994 21.71211
prop 0.3782009 0.2546308 0.3671683 1.00000
> # stochastic mapping results
> describe.simmap(mtrees,plot=TRUE)
100 trees with a mapped discrete character with states:
A, B, C
trees have 55.72 changes between states on average
changes are of the following types:
A,B A,C B,A B,C C,A C,B
x->y 7.43 9.67 8.48 12.57 8.79 8.78
mean total time spent in each state is:
A B C total
raw 7.7147098 6.9697895 7.0276153 21.71211
prop 0.3553182 0.3210092 0.3236725 1.00000
Code is here - plus new phytools build (phytools 0.2-32) to install from source.
Hello Liam,
ReplyDeleteThis function is great.
I wonder if it would be possible to apply make.simmap to an object multiphylo (to deal with phylogenetic uncertainty) and to summerize the outcome on a consensus tree.
Best
Julien
Hi Julien.
DeleteCheck out http://blog.phytools.org/2013/04/using-makesimmap-on-set-of-trees.html. Thanks for the feedback!
- Liam
This summarizes the transition rates but seems to stop short of summarizing the outcome on a consensus tree?
DeleteThis can be done. Check out my previous post on the subject.
DeleteChecking it out. Thanks.
DeleteChecking it out. Thanks.
Delete