Boxplot t-test significance different than ANOVA significance

I have run a one-way ANOVA and turkey post-hoc analysis using the aov() and TukeyHSD() functions. I want to display the p-values (i.e., p adj) associated with the TukeyHSD output on a boxplot. The only way I know how is with the stat_compare_means() function from the ggpubr package.

The issue is that I can only run t-tests to get p-values with the stat_compare_means() function and this returns different p-values compared to the TukeyHSD output. How can I make the boxplots show the p adj values from the TukeyHSD output?

Note: For the example below, I am only interested in comparing the mean price of diamond cuts Fair with all other diamond cuts (e.g., Good, Very Good, Premium, and Ideal).

library(ggplot2)
#> Warning: package 'ggplot2' was built under R version 4.2.3
library(ggpubr)
library(tidyverse)
#> Warning: package 'tidyverse' was built under R version 4.2.3
#> Warning: package 'tibble' was built under R version 4.2.3
#> Warning: package 'tidyr' was built under R version 4.2.3
#> Warning: package 'readr' was built under R version 4.2.3
#> Warning: package 'purrr' was built under R version 4.2.3
#> Warning: package 'dplyr' was built under R version 4.2.3
#> Warning: package 'stringr' was built under R version 4.2.3
#> Warning: package 'forcats' was built under R version 4.2.3
#> Warning: package 'lubridate' was built under R version 4.2.3

# Run the one-way ANOVA and Tukey post-hoc test
# Note the cut comparisons between Fair and all others (first 4 rows)
dat_aov <- aov(price~cut, data = diamonds)
TukeyHSD(dat_aov)
#>   Tukey multiple comparisons of means
#>     95% family-wise confidence level
#> 
#> Fit: aov(formula = price ~ cut, data = diamonds)
#> 
#> $cut
#>                          diff         lwr        upr     p adj
#> Good-Fair          -429.89331  -740.44880  -119.3378 0.0014980
#> Very Good-Fair     -376.99787  -663.86215   -90.1336 0.0031094
#> Premium-Fair        225.49994   -59.26664   510.2665 0.1950425
#> Ideal-Fair         -901.21579 -1180.57139  -621.8602 0.0000000
#> Very Good-Good       52.89544  -130.15186   235.9427 0.9341158
#> Premium-Good        655.39325   475.65120   835.1353 0.0000000
#> Ideal-Good         -471.32248  -642.36268  -300.2823 0.0000000
#> Premium-Very Good   602.49781   467.76249   737.2331 0.0000000
#> Ideal-Very Good    -524.21792  -647.10467  -401.3312 0.0000000
#> Ideal-Premium     -1126.71573 -1244.62267 -1008.8088 0.0000000

# Make a boxplot that displays the p-values of diamond cut
# comparisons with cut Fair
# Note that from TukeyHSD, the Premium-Fair p-val is 0.195 but
# using t.test below it is 0.019
diamonds %>%
  ggplot(aes(x = fct_rev(cut), y = price)) +
  geom_boxplot() +
  coord_flip() +
  xlab('Diamond Cut') +
  stat_compare_means(method = "t.test",
                     comparisons = list(c('Fair','Good'),
                                        c('Fair','Very Good'),
                                        c('Fair','Premium'),
                                        c('Fair','Ideal')),
                     label = "p.format",
                     tip.length = 0) +
  theme_bw()

^{Created on 2023-06-25 with reprex v2.0.2}

You can use base boxplots.

data('diamonds', package='ggplot2')
tuk <- TukeyHSD(aov(price ~ cut, data=diamonds))
par(mar=c(4, 7, 4, 2)+.1)
b <- boxplot(price ~ cut, data=diamonds, horizontal=TRUE, col=0, pch=20, cex=.8,
las=1, ylab='', ylim=c(0, max(diamonds$price)*1.5), border='grey15')
mtext('Diamonds cut', 2, 6)
mx <- max(diamonds$price)
for (i in sq <- seq_along(b$names)[-length(b$names)]) segments(mx + 2e3*i, 5, mx + 2e3*i, 5 - i)
text(mx + 2e3*sq + 750, seq.int(4.5, by=-.5, length.out=4), 
signif(tuk$cut[1:4, 4], 3), adj=.5, srt=270, cex=.9, col=c('red', 'blue')[(tuk$cut < .05) + 1])

Here is a workaround using ggplot2 that takes advantage of the annotate() function.

library(ggplot2)
#> Warning: package 'ggplot2' was built under R version 4.2.3
library(ggpubr)
library(tidyverse)
#> Warning: package 'tidyverse' was built under R version 4.2.3
#> Warning: package 'tibble' was built under R version 4.2.3
#> Warning: package 'tidyr' was built under R version 4.2.3
#> Warning: package 'readr' was built under R version 4.2.3
#> Warning: package 'purrr' was built under R version 4.2.3
#> Warning: package 'dplyr' was built under R version 4.2.3
#> Warning: package 'stringr' was built under R version 4.2.3
#> Warning: package 'forcats' was built under R version 4.2.3
#> Warning: package 'lubridate' was built under R version 4.2.3

# Run the one-way ANOVA and Tukey post-hoc test
# Note the cut comparisons between Fair and all others (first 4 rows)
dat_aov <- aov(price~cut, data = diamonds)
TukeyHSD(dat_aov)
#>   Tukey multiple comparisons of means
#>     95% family-wise confidence level
#> 
#> Fit: aov(formula = price ~ cut, data = diamonds)
#> 
#> $cut
#>                          diff         lwr        upr     p adj
#> Good-Fair          -429.89331  -740.44880  -119.3378 0.0014980
#> Very Good-Fair     -376.99787  -663.86215   -90.1336 0.0031094
#> Premium-Fair        225.49994   -59.26664   510.2665 0.1950425
#> Ideal-Fair         -901.21579 -1180.57139  -621.8602 0.0000000
#> Very Good-Good       52.89544  -130.15186   235.9427 0.9341158
#> Premium-Good        655.39325   475.65120   835.1353 0.0000000
#> Ideal-Good         -471.32248  -642.36268  -300.2823 0.0000000
#> Premium-Very Good   602.49781   467.76249   737.2331 0.0000000
#> Ideal-Very Good    -524.21792  -647.10467  -401.3312 0.0000000
#> Ideal-Premium     -1126.71573 -1244.62267 -1008.8088 0.0000000

# Make a boxplot that displays the p-values of diamond cut
# comparisons with cut Fair
diamonds %>%
  ggplot(aes(x = fct_rev(cut), y = price)) +
  geom_boxplot() +
  coord_flip(ylim = c(0,30000)) +
  xlab('Diamond Cut') +
  annotate('segment', x = c(5,5,5,5), xend = c(4,3,2,1), y = c(20000,23000,26000,29000), yend = c(20000,23000,26000,29000), color = 'black') +
  annotate('text', x = c(4.5,4,3.5,3), y = c(20900,23900,26900,29900),label = c('0.001','0.003','> 0.05', '< 0.001'), color = 'black', angle = -90) +
  theme_bw()

^{Created on 2023-06-26 with reprex v2.0.2}