R small pvalues

Question

I am calculating z-scores to see if a value is far from the mean/median of the distribution. I had originally done it using the mean, then turned these into 2-side pvalues. But now using the median I noticed that there are some Na's in the pvalues.

I determined this is occuring for values that are very far from the median. And looks to be related to the pnorm calculation. " 'qnorm' is based on Wichura's algorithm AS 241 which provides precise results up to about 16 digits. "

Does anyone know a way around this as I would like the very small pvalues. Thanks,

> z<- -12.5
> 2-2*pnorm(abs(z))
[1] 0
> z<- -10  
> 2-2*pnorm(abs(z))
[1] 0 
> z<- -8 
> 2-2*pnorm(abs(z))
[1] 1.332268e-15

Lesbianism

Answer 1

Intermediately, you are actually calculating very high p-values:

options(digits=22)
z <- c(-12.5,-10,-8)
pnorm(abs(z))
# [1] 1.0000000000000000000000 1.0000000000000000000000 0.9999999999999993338662
2-2*pnorm(abs(z))
# [1] 0.000000000000000000000e+00 0.000000000000000000000e+00 1.332267629550187848508e-15

I think you will be better off using the low p-values (close to zero) but I am not good enough at math to know whether the error at close-to-one p-values is in the AS241 algorithm or the floating point storage. Look how nicely the low values show up:

pnorm(z)
# [1] 3.732564298877713761239e-36 7.619853024160526919908e-24 6.220960574271784860433e-16

Keep in mind 1 - pnorm(x) is equivalent to pnorm(-x). So, 2-2*pnorm(abs(x)) is equivalent to 2*(1 - pnorm(abs(x)) is equivalent to 2*pnorm(-abs(x)), so just go with:

2 * pnorm(-abs(z))
#  [1] 7.465128597755427522478e-36 1.523970604832105383982e-23 1.244192114854356972087e-15

which should get more precisely what you are looking for.

Untie

Answer 2

One thought, you'll have to use an exp() with larger precision, but you might be able to use log(p) to get slightly more precision in the tails, otherwise you are effectively at 0 for the non-log p values in terms of the range that can be calculated:

> z<- -12.5
> pnorm(abs(z),log.p=T)
[1] -7.619853e-24

Converting back to the p value doesn't work well, but you could compare on log(p)...

> exp(pnorm(abs(z),log.p=T))
[1] 1

Sapphire

Answer 3

pnorm is a function which gives what P value is based on given x. If You do not specify more arguments, then default distribution is Normal with mean 0, and standart deviation 1.
Based on simetrity, pnorm(a) = 1-pnorm(-a).
In R, if you add positive numbers it will round them. But if you add negative no rounding is done. So using this formula and negative numbers you can calculate needed values.

> pnorm(0.25)
[1] 0.5987063
> 1-pnorm(-0.25)
[1] 0.5987063
> pnorm(20)
[1] 1
> pnorm(-20)
[1] 2.753624e-89

Velarde

To me it appears that this doesn't work for the scenario that is described in the question 2-2*pnorm(abs(z)). – Machute