Thursday, April 24, 2014

New print & plot methods for describe.simmap

I just added new print & plot methods to phytools for the object returned by describe.simmap, a function designed to summarize the results of stochastic mapping conducted with the phytools function make.simmap. Here's a quick demo:

> require(phytools)
Loading required package: phytools
Loading required package: ape
Loading required package: maps
Loading required package: rgl
> packageVersion("phytools")
[1] ‘0.4.2’
> data(anoletree)
> x<-getStates(anoletree,"tips")
> trees<-make.simmap(anoletree,x,nsim=100,model="ER",
  message=FALSE)
> obj<-describe.simmap(trees)
> ## print method
> obj
100 trees with a mapped discrete character with states:
 CG, GB, TC, TG, Tr, Tw

trees have 25.88 changes between states on average

changes are of the following types:
     CG,GB CG,TC CG,TG CG,Tr CG,Tw GB,CG GB,TC GB,TG GB,Tr
x->y  0.27  0.39  0.29  0.18  0.37  0.53  0.77  0.96  0.39
     GB,Tw TC,CG TC,GB TC,TG TC,Tr TC,Tw TG,CG TG,GB TG,TC
x->y  1.09  1.23  0.74  0.51  0.71  0.94  1.15  3.32  1.89
     TG,Tr TG,Tw Tr,CG Tr,GB Tr,TC Tr,TG Tr,Tw Tw,CG Tw,GB
x->y  1.05  1.78  0.21  0.29  0.21  0.21  0.28  0.86   1.8
     Tw,TC Tw,TG Tw,Tr
x->y  1.82  0.97  0.67

mean total time spent in each state is:
             CG         GB         TC        TG          Tr
raw  12.9528230 44.2907143 32.4717543 69.734876 12.85247436
prop  0.0629796  0.2153516  0.1578851  0.339067  0.06249168
            Tw   total
raw  33.364318 205.667
prop  0.162225   1.000

> ## plot method
> plot(obj,fsize=0.6,cex=c(0.6,0.3),ftype="i")
> ## add legend
> states<-sort(unique(getStates(trees[[1]],"tips")))
> add.simmap.legend(colors=
  setNames(palette()[1:length(states)],states),fsize=0.8)
Click where you want to draw the legend

That's not bad. The plotted pies at nodes (and tips) give the posterior probability of each node being in each state from (in this case) empirical Bayesian stochastic character mapping.

The code for these new methods, and the updated describe.simmap, is here; and a new phytools build with these methods can be downloaded from the phytools page.

1 comment:

  1. How to avoid the circles at the tips from being on top of the tip labels? I thought that was what the option offset did, but it's not working. E.g.,

    plot(tr, offset=1)

    gives exactly the same plot as

    plot(tr, offset=1000)

    Thanks!

    ReplyDelete