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