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)
> 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"
$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.
There was a small bug in this - fixed here.
ReplyDelete