Small update to phylosig

The following query was recently posted to the bioinformatics question & answer site, BioStar:

"I use phytools to analyse phylogenetic signal in different tree topologies. In likelihood ratio test with likelihoods for Pagel's lambda I recieve p-value=0 in all my topologies (1000 trees). Not 0.000, or something like that, just 0. I understand, that it indicates highly significant LRT result. But I'm in doubt, is it theoretically possible? And what does it mean from mathematical point of view?"

The posted answer (that is that the P-value is low enough to push the limit of its floating point representation in the method as implemented) is pretty much right on the money; however there is a little more nuance in that I for some reason have computed the P-value of the likelihood-ratio using:

1-pchisq(...)

instead of:

pchisq(...,lower.tail=FALSE)

This distinction is inconsequential until pchisq(...) approaches 1 (to a very high degree of numerical precision), at which point 1-pchisq(...) goes to 0. This can be true even if pchisq(...,lower.tail=FALSE) is numerically distinguishable from 0.

So, for instance:

> 1-pchisq(10,df=1)
[1] 0.001565402
> pchisq(10,df=1,lower.tail=FALSE)
[1] 0.001565402
> 1-pchisq(20,df=1)
[1] 7.744216e-06
> pchisq(20,df=1,lower.tail=FALSE)
[1] 7.744216e-06
> # but...
> 1-pchisq(80,df=1)
[1] 0
> pchisq(80,df=1,lower.tail=FALSE)
[1] 3.744097e-19

Does this distinction between P=3.74×10^-19 and P=0 "matter" (that is, even ignoring the fact that the χ² only approximates the probability distribution of the likelihood-ratio anyway)? Well, of course not. Nonetheless, this is easy enough to fix. Here's a link to the updated version. Let's load a new build of phytools containing this update and try it:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.34’
> # simulate tree & data
> tree<-pbtree(n=100)
> x<-fastBM(tree)
> phylosig(tree,x,method="lambda",test=TRUE)
$lambda
[1] 1.000031
$logL
[1] -135.1252
$logL0
[1] -201.9702
$P
[1] 0

> detach("package:phytools",unload=TRUE)
> install.packages("phytools_0.2-35.tar.gz",type="source", repos=NULL)
* installing *source* package 'phytools' ...
...
* DONE (phytools)
> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.35’
> phylosig(tree,x,method="lambda",test=TRUE)
$lambda
[1] 1.000031
$logL
[1] -135.1252
$logL0
[1] -201.9702
$P
[1] 6.387017e-31

That's it.