3

I'm computing a p-value from an F-test in R, and it seems to have trouble displaying anything lower than 1e-16. E.g., for F values from 80 to 90 with degrees of freedom 1,200: 

> 1-pf(80:90,1,200)
 [1] 2.220446e-16 2.220446e-16 1.110223e-16 1.110223e-16 0.000000e+00 0.000000e+00 0.000000e+00
 [8] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00

How can I increase the precision of the pf() function's calculations?

Stephen Turner
  • 2,574
  • 8
  • 31
  • 44
  • 2
    You mean you want more digit display? Try : options(digits=15) Or maybe you need to know that R uses finite precision arithmetic, your answers aren't accurate beyond 15 or 16 decimal places – ThiS Jan 31 '13 at 16:29
  • 1
    What is the practical purpose of this exercise? If the probability is extremely low, its "exact" value shouldn't matter. However, maybe `pf(80:90,1,200,log.p=TRUE)` is of interest. – Roland Jan 31 '13 at 16:29
  • 1
    Read http://stackoverflow.com/q/11328784/1412059 – Roland Jan 31 '13 at 16:37
  • @Roland: Needed exact precision to very low P. Performing millions of tests, a Bonferroni correction will necessarily mean needing very low p-values to call something statistically significant. – Stephen Turner Jan 31 '13 at 22:41
  • Performing millions of tests sounds dubious, even if you adjust p-values. – Roland Feb 01 '13 at 10:53
  • It's very common in genetics to perform many millions of tests. http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3023815/ – Stephen Turner Feb 01 '13 at 13:28

2 Answers2

12

p-values this low are meaningless anyway. Firstly, most calculations use slight approximations so the imprecision comes to dominate the result as you tend towards a zero p-value and secondly, and probably more importantly, any tiny deviation of your population from the modelled distribution will overwhelm the accuracy you desire.

Simply quote the p-values as 'p < 0.0001' and be done with it.

Jack Aidley
  • 19,439
  • 7
  • 43
  • 70
4

In addition to Jack's practical advice about p-values, the functions are well defined (if not practical) so I'll give the finite-precision mathematical reason this doesn't work.

.Machine$double.eps is 2.220446e-16 and that is the smallest number that you can add to 1 and get something different. So differencing from 1, that is the smallest value you get. 

> pf(80:90,1,200)
 [1] 1 1 1 1 1 1 1 1 1 1 1
> sprintf("%.17f",pf(80:90,1,200))
 [1] "0.99999999999999978" "0.99999999999999978" "0.99999999999999989"
 [4] "0.99999999999999989" "1.00000000000000000" "1.00000000000000000"
 [7] "1.00000000000000000" "1.00000000000000000" "1.00000000000000000"
[10] "1.00000000000000000" "1.00000000000000000"
> sprintf("%a", pf(80:90,1,200))
 [1] "0x1.ffffffffffffep-1" "0x1.ffffffffffffep-1" "0x1.fffffffffffffp-1"
 [4] "0x1.fffffffffffffp-1" "0x1p+0"               "0x1p+0"              
 [7] "0x1p+0"               "0x1p+0"               "0x1p+0"              
[10] "0x1p+0"               "0x1p+0"

However you can use the approximation $1-p = -\ln(p)$ and the fact that you can get the log of the p-values more precisely

> -pf(80:90,1,200,log.p=TRUE)
 [1] 2.540347e-16 1.770938e-16 1.236211e-16 8.640846e-17 6.047690e-17
 [6] 4.238264e-17 2.974043e-17 2.089598e-17 1.470045e-17 1.035491e-17
[11] 7.303070e-18
Brian Diggs
  • 57,757
  • 13
  • 166
  • 188