近日，在学习中遇到了一个问题，问题如下：
我想要绘制这样的图形，按分组在组内标注出差异显著性，要求字母标注，如下图。
碰到的问题是，箱线图展示出来的数据并非是平均值的值，而显著性字母标记是在分组统计到最终的平均值后得到的，与平均值一一对应，所以这里画完箱线图，想要直接将字母标注映射到图中是不可能，会报错，报错提示如下：
Error in `geom_text()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 2nd layer.
Caused by error in `check_aesthetics()`:
! Aesthetics must be either length 1 or the same as the data (18)
✖ Fix the following mappings: `label`
Run `rlang::last_trace()` to see where the error occurred.
虽然看不大懂，但是大概知道是映射出的问题，我的理解就是字母标注和数据未能做到一一对应。
摸索了好一阵子，通过网络搜索学习，最终解决方案如下：
首先，我构建一个示例数据，纯属乱编，仅为这个案例使用。
# 构建数据
group = c(rep("a", 9), rep("b", 9))
Treatment = rep(c("T1","T2","T3"), each = 3, 2)
Block = rep(c("B1","B2","B3"), 6)
V = c(1.5,1.8,1.1,2.6,3.2,3.8,3.7,4.1,4.9,
2.7,3.4,2.1,4.1,4.9,5.2,5.8,6.2,6.5)sample.boxplot <- data.frame(group, Treatment, Block, V)
sample.boxplot
## group Treatment Block V
## 1 a T1 B1 1.5
## 2 a T1 B2 1.8
## 3 a T1 B3 1.1
## 4 a T2 B1 2.6
## 5 a T2 B2 3.2
## 6 a T2 B3 3.8
## 7 a T3 B1 3.7
## 8 a T3 B2 4.1
## 9 a T3 B3 4.9
## 10 b T1 B1 2.7
## 11 b T1 B2 3.4
## 12 b T1 B3 2.1
## 13 b T2 B1 4.1
## 14 b T2 B2 4.9
## 15 b T2 B3 5.2
## 16 b T3 B1 5.8
## 17 b T3 B2 6.2## 18 b T3 B3 6.5library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors先进行方差分析。
sample.boxplot.aov <- aov(V ~ Block + group * Treatment, data = sample.boxplot)
anova(sample.boxplot.aov)## Analysis of Variance Table
##
## Response: V
## Df Sum Sq Mean Sq F value Pr(>F)
## Block 2 1.1378 0.5689 2.4568 0.1355
## group 1 11.2022 11.2022 48.3781 3.921e-05 ***
## Treatment 2 29.2311 14.6156 63.1190 2.131e-06 ***
## group:Treatment 2 0.3378 0.1689 0.7294 0.5062
## Residuals 10 2.3156 0.2316
## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1从结果可以看出主效应group和Treatment分别对V影响显著，无显著交互作用。进一步进行主效应的多重比较。
library(agricolae)
# 主效应group
sampleout1 <- LSD.test(sample.boxplot.aov,"group",p.adj="none")
sampleout1## $statistics
## MSerror Df Mean CV t.value LSD
## 0.2315556 10 3.755556 12.81308 2.228139 0.505433
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD none group 2 0.05
##
## $means
## V std r se LCL UCL Min Max Q25 Q50 Q75
## a 2.966667 1.296148 9 0.1604007 2.609272 3.324062 1.1 4.9 1.8 3.2 3.8
## b 4.544444 1.564538 9 0.1604007 4.187049 4.901840 2.1 6.5 3.4 4.9 5.8
##
## $comparison
## NULL
##
## $groups
## V groups
## b 4.544444 a
## a 2.966667 b
##
## attr(,"class")
## [1] "group"
# 主效应Treatment
sampleout2 <- LSD.test(sample.boxplot.aov,"Treatment",p.adj="none")
sampleout2## $statistics
## MSerror Df Mean CV t.value LSD
## 0.2315556 10 3.755556 12.81308 2.228139 0.6190265
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD none Treatment 3 0.05
##
## $means
## V std r se LCL UCL Min Max Q25 Q50 Q75
## T1 2.100000 0.8366600 6 0.19645 1.662282 2.537718 1.1 3.4 1.575 1.95 2.55
## T2 3.966667 0.9892758 6 0.19645 3.528949 4.404384 2.6 5.2 3.350 3.95 4.70
## T3 5.200000 1.1489125 6 0.19645 4.762282 5.637718 3.7 6.5 4.300 5.35 6.10
##
## $comparison
## NULL
##
## $groups
## V groups
## T3 5.200000 a
## T2 3.966667 b
## T1 2.100000 c
##
## attr(,"class")## [1] "group"因为最终成图的字母标注是在各组内标记的，所以这里我对各组分别进行单因素方差分析，并多重比较。
# a组组内处理差异
sample.boxplota <- sample.boxplot %>% filter(group == "a")
sample.boxplota## group Treatment Block V
## 1 a T1 B1 1.5
## 2 a T1 B2 1.8
## 3 a T1 B3 1.1
## 4 a T2 B1 2.6
## 5 a T2 B2 3.2
## 6 a T2 B3 3.8
## 7 a T3 B1 3.7
## 8 a T3 B2 4.1## 9 a T3 B3 4.9sample.boxplota.aov <- aov(V ~ Block + Treatment, data = sample.boxplota)
anova(sample.boxplota.aov)## Analysis of Variance Table
##
## Response: V
## Df Sum Sq Mean Sq F value Pr(>F)
## Block 2 0.6867 0.3433 1.3377 0.359067
## Treatment 2 11.7267 5.8633 22.8442 0.006481 **
## Residuals 4 1.0267 0.2567
## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sampleouta <- LSD.test(sample.boxplota.aov, "Treatment", p.adj="none")
sampleouta## $statistics
## MSerror Df Mean CV t.value LSD
## 0.2566667 4 2.966667 17.07717 2.776445 1.148493
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD none Treatment 3 0.05
##
## $means
## V std r se LCL UCL Min Max Q25 Q50 Q75
## T1 1.466667 0.3511885 3 0.2924988 0.6545598 2.278774 1.1 1.8 1.3 1.5 1.65
## T2 3.200000 0.6000000 3 0.2924988 2.3878931 4.012107 2.6 3.8 2.9 3.2 3.50
## T3 4.233333 0.6110101 3 0.2924988 3.4212264 5.045440 3.7 4.9 3.9 4.1 4.50
##
## $comparison
## NULL
##
## $groups
## V groups
## T3 4.233333 a
## T2 3.200000 a
## T1 1.466667 b
##
## attr(,"class")## [1] "group"# b组组内处理差异
sample.boxplotb <- sample.boxplot %>% filter(group == "b")
sample.boxplotb## group Treatment Block V
## 1 b T1 B1 2.7
## 2 b T1 B2 3.4
## 3 b T1 B3 2.1
## 4 b T2 B1 4.1
## 5 b T2 B2 4.9
## 6 b T2 B3 5.2
## 7 b T3 B1 5.8
## 8 b T3 B2 6.2## 9 b T3 B3 6.5sample.boxplotb.aov <- aov(V ~ Block + Treatment, data = sample.boxplotb)
anova(sample.boxplotb.aov)## Analysis of Variance Table
##
## Response: V
## Df Sum Sq Mean Sq F value Pr(>F)
## Block 2 0.6156 0.3078 1.0949 0.417616
## Treatment 2 17.8422 8.9211 31.7352 0.003515 **
## Residuals 4 1.1244 0.2811
## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sampleoutb <- LSD.test(sample.boxplotb.aov, "Treatment", p.adj="none")
sampleoutb## $statistics
## MSerror Df Mean CV t.value LSD
## 0.2811111 4 4.544444 11.66697 2.776445 1.201939
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD none Treatment 3 0.05
##
## $means
## V std r se LCL UCL Min Max Q25 Q50 Q75
## T1 2.733333 0.6506407 3 0.3061106 1.883434 3.583233 2.1 3.4 2.4 2.7 3.05
## T2 4.733333 0.5686241 3 0.3061106 3.883434 5.583233 4.1 5.2 4.5 4.9 5.05
## T3 6.166667 0.3511885 3 0.3061106 5.316767 7.016566 5.8 6.5 6.0 6.2 6.35
##
## $comparison
## NULL
##
## $groups
## V groups
## T3 6.166667 a
## T2 4.733333 b
## T1 2.733333 c
##
## attr(,"class")## [1] "group"进一步的对原始数据进行分组统计构建新的数据集sample.boxplot1。
# 分组统计数据，新建数据集sample.boxplot1
sample.boxplot1 <- sample.boxplot %>% group_by(group, Treatment) %>%
summarise(vmean = mean(V), vsd = sd(V))## `summarise()` has grouped output by 'group'. You can override using the## `.groups` argument.sample.boxplot1## # A tibble: 6 × 4
## # Groups: group [2]
## group Treatment vmean vsd
## <chr> <chr> <dbl> <dbl>
## 1 a T1 1.47 0.351
## 2 a T2 3.2 0.6
## 3 a T3 4.23 0.611
## 4 b T1 2.73 0.651
## 5 b T2 4.73 0.569## 6 b T3 6.17 0.351# 查看之前组内单因素方差分析结果的字母标注
sampleouta$groups## V groups
## T3 4.233333 a
## T2 3.200000 a## T1 1.466667 bsampleoutb$groups## V groups
## T3 6.166667 a
## T2 4.733333 b## T1 2.733333 c# 将字母标注添加到sample.boxplot1
sample.boxplot1 <- data.frame(sample.boxplot1,
vsig = c("b", "a", "a", "c", "b", "a"))
sample.boxplot1## group Treatment vmean vsd vsig
## 1 a T1 1.466667 0.3511885 b
## 2 a T2 3.200000 0.6000000 a
## 3 a T3 4.233333 0.6110101 a
## 4 b T1 2.733333 0.6506407 c
## 5 b T2 4.733333 0.5686241 b## 6 b T3 6.166667 0.3511885 a绘图
# 绘图并添加字母标注
sample.boxplot %>%
ggplot(aes(x = Treatment, y = V, fill = group)) +
geom_boxplot() +
geom_point(data = sample.boxplot1,
aes(x = Treatment, y = vmean),
position = position_dodge(0.8), color = "white") +
geom_text(data = sample.boxplot1, aes(x = Treatment, y = vmean + 0.8,
label = vsig),position = position_dodge(0.8))
调整主题美化图形。
# 主题调整美化
sample.boxplot %>%
ggplot(aes(x = Treatment, y = V, fill = group)) +
geom_boxplot() +
geom_point(data = sample.boxplot1,
aes(x = Treatment, y = vmean),
position = position_dodge(0.8), color = "white") +
geom_text(data = sample.boxplot1, aes(x = Treatment, y = vmean + 0.8, label = vsig),
position = position_dodge(0.8)) +
theme_bw() +
theme(text = element_text(size = 14, family = "serif"),
axis.text = element_text(family = "serif", colour = "black")) +
theme(legend.position = "top", legend.title = element_blank(),
text = element_text(size = 14, family = "serif"),
axis.text = element_text(size = 14, family = "serif", colour = "black")) 
本文使用 文章同步助手 同步
