Saturday, March 23, 2013

New version of describe.simmap with informative message

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

Code is here - plus new phytools build (phytools 0.2-32) to install from source.

6 comments:

  1. Hello Liam,

    This 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

    ReplyDelete
    Replies
    1. This summarizes the transition rates but seems to stop short of summarizing the outcome on a consensus tree?

      Delete
    2. This can be done. Check out my previous post on the subject.

      Delete