I have three vectors of ~8000 genes, each of which has an associated mutation frequency, e.g.:
alpha = c(0.84, 0.87, 0.91...)
beta = c(0.97, 0.94, 0.99...)
kappa = c(0.72, 0.68, 0.75...)
I am using the t.test function in R to compute a P-value between these vectors. Because of how many genes there are and how consistent the differences between each vector, my P-values are always super low. For example, when I run:
ab = t.test(x = alpha, y = beta)
ak = t.test(x = alpha, y = kappa)
bk = t.test(x = beta, y = kappa)
The summaries of all tests output the same value 2.2e-16, which I don't like. When I extract the P-value directly:
t.test(alpha, kappa)$p.value
R simply outputs a 0. I'm assuming this is due to some internal float/double size limit, but it doesn't look great. We are publishing this data, and whether we take the summary P-values that are all the same or the 0, neither looks good. Is there a way around this in R such that I can compute arbitrarily low P-values? If another tool would work better, I'd be interested in that too. Thank you!
If your t-statistic is really so enormous that your p-value is underflowing to zero (i.e. the p-value is less than approximately 1e-308 (!!), and thus can't be distinguished from zero in standard 64-bit floating-point precision), you can adapt the machinery from this answer to get the p-value ...
I can't (or won't) help commenting that such extreme p-values are certainly not meaningful in any usual statistical sense (Andrew Gelman gets snarky on his blog about p-values of 10^(-246), which are thousands of orders of magnitude larger than the ones you report in comments ...). While there is some value in conforming to the conventions of your field, it might be better to just report t-statistics, which would give you a reasonable measure of the difference between groups (or just the raw difference in the means), and leave out the (meaningless) p-values?
x <- rep(1:5, 100)
y <- rep(10001:10005, 100)
tt <- t.test(x,y)
Results:
Welch Two Sample t-test
data: x and y
t = -111692, df = 998, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-10000.176 -9999.824
sample estimates:
mean of x mean of y
3 10003
tt$p.value ## 0
pvalue.extreme <- function(t, df) {
log.pvalue <- log(2) + pt(abs(t), df = df,
lower.tail = FALSE, log.p = TRUE)
log10.pvalue <- log.pvalue/log(10) ## from natural log to log10
mantissa <- 10^(log10.pvalue %% 1)
exponent <- log10.pvalue %/% 1
## or return(c(mantissa,exponent))
return(sprintf("p value is %1.2f times 10^(%d)",mantissa,exponent))
}
pvalue.extreme(tt$statistic, tt$parameter)
[1] "p value is 1.11 times 10^(-3543)"