Thursday, February 28, 2013

Function to count transitions from a mapped tree

To address a user request I just posted a simple function, countSimmap (code here), to count the number of transitions (in total and by type) on a discrete character mapped tree such as a stochastically mapped (i.e., SIMMAP) tree.

Here's a quick demo:
> require(phytools)
> source("utilities.R")
> Q<-matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3)
> colnames(Q)<-rownames(Q)<-c("A","B","C")
> tree<-sim.history(pbtree(n=10,scale=1),Q)
> cols<-setNames(c("blue","red","green"),colnames(Q))
> par(col="white")
> plotTree(tree,ftype="i",lwd=5)
> par(col="black")
> plotSimmap(tree,cols,pts=F,lwd=3,ftype="i",add=T)
> countSimmap(tree,colnames(Q))
$N
[1] 13

$Tr
  A B C
A 0 6 1
B 4 0 2
C 0 0 0

$message
[1] "N is the total number of character changes on the tree"
[2] "Tr gives the number of transitions from row state->column state"

The function is even fast enough to run on quite large trees with lots of transitions, for instance:
> tree<-sim.history(pbtree(n=1000,scale=10),Q)
> system.time(XX<-countSimmap(tree,colnames(Q)))
  user  system elapsed
  0.22    0.00    0.22
> XX
$N
[1] 3762

$Tr
   A   B   C
A   0 623 646
B 605   0 623
C 656 609   0

$message
[1] "N is the total number of character changes on the tree"
[2] "Tr gives the number of transitions from row state->column state"

That's it.

No comments:

Post a Comment