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)  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  0.0001405038 > 2*ptukey(4.059, nmeans = 3, df = 62, lower=FALSE) # too high  0.03053126
How should I calculate the Tukey-adjusted p-value so that I can obtain more digits?