Friday, November 4, 2016

S3 density method for objects of class "multiSimmap"

I just pushed some new methods S3 methods for objects of class "multiSimmap" - that is a set of stochastic character mapped trees. In particular, I created a density method for computing the posterior distribution of changes on the tree of each type.

This is a job in progress - so far I have implemented this only for binary characters (they don't need to be coded 0 / 1, but not more than two states are allowed). There is also no documentation; however because they are generic methods, R will let us build & install the package anyway.

Here is a quick demo:

library(phytools)
## Loading required package: ape
## Loading required package: maps
packageVersion("phytools")
## [1] '0.5.57'
tree
## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.
x
## A B C D E F G H I J K L M N O P Q R S T U V W X Y Z 
## 0 0 0 0 0 0 0 1 0 1 1 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 
## Levels: 0 1
trees<-make.simmap(tree,x,model="ER",nsim=200)
## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##            0          1
## 0 -0.9946799  0.9946799
## 1  0.9946799 -0.9946799
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##   0   1 
## 0.5 0.5
## Done.
trees
## 200 phylogenetic trees with mapped discrete characters
obj<-density(trees)
## Loading required package: coda
obj
## 
## Distribution of changes from stochasting mapping:
##  0->1            1->0
##  Min.   :3       Min.   :0
##  Median :7       Median :3
##  Mean   :7.31    Mean   :3.45
##  Max.   :11      Max.   :13
## 
## 95% HPD interval(0->1): [4, 8]
## 95% HPD interval(1->0): [0, 8]
plot(obj)

plot of chunk unnamed-chunk-1

The density method for objects of class "multiSimmap" also wraps the standard densityMap method. For instance:

obj<-density(trees,method="densityMap",res=200)
## sorry - this might take a while; please be patient
obj
## Object of class "densityMap" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
## 
## (2) The mapped posterior density of a discrete binary character with states (0, 1).
plot(obj,lwd=6,outline=TRUE)

plot of chunk unnamed-chunk-2

This is basically what I have right now. More to come later!

No comments:

Post a Comment