# How to calculate Tukey-adjusted p-values for emmeans pairwise comparisons?

I would like to calculate Tukey-adjusted p-values for emmeans pairwise comparisons. I know that these can be obtained directly with functions like pairs() and CLD(). However, when there are three leading zeroes in the p-value, only one digit is displayed. I recognize that in this case the significance of the test statistic is not in question, but I like to consistently report p-values with two digits in papers.

In attempting to calculate a more precise p-value based on the output of pairs(), I have not been able to figure out how to do this when there are multiple comparisons. It's straightforward when there is just one comparison:

> pairs(emmeans(model1, "harvest"), details = T)
contrast              estimate        SE df t.ratio p.value
Spring - Spring/Fall 0.4521333 0.1006861 15   4.491  0.0004

> 2*pt(4.491, 15, lower=FALSE)
[1] 0.0004309609


However, when there are multiple comparisons, I can't figure out how to calculate the appropriate Tukey-adjusted p-value. An unadjusted p-value is too low and an adjusted p-value is too high (using the contrast between factor levels 15 and 61 as an example).

> pairs(emmeans(model2, "row.space"))
contrast  estimate        SE df t.ratio p.value
15 - 30  0.1979111 0.1034653 62   1.913  0.1436
15 - 61  0.4199143 0.1034653 62   4.059  0.0004
30 - 61  0.2220032 0.1034653 62   2.146  0.0890

P value adjustment: tukey method for comparing a family of 3 estimates

> 2*pt(4.059, 62, lower=FALSE)  # too low
[1] 0.0001405038

> 2*ptukey(4.059, nmeans = 3, df = 62, lower=FALSE)  # too high
[1] 0.03053126


How should I calculate the Tukey-adjusted p-value so that I can obtain more digits?

• If you save the output of res <- summary(pairs(emmeans(model2, "row.space")), you can access the exact $p$-values in this object (e.g. res$p.value). The$p\$-values themselves are already Tukey-adjusted, as the output states. – COOLSerdash Apr 8 '19 at 15:51
• I will maintain that a p value that small is NOT accurate to more than one digit, and probably not that. Such precision would require precise knowledge that the tails of the sampling distribution are exactly as specified in the model, and you simply cannot verify that, except perhaps with simulated data. – Russ Lenth Apr 8 '19 at 23:32

Tukey-adjusted P values are computed using the ptukey() function in R (Studentized range distribution). This is one of the toughest distributions to compute, among those in common use. The help page for ptukey states:

Note

A Legendre 16-point formula is used for the integral of ptukey. The computations are relatively expensive, especially for qtukey which uses a simple secant method for finding the inverse of ptukey. qtukey will be accurate to the 4th decimal place.

This means that the probabilities are computed using numerical integration. I really doubt that a P value way out in the tail can be computed to many digits' accuracy. I suggest that if you have a habit of reporting P values to 2 digits accuracy, it is a habit worth breaking. Nobody needs to distinguish between P = .00022 and P = .00024, even if 2 digits accuracy is believable.

However, as is suggested in a comment, if you do summary(pairs(...)) or as.data.frame(pairs(...)), you get something that inherits from data.frame, and every P value can be extracted to 9 or so digits. But it is false precision.

Also, here is the correct calculation of that P value:

> 1 - ptukey(4.059 * sqrt(2), 3, 62)
[1] 0.0004082944


Actually, the actual result is somewhere in a range depending on the rounding of $$t$$:

> 1 - ptukey(4.0595 * sqrt(2), 3, 62)
[1] 0.0004076112

> 1 - ptukey(4.0585 * sqrt(2), 3, 62)
[1] 0.0004089787


The $$\sqrt{2}$$ factor is needed because the standardization is based on the SE of one mean rather than the difference of two means; and the Studentized range test is a one-sided test based on the maximum minus the minimum.