How to run and annotate separate TukeyHSD for individual facets in ggplot2 boxplot?

Question

I am not very good at R and am trying to pull together this code that is not quite working out how I would like it to. I would really appreciate any help on this!

I would like to perform TukeyHSD test among treatment groups in individual facets in my ggplot boxplots. Currently though, my figures apply a single TukeyHSD across all the boxplots in the figure and this results in a huge number of groupings as you can see in the figure:

my current plot

As I mentioned, it would be preferable to have TukeyHSD run on the individual Depth separated "0" facet, then "5" facet, then "30" facet separately. Is this possible by modifying the code I have been using?

data1 <- read.delim(file="clipboard")

data1$Treatment <- as.factor(data1$Treatment)
data1$Depth <- as.factor(data1$Depth)

model<- aov(MBC~Treatment*Depth, data=data1)
model
library(emmeans)
library('multcomp')
cld_dat = as.data.frame( cld(emmeans(model,~Depth*Treatment),
                               Letters = letters ) )

ggplot(data1, aes(x=Treatment, y=MBC, fill=Treatment)) +
  geom_boxplot() + 
  ylab("MBC") +
  ggtitle("Melinis") +
  facet_wrap(~Depth,ncol=3) +
  geom_text(data = cld_dat, aes(y = 140, label = .group))

One more question, if this is possible: how would I add another y variable "CB" as a second row identical to how I have the first row variable "MBC"?

Thank you for any suggestions!

        Treatment Depth    MBC        CB

Answer 1

If I understand correctly which factor has those facet levels, what you need is

 cld(emmeans(model, ~ Treatment | Depth))