使用TukeyHSD输出自动在ggplot条形图中添加显著性字母

Question

使用此数据...

hogs.sample<-structure(list(Zone = c("B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D"), Levelname = c("Medium", "High", 
"Low", "Med.High", "Med.High", "Med.High", "Med.High", "Med.High", 
"Med.High", "Medium", "Med.High", "Medium", "Med.High", "High", 
"Medium", "High", "Low", "Med.High", "Low", "High", "Medium", 
"Medium", "Med.High", "Low", "Low", "Med.High", "Low", "Low", 
"High", "High", "Med.High", "High", "Med.High", "Med.High", "Medium", 
"High", "Low", "Low", "Med.High", "Low"), hogs.fit = c(-0.122, 
-0.871, -0.279, -0.446, 0.413, 0.011, 0.157, 0.131, 0.367, -0.23, 
0.007, 0.05, 0.04, -0.184, -0.265, -1.071, -0.223, 0.255, -0.635, 
-1.103, 0.008, -0.04, 0.831, 0.042, -0.005, -0.022, 0.692, 0.402, 
0.615, 0.785, 0.758, 0.738, 0.512, 0.222, -0.424, 0.556, -0.128, 
-0.495, 0.591, 0.923)), row.names = c(NA, -40L), groups = structure(list(
    Zone = c("B", "D"), .rows = structure(list(1:20, 21:40), ptype = integer(0), class = c("vctrs_list_of", 
    "vctrs_vctr", "list"))), row.names = c(NA, -2L), class = c("tbl_df", 
"tbl", "data.frame"), .drop = TRUE), class = c("grouped_df", 
"tbl_df", "tbl", "data.frame"))

我试图根据Tukey的HSD方法在下面的图中添加重要性字母...

library(agricolae)
library(tidyverse)
hogs.plot <- ggplot(hogs.sample, aes(x = Zone, y = exp(hogs.fit), 
                                     fill = factor(Levelname, levels = c("High", "Med.High", "Medium", "Low")))) +  
  stat_summary(fun = mean, geom = "bar", position = position_dodge(0.9), color = "black") +  
  stat_summary(fun.data = mean_se, geom = "errorbar", position = position_dodge(0.9), width = 0.2) + 
  labs(x = "", y = "CPUE (+/-1SE)", legend = NULL) + 
  scale_y_continuous(expand = c(0,0), labels = scales::number_format(accuracy = 0.1)) + 
  scale_fill_manual(values = c("midnightblue", "dodgerblue4", "steelblue3", 'lightskyblue')) + 
  scale_x_discrete(breaks=c("B", "D"), labels=c("Econfina", "Steinhatchee"))+
  scale_color_hue(l=40, c = 100)+
 # coord_cartesian(ylim = c(0, 3.5)) +
  labs(title = "Hogs", x = "", legend = NULL) + 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(),
        axis.text.x = element_text(), axis.title.x = element_text(vjust = 0),
        axis.title.y = element_text(size = 8))+
  theme(legend.title = element_blank(), 
        legend.position = "none")
hogs.plot

我的理想输出应该是这样的...

我不确定这些字母在我的样本图中是否完全准确，但它们表示哪些组之间存在显著差异。区域是独立的，所以我不想在两个区域之间进行比较，因此我使用以下代码分别运行它们。

hogs.aov.b <- aov(hogs.fit ~Levelname, data = filter(hogs.sample, Zone == "B"))
hogs.aov.summary.b <- summary(hogs.aov.b)
hogs.tukey.b <- TukeyHSD(hogs.aov.b)
hogs.tukey.b

hogs.aov.d <- aov(hogs.fit ~ Levelname, data = filter(hogs.sample, Zone == "D"))
hogs.aov.summary.d <- summary(hogs.aov.d)
hogs.tukey.d <- TukeyHSD(hogs.aov.d)
hogs.tukey.d

我尝试了这条路线，但除了猪之外，还有很多其他物种可以应用于此。 显示图表中的显著差异 我可以一次获取一个区域的字母，但我不确定如何将两个区域添加到图中。它们也是无序的。我修改了来自网页的此代码，字母不能很好地放置在条形图上方。

