欢迎关注R语言数据分析指南
本节来介绍如何使用R做三因素方差分析，绘制并排堆砌条形图并添加显著性字母标记，下面来看具体例子
清理变量
rm(list=ls())
安装并加载R包
package.list=c("tidyverse","ggthemes","multcompView","egg")
for (package in package.list) {
if (!require(package,character.only=T, quietly=T)) {
install.packages(package)
library(package, character.only=T)
}
}
加载数据
df <- read_tsv("co2.xls")
数据结构
Plant Type Treatment conc uptake
1 Qn1 Quebec nonchilled 95 16.0
2 Qn1 Quebec nonchilled 175 30.4
3 Qn1 Quebec nonchilled 250 34.8
4 Qn1 Quebec nonchilled 350 37.2
5 Qn1 Quebec nonchilled 500 35.3
6 Qn1 Quebec nonchilled 675 39.2
三因素方差分析
anova <- aov(uptake ~ factor(conc)*Type*Treatment,data=df)
summary(anova)
多重均值比较
Tukey <- TukeyHSD(anova)
cld <- multcompLetters4(anova, Tukey)
整理数据
dt <- group_by(CO2, conc, Type, Treatment) %>%
summarise(uptake_mean=mean(uptake), sd=sd(uptake)) %>%
arrange(desc(uptake_mean))
cld <- as.data.frame.list(cld$`factor(conc):Type:Treatment`)
dt$Tukey <- cld$Letters
> dt
# A tibble: 28 × 6
# Groups: conc, Type [14]
conc Type Treatment uptake_mean sd Tukey
<dbl> <fct> <fct> <dbl> <dbl> <chr>
1 1000 Quebec nonchilled 43.2 3.06 a
2 675 Quebec nonchilled 41.5 2.35 a
3 1000 Quebec chilled 40.8 1.91 ab
4 350 Quebec nonchilled 40.4 2.75 ab
5 500 Quebec nonchilled 39.6 3.90 abc
6 675 Quebec chilled 37.5 2.10 abcd
7 250 Quebec nonchilled 37.4 2.76 abcd
8 500 Quebec chilled 36.7 3.61 abcde
9 350 Quebec chilled 35.8 2.62 abcde
10 250 Quebec chilled 34.5 3.93 abcde
dt列出了三个因素的列，即uptake的平均值、标准差和表示平均值之间显着差异的字母
绘制基础图形
ggplot(dt,aes(x = factor(conc), y = uptake_mean, fill = Type:Treatment)) +
geom_bar(stat = "identity", position = "dodge") +
geom_errorbar(aes(ymax = uptake_mean + sd, ymin = uptake_mean - sd),
position = position_dodge(0.9), width = 0.25, color = "Gray25") +
xlab(expression(CO[2]~Concentration~'('~mL~L^-1~')')) +
ylab(expression(CO[2]~Uptake~'('~µmol~m^2~s^-1~')')) +
scale_fill_brewer(palette = "Greens") +
theme_few()
分面设置
ggplot(dt, aes(x = factor(conc), y = uptake_mean, fill = Treatment)) +
geom_bar(stat = "identity", position = "dodge") +
geom_errorbar(aes(ymax = uptake_mean + sd, ymin = uptake_mean - sd),
position = position_dodge(0.9), width = 0.25, color = "Gray25") +
xlab(expression(CO[2]~Concentration~'('~mL~L^-1~')')) +
ylab(expression(CO[2]~Uptake~'('~µmol~m^2~s^-1~')')) +
scale_fill_brewer(palette = "Greens") +
theme_few() +
facet_grid(.~Type, labeller = label_both)
调整细节
p <- ggplot(dt, aes(x = factor(conc), y = uptake_mean, fill = Treatment)) +
geom_bar(stat = "identity", position = "dodge") +
geom_errorbar(aes(ymax = uptake_mean + sd, ymin = uptake_mean - sd),
position = position_dodge(0.9), width = 0.25, color = "Gray25") +
xlab(expression(CO[2]~Concentration~'('~mL~L^-1~')')) +
ylab(expression(CO[2]~Uptake~'('~µmol~m^2~s^-1~')')) +
theme_few() +
theme(legend.position = c(0.95,0.98),legend.justification = c(1, 1),
legend.title = element_blank(),
axis.text=element_text(color="black"),
axis.title = element_text(color="black")) +
scale_fill_manual(values = c("#C1D5A5", "#84A17C")) +
scale_y_continuous(expand = expansion(0),limits = c(0,50),breaks = seq(0,50,5))+
facet_grid(.~Type, labeller = label_both) +
geom_text(aes(label=Tukey, y = uptake_mean + sd + 2), size = 3, color = "Gray25",
show.legend = FALSE,position = position_dodge(0.9))
tag_facet(p, fontface = 1, tag_pool = c("(a) Quebec",
"(b) Mississipi"),
open = NULL, close = NULL, hjust = -0.05)
