Двухсторонний ANOVA Tukey Test и боксплот в R - PullRequest
Новая
       41

Двухсторонний ANOVA Tukey Test и боксплот в R

0 голосов
/ 10 октября 2019

Я следую этому учебнику и этому вопросу 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]], ...): факторы недопустимы

Я знаю, что я должен поставить числа, не факторы в квантиле, но я неНе знаю, если я исправлю эту ошибку, на моем графике появятся буквы.

1 Ответ

0 голосов
/ 10 октября 2019 
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"))    

data$trat_adj<-paste(data$trat,data$adj,sep = "_")
data<-as.data.frame(data[-c(2,3)])
str(data)
#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)
}

p <- ggplot(data=data, aes(x=data$trat_adj, y=data$dw1, fill=data$trat_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)) 
p

for (facetk in as.character(unique(data$nozzle))) {   
  subdf <- subset(data, 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','data$trat_adj')
  yvalue <- aggregate(subdf$dw1, list(subdf$trat_adj), data=subdf, quantile, probs=.75)  

  final <- data.frame(labels, yvalue[,2])

  names(final)<-c("letters","trat_adj","dw1")
  final$nozzle <-  facetk

  p <- p + geom_text(data = final,  aes(x=trat_adj, y=dw1, fill=trat_adj,label=letters),
                     vjust=-1.5, hjust=-.5)
}
p

enter image description here

...