Thursday, December 1, 2011

Testing for phylogenetic signal (K) different than 1.0

I can see as administrator of this blog some of the google search terms that lead users to the site. Today, someone ending up on the phytools blog with the following search: "phylogenetic signal k significantly different than one."

Testing for whether K (Blomberg et al. 2003) is significantly different from 1.0 is not explicitly implemented in my function phylosig(), but it is straightforward to imagine how one would go about doing it with simulation.

First, let's simulate data with strong phylogenetic signal (expected K of 1.0):

> require(phytools)
> tree<-rtree(100)
> x<-fastBM(tree)
> K<-phylosig(tree,x)
> K
[1] 1.103667

Now, let's get a null distribution for K=1.0 via simulation on our phylogeny:

> nullK<-apply(fastBM(tree,n=1000,sig2=mean(pic(x,tree)^2)),2, phylosig,tree=tree)

Finally, let's count the number of simulated K values that are more extreme than our observed value for K. I would suggest that we do this by counting the number of times that the absolute value of the logarithm of our observed value for K is no smaller than the absolute values of our simulated values for K. This just makes K=1.5 & K=0.6667 equivalent departures from BM. I.e.,

> P<-mean(abs(log(c(K,nullK)))>=abs(log(K)))
> P
[1] 0.8011988

Ok, now let's contrast this with data generated with no signal.

> x<-rnorm(n=length(tree$tip)); names(x)<-tree$tip.label
> K<-phylosig(tree,x)
> K
[1] 0.1998674
> nullK<-apply(fastBM(tree,n=1000,sig2=mean(pic(x,tree)^2)),2, phylosig,tree=tree)
> P<-mean(abs(log(c(K,nullK)))>abs(log(K)))
> P
[1] 0.000999001

Makes sense.


  1. There is a little trick in here that I didn't mention. For calculation, R treats logical values TRUE and FALSE as 1 and 0, respectively. Thus:

    > x<-c(TRUE,TRUE,FALSE)
    > x
    > 2*x
    [1] 2 2 0
    > sum(x)
    [1] 2
    > mean(x)
    [1] 0.6666667

  2. Thanks, Liam. Very helpful post!