Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It only takes a minute to sign up.
Sign up to join this community
What test for pairwise comparisons does emmeans uses by default, when executing
emmeans(model, pairwise~predictor)?
As far as I can understand the Tukey method (Tukey HSD) is used by default just for p-values adjustment, not for pairwise comparisons by themselves. My R knowledge is too poor to deconstruct the raw code of emmeans on Github, so hope someone will shed light on the issue.
OK, I think I understand it now. The answer is that the term "adjusted P value" doesn't necessarily mean that a set of P values is obtained, and then we adjust those values somehow. There are some cases where this is true, but others where we compute them differently to begin with.
To illustrate, let's look at one of the built-in datasets in emmeans
> fiber.lm <- lm(strength ~ diameter + machine, data=fiber)
> library(emmeans)
> emm <- emmeans(fiber.lm, "machine")
> pairs(emm, adjust = "none")
contrast estimate SE df t.ratio p.value
A - B -1.04 1.01 11 -1.024 0.3280
A - C 1.58 1.11 11 1.431 0.1803
B - C 2.62 1.15 11 2.283 0.0433
Let's do these manually. First, get the t ratios:
> t.rat <- c(-1.024, 1.431, 2.283)
Calculate the unadjusted P values; these are twice the right-hand tail areas:
> (upv <- 2 * (1 - pt(abs(t.rat), 11)))
[1] 0.32782764 0.18022061 0.04330727
These match the results from pairs().
Now, if we want a Bonferroni adjustment, we adjust these by multiplying by the number of tests:
> 3 * upv
[1] 0.9834829 0.5406618 0.1299218
You can verify this using pairs(emm, adjust = "bonf") (results not shown).
If we want a Tukey adjustment, we use the Studentized Range distribution and compute the right-hand tail areas; we also have to multiply the t ratios by the square root of 2 because of how the distribution is scaled:
> 1 - ptukey(sqrt(2) * abs(t.rat), 3, 11)
[1] 0.5778266 0.3595113 0.1006032
You can verify this with the software:
> pairs(emm, adjust = "tukey")
contrast estimate SE df t.ratio p.value
A - B -1.04 1.01 11 -1.024 0.5781
A - C 1.58 1.11 11 1.431 0.3596
B - C 2.62 1.15 11 2.283 0.1005
P value adjustment: tukey method for comparing a family of 3 estimates
Note that the Bonferroni adjustment is indeed an adjustment to the regular P values upv. But the Tukey-adjusted P values are an entirely separate calculation based on the t ratios but not on upv. So there is no "default" set of P values underlying the Tukey-adjusted P values.
I hope this adds clarity.
