使用 TukeyHSD 的输出自动将重要字母添加到 ggplot 条形图中

Question

使用 TukeyHSD 的输出自动将重要字母添加到 ggplot 条形图中

Joh*_*ish 8 text label r significance ggplot2

使用这些数据...

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"))

Run Code Online (Sandbox Code Playgroud)

我正在尝试将基于 Tukeys 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

Run Code Online (Sandbox Code Playgroud)

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

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

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

Run Code Online (Sandbox Code Playgroud)

我尝试了这条路线，但除了猪之外，我还有很多物种可以应用此方法。在图表中显示统计上显着的差异

我可以一次获取一个区域的字母，但我不确定如何将两个区域添加到图中。它们也出故障了。我从网页修改了这段代码，但这些字母没有很好地放置在栏的顶部。

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")

Run Code Online (Sandbox Code Playgroud)

也有很多这样的例子https://www.staringatr.com/3-the-grammar-of-graphics/bar-plots/3_tukeys/但它们都只是手动添加文本。如果代码能够获取 Tukey 输出并自动将重要性字母添加到图中，那就太好了。

谢谢

Answer 1

jar*_*rot 10

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

\n

正确示例：

\n

## Source: Rosane Rech\n## https://statdoe.com/cld-customisation/#adding-the-letters-indicating-significant-differences\n## https://www.youtube.com/watch?v=Uyof3S1gx3M\n\nlibrary(tidyverse)\nlibrary(ggthemes)\nlibrary(multcompView)\n\n# analysis of variance\nanova <- aov(weight ~ feed, data = chickwts)\n\n# Tukey\'s test\ntukey <- TukeyHSD(anova)\n\n# compact letter display\ncld <- multcompLetters4(anova, tukey)\n\n# table with factors and 3rd quantile\ndt <- group_by(chickwts, feed) %>%\n  summarise(w=mean(weight), sd = sd(weight)) %>%\n  arrange(desc(w))\n\n# extracting the compact letter display and adding to the Tk table\ncld <- as.data.frame.list(cld$feed)\ndt$cld <- cld$Letters\n\nprint(dt)\n#> # A tibble: 6 \xc3\x97 4\n#>   feed          w    sd cld  \n#>   <fct>     <dbl> <dbl> <chr>\n#> 1 sunflower  329.  48.8 a    \n#> 2 casein     324.  64.4 a    \n#> 3 meatmeal   277.  64.9 ab   \n#> 4 soybean    246.  54.1 b    \n#> 5 linseed    219.  52.2 bc   \n#> 6 horsebean  160.  38.6 c\n\nggplot(dt, aes(feed, w)) + \n  geom_bar(stat = "identity", aes(fill = w), show.legend = FALSE) +\n  geom_errorbar(aes(ymin = w-sd, ymax=w+sd), width = 0.2) +\n  labs(x = "Feed Type", y = "Average Weight Gain (g)") +\n  geom_text(aes(label = cld, y = w + sd), vjust = -0.5) +\n  ylim(0,410) +\n  theme_few()\n

Run Code Online (Sandbox Code Playgroud)\n

\n

我对你的数据的“最佳猜测”：

\n