library(agricolae)
library(tidyverse)
# get highest point overall
abs_max <- max(bass.dat.d$bass.fit)
# get the highest point for each class
maxs <- bass.dat.d %>%
  group_by(Levelname) %>%
  # I like to add a little bit to each value so it rests above
  # the highest point. Using a percentage of the highest point
  # overall makes this code a bit more general
  summarise(bass.fit=max(mean(exp(bass.fit))))
# get Tukey HSD results
Tukey_test <- aov(bass.fit ~ Levelname, data=bass.dat.d) %>%
  HSD.test("Levelname", group=TRUE) %>%
  .$groups %>%
  as_tibble(rownames="Levelname") %>%
  rename("Letters_Tukey"="groups") %>%
  select(-bass.fit) %>%
  # and join it to the max values we calculated -- these are
  # your new y-coordinates
  left_join(maxs, by="Levelname")

还有很多类似的例子，比如https://www.staringatr.com/3-the-grammar-of-graphics/bar-plots/3_tukeys/，但它们都是手动添加文本。有一个能够自动将Tukey输出添加到图表中的代码会很好。

谢谢

- Johnny5ish

Answer 1

我不理解你的数据/分析（例如为什么要在hogs.fit上使用exp()以及这些字母应该代表什么），所以我不确定这是否正确，但是因为没有其他人回答，所以这是我最好的猜测：

正确的示例：

## Source: Rosane Rech
## https://statdoe.com/cld-customisation/#adding-the-letters-indicating-significant-differences
## https://www.youtube.com/watch?v=Uyof3S1gx3M

library(tidyverse)
library(ggthemes)
library(multcompView)

# analysis of variance
anova <- aov(weight ~ feed, data = chickwts)

# Tukey's test
tukey <- TukeyHSD(anova)

# compact letter display
cld <- multcompLetters4(anova, tukey)

# table with factors and 3rd quantile
dt <- group_by(chickwts, feed) %>%
  summarise(w=mean(weight), sd = sd(weight)) %>%
  arrange(desc(w))

# extracting the compact letter display and adding to the Tk table
cld <- as.data.frame.list(cld$feed)
dt$cld <- cld$Letters

print(dt)
#> # A tibble: 6 × 4
#>   feed          w    sd cld  
#>   <fct>     <dbl> <dbl> <chr>
#> 1 sunflower  329.  48.8 a    
#> 2 casein     324.  64.4 a    
#> 3 meatmeal   277.  64.9 ab   
#> 4 soybean    246.  54.1 b    
#> 5 linseed    219.  52.2 bc   
#> 6 horsebean  160.  38.6 c

ggplot(dt, aes(feed, w)) + 
  geom_bar(stat = "identity", aes(fill = w), show.legend = FALSE) +
  geom_errorbar(aes(ymin = w-sd, ymax=w+sd), width = 0.2) +
  labs(x = "Feed Type", y = "Average Weight Gain (g)") +
  geom_text(aes(label = cld, y = w + sd), vjust = -0.5) +
  ylim(0,410) +
  theme_few()

根据您的数据，我最好的猜测：

hogs.sample <- structure(list(Zone = c("B", "B", "B", "B", "B", "B", "B", "B", 
                                       "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "D", 
                                       "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
                                       "D", "D", "D", "D", "D", "D"), Levelname = c("Medium", "High", 
                                                                                    "Low", "Med.High", "Med.High", "Med.High", "Med.High", "Med.High", 
                                                                                    "Med.High", "Medium", "Med.High", "Medium", "Med.High", "High", 
                                                                                    "Medium", "High", "Low", "Med.High", "Low", "High", "Medium", 
                                                                                    "Medium", "Med.High", "Low", "Low", "Med.High", "Low", "Low", 
                                                                                    "High", "High", "Med.High", "High", "Med.High", "Med.High", "Medium", 
                                                                                    "High", "Low", "Low", "Med.High", "Low"), hogs.fit = c(-0.122, 
                                                                                                                                           -0.871, -0.279, -0.446, 0.413, 0.011, 0.157, 0.131, 0.367, -0.23, 
                                                                                                                                           0.007, 0.05, 0.04, -0.184, -0.265, -1.071, -0.223, 0.255, -0.635, 
                                                                                                                                           -1.103, 0.008, -0.04, 0.831, 0.042, -0.005, -0.022, 0.692, 0.402, 
                                                                                                                                           0.615, 0.785, 0.758, 0.738, 0.512, 0.222, -0.424, 0.556, -0.128, 
                                                                                                                                           -0.495, 0.591, 0.923)), row.names = c(NA, -40L), groups = structure(list(
                                                                                                                                             Zone = c("B", "D"), .rows = structure(list(1:20, 21:40), ptype = integer(0), class = c("vctrs_list_of", 
                                                                                                                                                                                                                                    "vctrs_vctr", "list"))), row.names = c(NA, -2L), class = c("tbl_df", 
                                                                                                                                                                                                                                                                                               "tbl", "data.frame"), .drop = TRUE), class = c("grouped_df", 
                                                                                                                                                                                                                                                                                                                                              "tbl_df", "tbl", "data.frame"))

