A phytools user just commented that it would be cool to be able to automate the calculation of Pybus & Harvey's γ statistic for a set of phylogenies (say, a posterior sample of trees from a Bayesian analysis). This is actually already pretty easy - just do the following (trees is a set of trees stored in an object of class "multiPhylo"):
g<-sapply(trees,function(x) ltt(x)$gamma)
For example, let's try it on a set of pure-birth phylogenies simulated using phytools::pbtree:
> library(phytools)
> trees<-pbtree(n=100,nsim=1000,scale=1)
> g<-sapply(trees,function(x) ltt(x)$gamma)
> # gamma should be standard normal for Yule trees
> mean(g)
[1] -0.04430857
> var(g)
[1] 0.9603543
> hist(g,main="Pybus' gamma on pure-birth trees")
Cool.
Running ltt on a set of trees also incidentally creates a neat animation. Here's what it looked like for my 1000 simulated trees (more or less - I added the counter & diagonal reference line for effect):
That's it!
I just realized that the diagonal dashed line (the expectation for a pure-birth tree) should connect log(2) with log(100) [not log(1)=0 with log(100)], since the trees start with two lineages, not one. Oh well.
ReplyDelete