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.

2 comments:

  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
    [1] TRUE TRUE FALSE
    > 2*x
    [1] 2 2 0
    > sum(x)
    [1] 2
    > mean(x)
    [1] 0.6666667

    ReplyDelete
  2. Thanks, Liam. Very helpful post!

    ReplyDelete