I'm struggling to add the pvalues to a multiples ggboxplots. I have this type of table:

# A tibble: 6 x 6
   INFy IL.10  IL.8 MCP.1  TNFa group        
  <dbl> <dbl> <dbl> <dbl> <dbl> <fct>        
1  34.9  48.8  9.51  159.  53.9 Healthy_donor
2  NA    31.5  2.76  151.  35.1 Healthy_donor
3  19.4  52.4  3.77  156.  42.1 Healthy_donor
4  11.2  39.5 55.9   203.  60.4 Healthy_donor

levels(df$group)
[1] "Healthy_donor" "Severe"        "Convalescent" 

I pass column names and y-labels separately

cols <- setdiff(names(df), 'group')
y_labels <- c("INFy", 
              "IL.10",
              "IL.8",
              "MCP.1",
              "TNFa")

Then I used a map function and ggarrange to bind them together

Map(function(x, y) {
  ggboxplot(data = df, x = "group", y = x,
            fill = "group",
            palette = "jcb",
            ylab = y, xlab = "Group",
            add = "jitter",
            shape = "group")
}, cols, y_labels) -> list_plots

ggarrange(plotlist = list_plots, common.legend = TRUE)

The result is this:

I like the result, but I can't find how to add the pvalues, I tried with stat_pvalue_manual() a list of pairwise comparisons using pairwise_t_test() from rstatix package. Other approach I tried was reorder my table and use the facet.by parameter of ggboxplot()

    > df %>% sample_n_by(group, size = 2)
    # A tibble: 6 x 3
      X.pg.ml. cytokine group        
         <dbl> <chr>    <fct>        
    1     33.5 IL-10    Convalescent 
    2     10.5 IL-10    Convalescent 
    3     42.1 TNFa     Healthy_donor
    4     NA   IFNy     Healthy_donor
    5     54.8 IL-8     Severe       
    6   1203.  MCP-1    Severe  

stat.test <- df %>%
  group_by(cytokine) %>%
  pairwise_t_test(X.pg.ml. ~ group, 
                  pool.sd = TRUE,
                  p.adjust.method = "bonferroni")
stat.test <- stat.test %>% add_xy_position(x = "group")

ggboxplot(df, x = "group", 
          y = "X.pg.ml.", 
          facet.by = "cytokine",
          fill = "group",
          palette = "jcb",
          add = "jitter",
          shape = "group") +
  stat_pvalue_manual(stat.test, label = "p.adj.signif", tip.length = 0.01)

The result is this:

I don't like it, but I'm not sure if I can change the scale for each plot. I've read the documentation from ggpubr to try to solve this but I'm stuck. Is there a way to place the pvalues within the list of plots and ggarrange() (First figure)??