# anova
anova <- aov(hogs.fit ~ Levelname * Zone, data = hogs.sample)

# Tukey's test
tukey <- TukeyHSD(anova)

# compact letter display
cld <- multcompLetters4(anova, tukey)

# table with factors and 3rd quantile
dt <- hogs.sample %>% 
  group_by(Zone, Levelname) %>%
  summarise(w=mean(exp(hogs.fit)), sd = sd(exp(hogs.fit)) / sqrt(n())) %>%
  arrange(desc(w)) %>% 
  ungroup() %>% 
  mutate(Levelname = factor(Levelname,
                            levels = c("High",
                                       "Med.High",
                                       "Medium",
                                       "Low"),
                            ordered = TRUE))

# extracting the compact letter display and adding to the Tk table
cld2 <- data.frame(letters = cld$`Levelname:Zone`$Letters)
dt$cld <- cld2$letters

print(dt)
#> # A tibble: 8 × 5
#>   Zone  Levelname     w     sd cld  
#>   <chr> <ord>     <dbl>  <dbl> <chr>
#> 1 D     High      1.97  0.104  a    
#> 2 D     Med.High  1.69  0.206  ab   
#> 3 D     Low       1.36  0.258  abc  
#> 4 B     Med.High  1.14  0.0872 abc  
#> 5 B     Medium    0.875 0.0641 bcd  
#> 6 D     Medium    0.874 0.111  bcd  
#> 7 B     Low       0.696 0.0837 cd   
#> 8 B     High      0.481 0.118  d

ggplot(dt, aes(x = Levelname, y = w)) + 
  geom_bar(stat = "identity", aes(fill = Levelname), show.legend = FALSE) +
  geom_errorbar(aes(ymin = w - sd, ymax = w + sd), width = 0.2) +
  labs(x = "Levelname", y = "Average hogs.fit") +
  geom_text(aes(label = cld, y = w + sd), vjust = -0.5) +
  facet_wrap(~Zone)

^{这是由 reprex 包 (v2.0.1) 在 2021-10-01 创建的。}

- jared_mamrot

---

两个示例都在2023/01/23之后无法工作，可能是由于group_by中的显式更改导致的：Error in n(): ! 必须在dplyr动词内部使用。 运行 rlang::last_error() 以查看出错的位置。 - nouse

Answer 2

我认为我可以在jared_mamrot的回答上进行扩展！

首先，您将找到一个准备好被复制粘贴的reprex。以下是我对此的一些额外评论。

reprex

hogs.sample <- data.frame(
  stringsAsFactors = FALSE,
  Zone = c("B","B","B","B","B","B",
           "B","B","B","B","B","B","B","B","B","B","B","B",
           "B","B","D","D","D","D","D","D","D","D","D",
           "D","D","D","D","D","D","D","D","D","D","D"),
  Levelname = c("Medium","High","Low",
                "Med.High","Med.High","Med.High","Med.High","Med.High",
                "Med.High","Medium","Med.High","Medium","Med.High",
                "High","Medium","High","Low","Med.High","Low","High",
                "Medium","Medium","Med.High","Low","Low","Med.High",
                "Low","Low","High","High","Med.High","High",
                "Med.High","Med.High","Medium","High","Low","Low",
                "Med.High","Low"),
  hogs.fit = c(-0.122,-0.871,-0.279,-0.446,
               0.413,0.011,0.157,0.131,0.367,-0.23,0.007,0.05,
               0.04,-0.184,-0.265,-1.071,-0.223,0.255,-0.635,
               -1.103,0.008,-0.04,0.831,0.042,-0.005,-0.022,0.692,
               0.402,0.615,0.785,0.758,0.738,0.512,0.222,-0.424,
               0.556,-0.128,-0.495,0.591,0.923)
)