hogs.sample <- structure(list(Zone = c("B", "B", "B", "B", "B", "B", "B", "B", \n                                       "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "D", \n                                       "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", \n                                       "D", "D", "D", "D", "D", "D"), Levelname = c("Medium", "High", \n                                                                                    "Low", "Med.High", "Med.High", "Med.High", "Med.High", "Med.High", \n                                                                                    "Med.High", "Medium", "Med.High", "Medium", "Med.High", "High", \n                                                                                    "Medium", "High", "Low", "Med.High", "Low", "High", "Medium", \n                                                                                    "Medium", "Med.High", "Low", "Low", "Med.High", "Low", "Low", \n                                                                                    "High", "High", "Med.High", "High", "Med.High", "Med.High", "Medium", \n                                                                                    "High", "Low", "Low", "Med.High", "Low"), hogs.fit = c(-0.122, \n                                                                                                                                           -0.871, -0.279, -0.446, 0.413, 0.011, 0.157, 0.131, 0.367, -0.23, \n                                                                                                                                           0.007, 0.05, 0.04, -0.184, -0.265, -1.071, -0.223, 0.255, -0.635, \n                                                                                                                                           -1.103, 0.008, -0.04, 0.831, 0.042, -0.005, -0.022, 0.692, 0.402, \n                                                                                                                                           0.615, 0.785, 0.758, 0.738, 0.512, 0.222, -0.424, 0.556, -0.128, \n                                                                                                                                           -0.495, 0.591, 0.923)), row.names = c(NA, -40L), groups = structure(list(\n                                                                                                                                             Zone = c("B", "D"), .rows = structure(list(1:20, 21:40), ptype = integer(0), class = c("vctrs_list_of", \n                                                                                                                                                                                                                                    "vctrs_vctr", "list"))), row.names = c(NA, -2L), class = c("tbl_df", \n                                                                                                                                                                                                                                                                                               "tbl", "data.frame"), .drop = TRUE), class = c("grouped_df", \n                                                                                                                                                                                                                                                                                                                                              "tbl_df", "tbl", "data.frame"))\n\n# anova\nanova <- aov(hogs.fit ~ Levelname * Zone, data = hogs.sample)\n\n# Tukey\'s test\ntukey <- TukeyHSD(anova)\n\n# compact letter display\ncld <- multcompLetters4(anova, tukey)\n\n# table with factors and 3rd quantile\ndt <- hogs.sample %>% \n  group_by(Zone, Levelname) %>%\n  summarise(w=mean(exp(hogs.fit)), sd = sd(exp(hogs.fit)) / sqrt(n())) %>%\n  arrange(desc(w)) %>% \n  ungroup() %>% \n  mutate(Levelname = factor(Levelname,\n                            levels = c("High",\n                                       "Med.High",\n                                       "Medium",\n                                       "Low"),\n                            ordered = TRUE))\n\n# extracting the compact letter display and adding to the Tk table\ncld2 <- data.frame(letters = cld$`Levelname:Zone`$Letters)\ndt$cld <- cld2$letters\n\nprint(dt)\n#> # A tibble: 8 \xc3\x97 5\n#>   Zone  Levelname     w     sd cld  \n#>   <chr> <ord>     <dbl>  <dbl> <chr>\n#> 1 D     High      1.97  0.104  a    \n#> 2 D     Med.High  1.69  0.206  ab   \n#> 3 D     Low       1.36  0.258  abc  \n#> 4 B     Med.High  1.14  0.0872 abc  \n#> 5 B     Medium    0.875 0.0641 bcd  \n#> 6 D     Medium    0.874 0.111  bcd  \n#> 7 B     Low       0.696 0.0837 cd   \n#> 8 B     High      0.481 0.118  d\n\nggplot(dt, aes(x = Levelname, y = w)) + \n  geom_bar(stat = "identity", aes(fill = Levelname), show.legend = FALSE) +\n  geom_errorbar(aes(ymin = w - sd, ymax = w + sd), width = 0.2) +\n  labs(x = "Levelname", y = "Average hogs.fit") +\n  geom_text(aes(label = cld, y = w + sd), vjust = -0.5) +\n  facet_wrap(~Zone)\n

Run Code Online (Sandbox Code Playgroud)\n

\n

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

\n

Answer 2

Pau*_*idt 7

我想我可以扩展 jared_mamrot 的答案！

\n

首先，您会发现一个可以复制粘贴的表示。下面，我对此有一些补充意见。

\n

代表

\n

hogs.sample <- data.frame(\n  stringsAsFactors = FALSE,\n  Zone = c("B","B","B","B","B","B",\n           "B","B","B","B","B","B","B","B","B","B","B","B",\n           "B","B","D","D","D","D","D","D","D","D","D",\n           "D","D","D","D","D","D","D","D","D","D","D"),\n  Levelname = c("Medium","High","Low",\n                "Med.High","Med.High","Med.High","Med.High","Med.High",\n                "Med.High","Medium","Med.High","Medium","Med.High",\n                "High","Medium","High","Low","Med.High","Low","High",\n                "Medium","Medium","Med.High","Low","Low","Med.High",\n                "Low","Low","High","High","Med.High","High",\n                "Med.High","Med.High","Medium","High","Low","Low",\n                "Med.High","Low"),\n  hogs.fit = c(-0.122,-0.871,-0.279,-0.446,\n               0.413,0.011,0.157,0.131,0.367,-0.23,0.007,0.05,\n               0.04,-0.184,-0.265,-1.071,-0.223,0.255,-0.635,\n               -1.103,0.008,-0.04,0.831,0.042,-0.005,-0.022,0.692,\n               0.402,0.615,0.785,0.758,0.738,0.512,0.222,-0.424,\n               0.556,-0.128,-0.495,0.591,0.923)\n)\n\nlibrary(tidyverse)\n\n# {emmeans}, {multcomp} & {multcompView} ----------------------------------\nlibrary(emmeans)\nlibrary(multcomp)\nlibrary(multcompView)\n\n# set up model\nmodel <- lm(exp(hogs.fit) ~ Levelname * Zone, data = hogs.sample)\n\n# get (adjusted) means\nmodel_means <- emmeans(object = model,\n                       specs = ~ Levelname | Zone) \n\n# add letters to each mean\nmodel_means_cld <- cld(object = model_means,\n                       adjust = "Tukey",\n                       Letters = letters,\n                       alpha = 0.05)\n# show output\nmodel_means_cld\n#> Zone = B:\n#>  Levelname emmean    SE df lower.CL upper.CL .group\n#>  High       0.481 0.199 32  -0.0445     1.01  a    \n#>  Low        0.696 0.230 32   0.0884     1.30  ab   \n#>  Medium     0.875 0.199 32   0.3488     1.40  ab   \n#>  Med.High   1.139 0.133 32   0.7889     1.49   b   \n#> \n#> Zone = D:\n#>  Levelname emmean    SE df lower.CL upper.CL .group\n#>  Medium     0.874 0.230 32   0.2673     1.48  a    \n#>  Low        1.362 0.151 32   0.9650     1.76  ab   \n#>  Med.High   1.688 0.163 32   1.2592     2.12   b   \n#>  High       1.969 0.199 32   1.4436     2.50   b   \n#> \n#> Results are given on the exp (not the response) scale. \n#> Confidence level used: 0.95 \n#> Conf-level adjustment: sidak method for 4 estimates \n#> P value adjustment: tukey method for comparing a family of 4 estimates \n\n# format output for ggplot\nmodel_means_cld <- model_means_cld %>% \n  as.data.frame() %>% \n  mutate(Zone = case_when(\n    Zone == "B" ~ "Econfina",\n    Zone == "D" ~ "Steinhatchee"\n  ))\n\n# get ggplot\nggplot(data = model_means_cld,\n       aes(x = Levelname, y = emmean, fill = Levelname)) +\n  facet_grid(cols = vars(Zone)) +\n  geom_bar(stat = "identity", color = "black", show.legend = FALSE) +\n  geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2) +\n  scale_y_continuous(\n    name = "CPUE (adj. mean \xc2\xb1 1 std. error)",\n    expand = expansion(mult = c(0, 0.1)),\n    labels = scales::number_format(accuracy = 0.1)\n  ) +\n  xlab(NULL) +\n  labs(title = "Hogs",\n       caption = "Separaetly per Zone, means followed by a common letter are not significantly different according to the Tukey-test") +\n  geom_text(aes(label = str_trim(.group), y = emmean + SE), vjust = -0.5) +\n  scale_fill_manual(values = c("midnightblue", "dodgerblue4", "steelblue3", \'lightskyblue\')) +\n  theme_classic() +\n  theme(\n    panel.border = element_blank(),\n    panel.grid.major = element_blank(),\n    panel.background = element_blank(),\n    panel.grid.minor = element_blank(),\n    axis.line = element_line(),\n    axis.text.x = element_text(),\n    axis.title.x = element_text(vjust = 0),\n    axis.title.y = element_text(size = 8)\n  ) +\n  theme(legend.title = element_blank(),\n        legend.position = "none")\n

Run Code Online (Sandbox Code Playgroud)\n

\n

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

\n

您可能会对本章有关使用紧凑字母显示器感兴趣。请注意，这尤其涉及为什么我将标题放在 ggplot 下方。
此外，我相信与您所要求的相比，jared_mamrot 改变了一件重要的事情。您可以选择对每个区域的所有 8 个平均值进行比较，或者对所有 4 个级别名称平均值进行比较。从您展示的情节来看，您选择了第二个选项，我通过specs = ~ Levelname | Zonein在这里复制了它emmeans()。您可以选择选项 1，并通过将其更改为来找到与 jared_mamrot 相同的字母specs = ~ Levelname * Zone。两个选项都有效，但结果不同，您必须选择您想要的。编辑：我在这个更新的答案中对此进行了相当详细的介绍。
emmeans()最后请注意，如果您使用函数（如）来计算和比较这些平均值，则无需单独计算每个因子级别（组合）的算术平均值。此外，在更复杂的情况下，例如丢失数据，您不应该显示那些简单的算术平均值及其标准误差，而应该直接获取估计的边际平均值，又称为最小二乘平均值，又称为调整平均值，又称为基于模型的平均值。然而，在简单的情况下，它们与简单的算术平均值相同。

\n

归档时间：	3 年，11 月前
查看次数：	21823 次
最近记录：	1 年，8 月前

使用 TukeyHSD 的输出自动将重要字母添加到 ggplot 条形图中

正确示例：

我对你的数据的“最佳猜测”：

代表

评论