By running our ANOVA, checking the p-value, and finding the critical f-value, we can determine that there is a significant difference in the body mass of penguins across species. The problem is that we don’t know which species differ from each other. To find out, we need to perform a post-hoc test. In this case, we will use Tukey’s HSD Test.
Tukey’s HSD (Honestly Significant Difference) test is a post-hoc test used after ANOVA to determine which specific groups differ from each other.
Why Tukey’s HSD?
Tukey’s HSD test performs pairwise comparisons between all groups. It also controls for Type I error (false positives) when doing multiple comparisons. Additionally, it provides confidence intervals for the differences between means and is considered one of the most powerful post-hoc tests.
Setting Up Our Analysis
Let’s reload the penguins dataset and perform ANOVA on body mass by species:
import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport seaborn as snsimport scipy.stats as stats# Load the penguins dataseturl ="https://raw.githubusercontent.com/dataprofessor/data/refs/heads/master/penguins_cleaned.csv"df = pd.read_csv(url)# Perform ANOVA# adelie = df.loc[df["species"] == "Adelie", "body_mass_g"]# chinstrap = df.loc[df["species"] == "Chinstrap", "body_mass_g"]# gentoo = df.loc[df["species"] == "Gentoo", "body_mass_g"]vars= [ df.loc[df["species"] == species, "body_mass_g"]for species in df["species"].unique()] # this loop is equivalent to the three lines of code above, but it is more efficient and scalable if we have more than three groupsf, p = stats.f_oneway(*vars)alpha =0.05dfb =len(vars) -1# degrees of freedom between groupsdfw =len(df) -len(vars) # degrees of freedom within groupsf_crit = stats.f.ppf(1- alpha, dfb, dfw)print(f"ANOVA Results:")print(f"F-statistic: {f:.4f}")print(f"Critical value: {f_crit:.4f}")print(f"p-value: {p:.6f}")
Since our ANOVA results indicate a significant difference (p < 0.05) and the F-statistic is greater than the critical value, we can proceed with Tukey’s HSD test to identify which specific groups differ from each other.
Performing Tukey’s HSD Test
We will use both scipy.stats and statsmodels to perform Tukey’s HSD test. You can install statsmodels via uv if you haven’t already. See the commented line in the code below. The two libraries provide different ways to perform Tukey’s HSD. The statsmodels library is more user-friendly for this purpose, but it calculates the mean differences slightly differently than scipy.stats and ttest_ind. Whereas a positive mean difference in scipy.stats indicates that the first group has a higher mean than the second group, in statsmodels, a positive mean difference indicates that the second group has a higher mean than the first group. This variation is because statsmodels calculates the mean difference as (group2 - group1), while scipy.stats calculates it as (group1 - group2). Keep this in mind when interpreting the results. We will show you how to use and interpret both methods.
With scipy.stats
Let’s start with scipy.stats, which operates similarly as the scipy.stats.f_oneway function we used for ANOVA. We will manually specify each group for comparison. We have already done this for ANOVA, so we can reuse the vars we created earlier.
from scipy.stats import tukey_hsdscipy_tukey = tukey_hsd(*vars)print(scipy_tukey)print(df["species"].unique()) # to know the order of groups
Now let’s interpret the results from scipy.stats. The output table contains several important columns:
group1, group2: The two groups being compared.
Here we want to reference our df["species"].unique() output to know which species correspond to which group number. Since we looped over this unique list to create vars, the order is preserved.
meandiff: The difference between group means (group1 - group2).
Note that in scipy.stats, a positive mean difference indicates that group1 has a higher mean than group2.
p-adj: The adjusted p-value for multiple comparisons.
lower, upper: The 95% confidence interval for the mean difference.
In order to reject the null hypothesis (indicate that there is a significant difference between the two groups), we need to check two things. First, the p-value (p-adj) must be less than 0.05. Second, both the lower and upper confidence intervals must have the same sign (both positive or both negative). If the confidence interval crosses zero (one bound is negative, and the other is positive), we cannot conclude that there is a significant difference between the two groups. Then the signs indicate which group has the higher mean. If the lower and upper bounds are both positive, group1 has a higher mean than group2. If both bounds are negative, group1 has a lower mean than group2.
Looking at our output, we can conclude that Gentoo (variable index 1) penguins are significantly bigger than both Adelie (variable index 0) and Chinstrap (variable index 2) penguins. The p-adj values for comparisons of the Adelie-Chinstrap comparison are both over 0.05. In addition, the lower and upper limits of the confidence interval for these comparisons cross zero. It is a clear indication that there is no significant difference in body mass between Adelie and Chinstrap penguins.
On the contrary, Gentoo penguins compared to both Adelie and Chinstrap penguins have p-adj values effectively equivalent to zero, and both the lower and upper confidence intervals are the same sign. The signs indicate that Gentoo penguins’ mean difference is greater than both Adelie and Chinstrap penguins.
With statsmodels
If we look at the implementation of Tukey’s HSD test using statsmodels, we can see that it is more easily but reverses the order of group1 and group2.
Tukey’s HSD test is very easy to implement. To perform Tukey’s HSD test, we need to specify the dependent variable (body mass) into the endog parameter and the grouping variable (species) into the groups parameter. We will not have to manually specify each group as we did in ANOVA; instead, we can pass the entire DataFrame columns directly.
In your terminal, use uv to install statsmodels if you haven’t already:
uv add statsmodels
Then you can accomplish the following code.
from statsmodels.stats.multicomp import pairwise_tukeyhsd# Perform Tukey's HSD testtukey_result = pairwise_tukeyhsd( endog=df["body_mass_g"], # The dependent variable (numeric) groups=df["species"], # The grouping variable (categorical) alpha=0.05, # Significance level)# Display the resultstukey_result.summary()
Multiple Comparison of Means - Tukey HSD, FWER=0.05
group1
group2
meandiff
p-adj
lower
upper
reject
Adelie
Chinstrap
26.9239
0.9164
-132.3528
186.2005
False
Adelie
Gentoo
1386.2726
0.0
1252.2897
1520.2554
True
Chinstrap
Gentoo
1359.3487
0.0
1194.4304
1524.2671
True
As you can see, statsmodels makes it very straightforward to perform Tukey’s HSD test. It is also easy to interpret the results. It indicates the name of the variables for group1 and group2, along with the mean difference between the two, the p-value adjusted for multiple comparisons, the confidence interval for the mean difference, and whether to reject the null hypothesis. To read the results, we first focus on the p-value (p-adj). If the p-value is less than 0.05, we know the findings are statistically significant. Then we look at the lower and upper (these are the confidence intervals). They mean that if we conducted this experiment 100 times, 95 times out of 100, the mean difference would fall within this range. If both the lower and the upper bounds have the same sign (both positive or both negative), we can conclude that there is a significant difference between the two groups. If the confidence interval crosses zero (one bound is negative, and the other is positive), we cannot conclude that there is a significant difference between the two groups.
Then the signs indicate which group has the higher mean. This is where it becomes confusing because it doesn’t follow the same convention as scipy.stats. In statsmodels, a positive mean difference indicates that group2 has a higher mean than group1. Therefore:
- Gentoo penguins (group2) have a significantly higher body mass than both Adelie and Chinstrap penguins (group1),
- This is indicated by the positive mean differences and the rejection of the null hypothesis.
- The comparison between Adelie and Chinstrap penguins does not show a significant difference, as indicated by the p-value greater than 0.05 and the confidence interval crossing zero.