Can't place correct order of significant values over my bar graph (emmeans_test rearranges factors)

Question

I'm doing a one way ANOVA using rstatix and want to put significant pairwise test values (using emmeans_test) over my plot.

However, the output from emmeans_test reorders my factors. So the significance values are not correctly placed on the right bars.

Here is example data:

library(rstatix)
library(ggpubr)
library(dplyr)

#dataframe
ER <- read.table(header=TRUE, text="
  Sex  Group ER
  M    V     1046
  M    V     1290
  M    Z     1202
  M    Z     1056
  F    V     8000
  F    V     7859
  F    Z     4000
  F    Z     3409
")

ER <- ER %>%
  set_ref_level("Sex", ref = "M") #set males as my reference

ER$Sex shows my factors in the correct order, males first

However, after I do emmeans_test, the output changes. Maybe in alphabetical order?

res.aov <- ER %>% anova_test(ER ~ Sex * Group)
res.aov

pwc <- ER %>% 
  group_by(Sex) %>%
  emmeans_test(ER ~ Group, p.adjust.method = "bonferroni") #output now has females first. 

pwc

I tried to specify the order, in the way I want it to be in the plot, males first, then females, last females + T. Doing this does not change anything.

pwc1 <- pwc[c(2,1),]
pwc1

#plot data 
e <- ggboxplot(ER, x = "Sex", y = "ER", color = "Group",
               palette = "jco")
print(e)

pwc1 <- pwc1 %>% add_xy_position(x = "Sex")

e + stat_pvalue_manual(pwc1) + labs(subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc1))

My plot is in the correct order (Males, Females) but the significance values are switched. Please let me know what I'm missing! I've been looking for a fix online for days and trying different things to no avail.

(Thank you in advance, also I'm an R noob, so I apologize if I made a mistake in the code/post above)

Answer 1

I'm not very familiar with the rstatix and the ggpubr package. However, as you alreday figured out by your self, after doing the emmeans_test the order of your Sex variable gets lost, i.e. it gets converted to a character and hence you get the alphabetical order. A solution to this issue is to convert Sex in your pwc data back to a factor and set the levels as in your ER dataset using factor(Sex, levels = levels(ER$Sex)):

library(rstatix)
library(ggpubr)
library(dplyr)

#dataframe
ER <- read.table(header=TRUE, text="
  Sex  Group ER
  M    V     1046
  M    V     1290
  M    Z     1202
  M    Z     1056
  F    V     8000
  F    V     7859
  F    Z     4000
  F    Z     3409
")

ER <- ER %>%
  mutate(Sex = factor(Sex)) %>% 
  set_ref_level("Sex", ref = "M") #set males as my reference

res.aov <- ER %>% 
  anova_test(ER ~ Sex * Group)
#> Coefficient covariances computed by hccm()

pwc <- ER %>% 
  group_by(Sex) %>%
  emmeans_test(ER ~ Group, p.adjust.method = "bonferroni") %>% #output now has females first.
  # Set levels of Sex as in your ER dataset
  mutate(Sex = factor(Sex, levels = levels(ER$Sex)))

pwc1 <- pwc[c(2,1),]

#plot data 
e <- ggboxplot(ER, x = "Sex", y = "ER", color = "Group",
               palette = "jco")

pwc1 <- pwc1 %>% add_xy_position(x = "Sex")

e + stat_pvalue_manual(pwc1) + 
  labs(subtitle = get_test_label(res.aov, detailed = TRUE),
       caption = get_pwc_label(pwc1))