Thursday, March 21, 2013

Function to summarize the results of stochastic mapping

To address a user request I just quickly wrote a new utility function, describe.simmap, to summary the results from stochastic maps on one or multiple trees. It computes the total number of transitions and the number of each type (using countSimmap; and it computes the posterior probabilities for each node. Finally it can optionally plot those probabilities on a tree. Basically, it pulls together some different things that I've been doing for analyses and visualization in some recent posts (e.g., 1, 2, 3, 4, 5, 6). Code is here.

Here's a demo:

> library(phytools)
> source("describe.simmap.R")
> tree<-pbtree(n=60,scale=1)
> Q<-matrix(c(-3,1,1,1,1,-3,1,1,1,1,-3,1,1,1,1,-3),4,4)
> rownames(Q)<-colnames(Q)<-c("a","c","g","t")
> x<-sim.history(tree,Q)$states
> mtrees<-make.simmap(tree,x,model="ER",nsim=100)
make.simmap is sampling character histories conditioned on the transition matrix
Q =
           a          c          g          t
a -2.8373734  0.9457911  0.9457911  0.9457911
c  0.9457911 -2.8373734  0.9457911  0.9457911
g  0.9457911  0.9457911 -2.8373734  0.9457911
t  0.9457911  0.9457911  0.9457911 -2.8373734
(estimated using likelihood);
and state frequencies (used to sample the root node following Bollback 2006)
pi =
   a    c    g    t
0.25 0.25 0.25 0.25
> XX<-describe.simmap(mtrees,plot=TRUE,cex=0.7)
> XX
$count
       N a,c a,g a,t c,a c,g c,t g,a g,c g,t t,a t,c t,g
 [1,] 64   4   8   6   3   5   1  12   6   5   6   3   5
 [2,] 46   4   3   2   4   5   3   9   6   4   1   3   2
 [3,] 46   4   4   0   3   3   3   4   8   5   4   3   5
 [4,] 55   7   8   3   5   7   3   5   5   1   5   5   1
 [5,] 59 ...

$ace
       a    c    g    t
61  0.20 0.31 0.31 0.18
62  0.24 0.32 0.29 0.15
63  0.12 0.87 0.00 0.01
64  0.12 0.87 0.00 0.01
65  ...

$legend
black    red green3   blue
  "a"    "c"    "g"    "t"

That's it.

1 comment: