# moreover if you look at the numbers closely, a <- c(0.998, 2.022, 3.011, 3.993, 5.006, 5.983, 6.994, 8.017, 8.995, 9.984) b <- 1:10 k1 <- cor.test(a, b, method="pearson") k1$p.value # value is 0 k2 <- cor.test(a, -b, method="pearson") k2$p.value # p value is 1.7xe-20 ## The two numbers should be the sample a <- c(1.059,1.980,2.932,4.042,5.015,5.892,7.019,7.954,8.990,10.032) b <- 1:10 k1 <- cor.test(a, b, method="pearson") k1$p.value # 2.220446e-15 k2 <- cor.test(a, -b, method="pearson") k2$p.value # 2.272542e-15 # The correct value should be pt(k1$statistic, k1$parameter, lower.tail = F)*2 # 2.272542e-15 pt(k2$statistic, k2$parameter, lower.tail = T)*2 # 2.272542e-15 # Therefore, when the correlation is positive, the calculation has an incorrect lower bound and the calculation is not correct. # Possible reason: 2.2e-16 is the value of .Machine$double.eps - the smallest floating point number such that 1 + x != 1 # Is it because in your code somewhere you use 1-pvalue? (1-pt(k2$statistic, k2$parameter, lower.tail = F))*2 # 2.220446e-15 # which is the p-value calculated by the current function cor.test # Thank you for your attention. # Best wishes! # any problem? Feel free to contact lichuan@umich.edu

indeed, the small accuracy losses happen because 1 - pt() or 1 - pnorm() are used instead of p*(.., lower.tail=FALSE). I'll commit the fixes within a day or so.