library(tidyverse)

# {emmeans}, {multcomp} & {multcompView} ----------------------------------
library(emmeans)
library(multcomp)
library(multcompView)

# set up model
model <- lm(exp(hogs.fit) ~ Levelname * Zone, data = hogs.sample)

# get (adjusted) means
model_means <- emmeans(object = model,
                       specs = ~ Levelname | Zone) 

# add letters to each mean
model_means_cld <- cld(object = model_means,
                       adjust = "Tukey",
                       Letters = letters,
                       alpha = 0.05)
# show output
model_means_cld
#> Zone = B:
#>  Levelname emmean    SE df lower.CL upper.CL .group
#>  High       0.481 0.199 32  -0.0445     1.01  a    
#>  Low        0.696 0.230 32   0.0884     1.30  ab   
#>  Medium     0.875 0.199 32   0.3488     1.40  ab   
#>  Med.High   1.139 0.133 32   0.7889     1.49   b   
#> 
#> Zone = D:
#>  Levelname emmean    SE df lower.CL upper.CL .group
#>  Medium     0.874 0.230 32   0.2673     1.48  a    
#>  Low        1.362 0.151 32   0.9650     1.76  ab   
#>  Med.High   1.688 0.163 32   1.2592     2.12   b   
#>  High       1.969 0.199 32   1.4436     2.50   b   
#> 
#> Results are given on the exp (not the response) scale. 
#> Confidence level used: 0.95 
#> Conf-level adjustment: sidak method for 4 estimates 
#> P value adjustment: tukey method for comparing a family of 4 estimates 

# format output for ggplot
model_means_cld <- model_means_cld %>% 
  as.data.frame() %>% 
  mutate(Zone = case_when(
    Zone == "B" ~ "Econfina",
    Zone == "D" ~ "Steinhatchee"
  ))

# get ggplot
ggplot(data = model_means_cld,
       aes(x = Levelname, y = emmean, fill = Levelname)) +
  facet_grid(cols = vars(Zone)) +
  geom_bar(stat = "identity", color = "black", show.legend = FALSE) +
  geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2) +
  scale_y_continuous(
    name = "CPUE (adj. mean ± 1 std. error)",
    expand = expansion(mult = c(0, 0.1)),
    labels = scales::number_format(accuracy = 0.1)
  ) +
  xlab(NULL) +
  labs(title = "Hogs",
       caption = "Separaetly per Zone, means followed by a common letter are not significantly different according to the Tukey-test") +
  geom_text(aes(label = str_trim(.group), y = emmean + SE), vjust = -0.5) +
  scale_fill_manual(values = c("midnightblue", "dodgerblue4", "steelblue3", 'lightskyblue')) +
  theme_classic() +
  theme(
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.background = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line = element_line(),
    axis.text.x = element_text(),
    axis.title.x = element_text(vjust = 0),
    axis.title.y = element_text(size = 8)
  ) +
  theme(legend.title = element_blank(),
        legend.position = "none")

^{2021年10月18日由reprex包 (v2.0.1)创建}

评论

• 阅读这一章节，了解如何使用紧凑型字母显示法。请注意，该章节特别介绍了为什么将标题放在ggplot下面。
• 此外，我认为jared_mamrot相对于你的要求做出了一个重要的更改。你可以选择比较所有8个均值或分别比较每个区域内的4个Levelname均值。根据您展示的图表，您选择了第二个选项，并通过emmeans()中的 specs = ~ Levelname | Zone来复制它。如果将其更改为specs = ~ Levelname * Zone，则可以选择选项1并找到jared_mamrot找到的相同字母。两种选项都是有效的，但结果不同，您必须选择您想要的结果。编辑：我在最近的回答中对此进行了详细说明。
• 最后，请注意，如果您正在使用函数（如emmeans()）计算和比较这些均值，则无需单独计算因素水平（组合）的算术平均值。此外，在更复杂的情况下（例如缺少数据），您不应显示这些简单的算术平均值及其标准误差，而应直接使用估计的边际均值，又称最小二乘均值、调整后均值或基于模型的均值。然而，在简单情况下，它们与简单算术平均数相同。

- Paul Schmidt