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:

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] 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 χ^{2} 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:

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.

## No comments:

## Post a Comment