# Tukey's HSD after permutation tests. Theoretically possible?

I'm interested in making more use of tidymodels' infer package, which lets one perform a variety of statistical tests by a simple algorithm. Here's an example pipeline for an F-test, where we aim to compare yields (as measured by dried weight of plants) obtained under a control and two different treatment conditions.

library(infer)

observed_f <- PlantGrowth |>
specify(weight ~ group) |>
hypothesise(null = "independence") |>
calculate(stat = "F")

null_distribution <- PlantGrowth |>
specify(weight ~ group) |>
hypothesise(null = "independence") |>
generate(reps = 1000, type = "permute") |>
calculate(stat = "F")

null_distribution |>
visualise() +
shade_p_value(obs_stat = observed_f, direction = "greater")


null_distribution |>
get_p_value(obs_stat = observed_f, direction = "greater")

# A tibble: 1 × 1
p_value
<dbl>
1   0.025


To summarize, it creates an experimental F-distribution by calculating 1000 F-statistics by permutation and gives the probability of observing the F-statistic an F-statistic greater than the one coming from our samples.

Unfortunately, the authors of the package do not seem interested in implementing post-hoc tests, the results of which is almost always needed in my field. Of course I can go and do the equivalent of the t-test with infer for each pair, and then adjust the p-values with p.adjust() with a logic like below:

library(tidyverse)

p_vals <- PlantGrowth |>
distinct(group) |>
pull(group) |>
combn(m = 2, simplify = FALSE) |>
map(\(.x) {
dat <- PlantGrowth |>
filter(group %in% .x)
observed_diff <- dat |>
specify(weight ~ group) |>
calculate(stat = "diff in means")
null_distribution <- dat |>
specify(weight ~ group) |>
hypothesise(null = "independence") |>
generate(reps = 1000, type = "permute") |>
calculate(stat = "diff in means")
p_val <- null_distribution |>
get_p_value(obs_stat = observed_diff, direction = "both")
return(p_val)
})

p_vals |>
unlist() |>

p_value p_value p_value
0.798   0.126   0.036


p.adjust() lets me do the following corrections: "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr" but the problem is that Tukey's HSD is the one that is extensively used in my field, a test which the p.adjust() function does not provide.

I've been trying to understand how Tukey's HSD works by looking at the source code of TukeyHSD(). Apparently the p-values are produced by first calculating the differences between the means of each group. Then the differences are divided by the "standard error", which can be formulated as

The absolute value of the resulting "q-statistic" is then used to find the probability of observing the said statistic a greater q-statistic on a studentized range distribution (ptukey()) with parameters nmeans, which is the number of means and df, which is the residual degrees of freedom from the fitted ANOVA model. Overall, it looks like we make use of residual sum of squares, residual degrees of freedom, sample sizes, number of groups and the studentized range distribution to perform Tukey's HSD. Producing F-statistics or differences in means by resampling makes sense and looks quite straightforward, but I can't really link the ideas of resampling and Tukey's HSD together. Should all parameters utilized in Tukey's HSD be produced by resampling, and then what? I can imagine producing 1000 q-statistics by permutation but what about the studentized range distribution? Can I just produce an experimental studentized range distribution by using those 1000 q-statistics and calculate the probability of observing the q-statistic a greater q-statistic for each comparison? What type of an approach could be taken if one wants to perform Tukey's HSD after a permutation test?

• Why? The big advantage of permutation tests is that you don't need to assume a distribution model, which is something Tukey's HSD test does. If this test is "extensively used" in your field, you are doing something unusual already by doing permutation tests. I recommend you use p.adjust. However, the documention helpfully points out that Holm's method is superior to the Bonferroni correction. Aug 10, 2023 at 14:08
• Dec 6, 2023 at 15:08

Tukey's HSD test is a parametric method that uses the distribution of the range of the (standardized) differences in the means.

The permutation approach would be to compute this distribution of the range based on permuted samples.

It is a bit tricky to use the distribution obtained from permutations. It is not a simulation of the distribution given the null hypothesis being true$$^1$$.

1: This is also true for the F-statistic in ANOVA. The obtained p-value in a permutation test does not give the probability that you get the observed value of the statistic or larger given the null hypothesis is true. However, it is a valid p-value in the sense that: given the hypothesis being true, it's distribution is uniform.

• In a way the use of Tukey's HSD test is just like an alternative to ANOVA, but with a different rejection boundary.

Below is an image that shows this for three groups by plotting the rejection boundary as function of two of the three t-statistics (the third is dependent on the other two)

And the approach behind PERMANOVA can be copied to Tukey's HSD.

One issue might be that rejection is not neccesarily because of different means, but it can also because of different variances. Not just the null hypothesis, but the entire model (also thebimplicit assumptions) is being 'tested'. A rejection can mean that the null is wrong, but it can also mean that the assumptions are wrong.

I do not know whether there is some test based on this principle, but I wonder how sensitive it could be to different distributions for the different groups.

• Possibly something like a rank test could be used as alternative instead and for that you can compute/simulate a distribution of the range of the average ranks under the assumption that the null hypothesis is true.