title: How can I add Tukey HSD p-value brackets from a data frame to ggplots from a list of lists?
body: <p>I have a loop that produces a list of ggplots for any given number of analytes (in this example, 3). An ANOVA is performed and p-value generated and mapped for each ggplot.</p>
  <p>However, <strong>I also want to annotate these plots with Tukey HSD p-value brackets</strong>. If possible, I only want to visualize brackets that are below an adjusted p &lt; 0.05.</p>
  <p>I have written two scripts that perform separately these actions; the output for my loop is a list of lists (ggplot elements), while the output from my Tukey HSD script is a dataframe.</p>
  <p>The brackets that I am referring to would look something like in this example ggplot:</p>
  <p><a href="https://i.stack.imgur.com/zySpp.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/zySpp.png" alt="(View example pairwise comparison ggplot brackets here.)" /></a></p>
  <p>The following is my code. First, the data:</p>
  <pre><code>set.seed(10)
  Label&lt;-as.data.frame(c(&quot;Baseline&quot;,&quot;Baseline&quot;,&quot;Baseline&quot;,&quot;Baseline&quot;,&quot;A&quot;,&quot;B&quot;,&quot;B&quot;,&quot;C&quot;,&quot;D&quot;))
  Label&lt;-do.call(&quot;rbind&quot;, replicate(10, Label, simplify = FALSE))
  colnames(Label)&lt;-&quot;Label&quot;
  Analyte1&lt;- rnorm(90, mean = 1, sd = 1)
  Analyte2&lt;- rnorm(90, mean = 1, sd = 1)
  Analyte3&lt;- rnorm(90, mean = .2, sd = 1)
  df&lt;-cbind(Label,Analyte1,Analyte2,Analyte3)
  </code></pre>
  <p>The following is my ANOVA loop:</p>
  <pre><code>library(dplyr)
  library(ggplot2)
  library(ggpubr)
  library(rstatix)
  library(rlist)
  
  # Select numeric columns to obtain df length.
  sample_list &lt;-
    colnames(select_if(df, is.numeric))
  sample_list
  
  # Create a list where the plots will be saved.
  ANOVA.plots &lt;- list()
  
  # The for loop.
  for (i in 1:length(sample_list)) {
    ANOVA.ggplot &lt;-
      ggplot(
        df,
        aes_string(
          x = &quot;Label&quot;,
          y = sample_list[i],
          fill = &quot;Label&quot;,
          title = sample_list[i],
          outlier.shape = NA
        )
      ) +
      border(color = &quot;black&quot;, size = 2.5) +
      geom_jitter(
        aes(fill = Label),
        shape = 21,
        size = 4,
        color = &quot;black&quot;,
        stroke = 1,
        position = position_jitter()
      ) +
      geom_boxplot(alpha = 0.7,
                   linetype = 1,
                   size = 1.1) +
      theme(legend.position = &quot;none&quot;) +
      scale_fill_brewer(palette = &quot;Set2&quot;) +
      stat_boxplot(geom = &quot;errorbar&quot;,
                   size = 1.1,
                   width = 0.4) +
      theme (axis.title.y = element_blank()) +
      theme (axis.title.x = element_blank()) +
      stat_compare_means(
        inherit.aes = TRUE,
        data = df,
        method =  &quot;anova&quot;,
        paired = FALSE ,
        method.args = list(var.equal = TRUE),
        fontface = &quot;bold.italic&quot;,
        size = 5,
        vjust = 0,
        hjust = 0
      )
    ANOVA.plots[[i]] &lt;-  ANOVA.ggplot
  }
  
  # Print each ggplot element from the list of lists.
  ANOVA.plots
  </code></pre>
  <p>Currently, this chunk generates plots that look similar to this one of Analyte 2.</p>
  <p><a href="https://i.stack.imgur.com/Bae2p.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Bae2p.png" alt="(Analyte 2 ANOVA ggplot available here.)" /></a></p>
  <p>Lastly, this chunk is for obtaining the Tukey HSD data frame.</p>
  <pre><code>TUKEY.length&lt;-length(sample_list)
  
  TUKEY.list &lt;-
    vector(mode = &quot;list&quot;, length = TUKEY.length) # Empty list for looping, where Tukey HSD results are stored.
  
  for (i in 1:TUKEY.length + 1) {
    TUKEY.list[[i]] &lt;-
      aov(df[[i]] ~ df[[1]]) %&gt;% tukey_hsd() # Obtain Tukey HSD p-values from comparing groups.
  }
  
  TUKEY.list &lt;-
    TUKEY.list[lengths(TUKEY.list) &gt; 0L] # Remove empty lists, artifacts from piping.
  
  p.value.threshold &lt;- 0.05
  
  TUKEY.df&lt;-list.rbind(TUKEY.list)
  Analytes &lt;-
    colnames(df[, sapply(df, class) %in% c('integer', 'numeric')]) 
  # Stores the analyte names for Tukey. Only pulls numeric colnames.
  
  # The following was done to properly bind the analyte IDs to the Tukey data frame.
  Analytes.df&lt;-as.data.frame(Analytes)
  Analytes.df$Value&lt;-c(1:nrow(Analytes.df))
  Analytes.df&lt;-do.call(&quot;rbind&quot;, replicate(10, Analytes.df, simplify = FALSE))
  
  Analytes.df &lt;-
    Analytes.df[order(Analytes.df$Value), ] 
  # Orders the dataframe so that the rownames can be assigned to the TUKEY HSD.
  Analytes.df&lt;-as.data.frame(Analytes.df[,-2])
  
  toDrop &lt;- c(&quot;^term&quot;, &quot;^null&quot;) # Columns to drop, left over from Tukey loop.
  
  TUKEY.df&lt;-TUKEY.df[,!grepl(paste(toDrop, collapse = &quot;|&quot;),names(TUKEY.df))]
  TUKEY.df&lt;-cbind(Analytes.df, TUKEY.df)
  TUKEY.df&lt;-filter(TUKEY.df, p.adj &lt;= p.value.threshold)
  </code></pre>
  <p>The desired output would still show the ANOVA p-value, and <strong>have brackets for those groups who have an adjusted p-value below 0.05 following Tukey HSD</strong>.</p>
  <p>In this example, <strong>only Analyte 2, groups A and B</strong>, has an adjusted p-value below 0.05 (i.e. p = 0.000754), therefore a bracket should appear above these two plots. However, the code for mapping these Tukey-HSD brackets should be able to work for more variables (i.e. 20+ Analytes) and capture all pairwise comparisons that are significant for each plot, not just the one I have given.</p>