R中的双向ANOVA Tukey测试和箱线图

我下面这个教程和这个问题1，问题2绘制一个双向箱线图与面和字母表示处理间的差异静。有谁知道如何做这种图？

有关我的尝试的一些代码：

    library(multcompView)
    library(ggplot2)
    #sample dataframe
    data <- structure(list(nozzle = c("XR", "XR", "XR", "XR", "XR", "XR", "XR", "XR", 
                                      "XR", "XR", "XR", "XR", "XR", "XR", "XR", "XR", 
                                      "AIXR", "AIXR", "AIXR", "AIXR", "AIXR", "AIXR", 
                                      "AIXR", "AIXR", "AIXR", "AIXR", "AIXR", "AIXR", 
                                      "AIXR", "AIXR", "AIXR", "AIXR"), 
                           trat = c("Cle 12.8", "Cle 12.8", "Cle 12.8", "Cle 12.8", 
                                    "Cle 34", "Cle 34", "Cle 34", "Cle 34", "Cle 12.8", 
                                    "Cle 12.8", "Cle 12.8", "Cle 12.8", "Cle 34", "Cle 34", 
                                    "Cle 34", "Cle 34", "Cle 12.8", "Cle 12.8", "Cle 12.8", 
                                    "Cle 12.8", "Cle 34", "Cle 34", "Cle 34", "Cle 34", 
                                    "Cle 12.8", "Cle 12.8", "Cle 12.8", "Cle 12.8", "Cle 34", 
                                    "Cle 34", "Cle 34", "Cle 34"), 
                           adj = c("Without", "Without", "Without", "Without", "Without", 
                                   "Without", "Without", "Without", "With", "With", "With", 
                                   "With", "With", "With", "With", "With", "Without", "Without", 
                                   "Without", "Without", "Without", "Without", "Without", "Without", 
                                   "With", "With", "With", "With", "With", "With", "With", "With"), 
                           dw1 = c(3.71, 5.87, 6.74, 1.65, 0.27, 0.4, 0.37, 0.34, 0.24, 0.28, 0.32, 
                                   0.38, 0.39, 0.36, 0.32, 0.28, 8.24, 10.18, 11.59, 6.18, 0.2, 0.23, 
                                   0.2, 0.31, 0.28, 0.25, 0.36, 0.27, 0.36, 0.37, 0.34, 0.19)), row.names = c(NA, -32L), class = c("tbl_df", "tbl", "data.frame"))    

    #function 
    generate_label_df <- function(TUKEY, variable){

      # Extract labels and factor levels from Tukey post-hoc 
      Tukey.levels <- variable[,4]
      Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])

      #I need to put the labels in the same order as in the boxplot :
      Tukey.labels$treatment=rownames(Tukey.labels)
      Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
      return(Tukey.labels)
    }

    # What is the effect of the treatment on the value ?
model=lm( data$dw1~ data$trat:data$adj )
ANOVA=aov(model)

# Tukey test to study each pair of treatment :
TUKEY <- TukeyHSD(x=ANOVA, 'data$trat', conf.level=0.95)

# Tuckey test representation :
plot(TUKEY , las=1 , col="brown")

p <- ggplot(data=data, aes(x=trat , y=dw1, fill=adj)) +  
      geom_boxplot(outlier.shape=NA) +
      facet_grid(~nozzle) +
      scale_fill_brewer(palette="Reds") +
      theme_minimal() +
      theme(legend.position="none") +
      theme(axis.text.x=element_text(angle=45, hjust=1)) 

    for (facetk in as.character(unique(data$nozzle))) {   
      subdf <- subset(data, as.array(nozzle==facetk))
      model=lm(dw1 ~ trat:adj, data=subdf)
      ANOVA=aov(model)
      TUKEY <- TukeyHSD(ANOVA)

      labels <- generate_label_df(TUKEY , TUKEY$`trat:adj`)
      names(labels) <- c('Letters','trat')
      yvalue <- aggregate(.~nozzle, data=subdf, quantile, probs=0.75)  
      final <- merge(labels, yvalue)
      final$nozzle <-  facetk

      p <- p + geom_text(data = final,  aes(x=trat, y=value_y, label=Letters), 
                         vjust=-1.5, hjust=-.5)
    }
    p

我收到这个错误

分位数中的错误.default（X [[i]]，...）：不允许使用因子

我知道我必须输入数字，而不是分位数中的因素，但是我不知道是否解决该错误，否则它将出现在图形中。