EDIT: see warning below about resolution of Tukey p-values!!
dd <- data.frame(y=c(1:10,1001:1010),f=rep(c("A","B"),each=10))
fit <- aov(y~f,data=dd)
The printed p-value is zero:
(tt <- TukeyHSD(fit))
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = y ~ f, data = dd)
##
## $f
## diff lwr upr p adj
## B-A 1000 997.1553 1002.845 0
But looking at the (abbreviated) output of str()
shows there's more information there:
str(tt)
## List of 1
## $ f: num [1, 1:4] 1.00e+03 9.97e+02 1.00e+03 2.62e-14
## ..- attr(*, "dimnames")=List of 2
##
You can extract the value yourself:
tt$f[,"p adj"]
## [1] 2.620126e-14
Or as noted in the comments, print(tt,digits=15)
will work ...
WARNING
I decided to dig a little deeper and noticed in digging through the code of TukeyHSD.aov()
that it relies on ptukey()
, which in its "Examples" section warns that "the precision may not be more than about 8 digits". In particular, once the t-statistic is above about 30, the p-value maxes out (mins out?) at 2.62e-14
...
zval <- 10^seq(1,6,length=100)
pval <- ptukey(zval,2,18,lower.
par(las=1,bty="l")
plot(zval,pval,log="xy",type="l")
The bottom line is that you can't distinguish among p-values this small at all. You may need to rethink your strategy ...