> library(multcompView)
> library(stats)
> library(tidyverse)
> library(ggsci)
>
> data2 <- ToothGrowth
> data2$dose = as.factor(data2$dose)
>
> anova <- aov(len ~ supp*dose, data = data2) #two way anova
> summary(anova) #represents ANOVA
Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.4 205.4 15.572 0.000231 ***
dose 2 2426.4 1213.2 92.000 < 2e-16 ***
supp:dose 2 108.3 54.2 4.107 0.021860 *
Residuals 54 712.1 13.2
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> tukey <- TukeyHSD(anova)
>
> group_lettering <- multcompLetters4(anova, tukey) #
> group_lettering
$supp
OJ VC
"a" "b"
$dose
2 1 0.5
"a" "b" "c"
$`supp:dose`
VC:2 OJ:2 OJ:1 VC:1 OJ:0.5 VC:0.5
"a" "a" "a" "b" "b" "c"
> group_lettering2 <- data.frame(group_lettering$`supp:dose`$Letters)
> group_lettering2
group_lettering..supp.dose..Letters
VC:2 a
OJ:2 a
OJ:1 a
VC:1 b
OJ:0.5 b
VC:0.5 c
>
> mean_data2 <- data2 %>%
+ group_by(supp, dose) %>%
+ summarise(len_mean=mean(len), sd = sd(len)) %>%
+ arrange(desc(len_mean))
`summarise()` has grouped output by 'supp'. You can override using
the `.groups` argument.
>
> mean_data2$group_lettering <- group_lettering2$group_lettering..supp.dose..Letters
>
>
> ggplot(mean_data2,aes(x = dose, y = len_mean,group=supp))+
+ geom_bar(position=position_dodge(0.9),stat = "identity", aes(fill = supp),
+ show.legend = TRUE)
>
> ggplot(mean_data2, aes(x = dose, y = len_mean,group=supp))+
+ geom_bar(position=position_dodge(0.9),stat = "identity",
+ aes(fill = supp),show.legend = TRUE) +
+ geom_errorbar(aes(ymin = len_mean-sd, ymax=len_mean+sd),width = 0.1,
+ position=position_dodge(0.9))
>
> ggplot(mean_data2, aes(x = dose, y = len_mean,group=supp))+
+ geom_bar(position=position_dodge(0.9),stat = "identity",
+ aes(fill = supp), show.legend = TRUE) +
+ geom_errorbar(aes(ymin = len_mean-sd, ymax=len_mean+sd),
+ width = 0.1, position=position_dodge(0.9))+
+ geom_text(aes(label = group_lettering, y = len_mean + sd),
+ vjust=-0.4,position=position_dodge(0.9))
>
>
> ggplot(mean_data2, aes(x = dose, y = len_mean,group=supp))+
+ geom_bar(position=position_dodge(0.9),stat = "identity",
+ aes(fill = supp), show.legend = TRUE) +
+ geom_errorbar(aes(ymin = len_mean-sd, ymax=len_mean+sd),
+ width = 0.1, position=position_dodge(0.9)) +
+ geom_text(aes(label = group_lettering, y = len_mean + sd),
+ vjust=-0.4, position=position_dodge(0.9)) +
+ scale_y_continuous(expand = expansion(0),limits = c(0,35),
+ breaks = seq(0,35,5))+
+ labs(x=NULL,y=NULL)+
+ theme_minimal()+
+ theme(
+ plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
+ panel.background = element_blank(),
+ axis.line = element_line(color = "black"),
+ axis.title = element_text(size = 10, color = "black",face = "bold"),
+ axis.text = element_text(size = 10,color = "black"),
+ axis.text.x = element_text(margin=margin(t =3)),
+ axis.text.y = element_text(size = 10),
+ axis.title.y = element_text(margin = margin(r = 10)),
+ axis.ticks.x = element_blank())+
+ scale_fill_jco()
> library(tidyverse)
> library(ggthemes)
> library(multcompView)
> library(egg)
> df <- read_tsv("co2.xls")
Rows: 84 Columns: 5
── Column specification ───────────────────────────────────────────
Delimiter: "\t"
chr (3): Plant, Type, Treatment
dbl (2): conc, uptake
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
> anova <- aov(uptake ~ factor(conc)*Type*Treatment, data = df)
> summary(anova)
Df Sum Sq Mean Sq F value Pr(>F)
factor(conc) 6 4069 678 80.548 < 2e-16 ***
Type 1 3366 3366 399.758 < 2e-16 ***
Treatment 1 988 988 117.368 2.32e-15 ***
factor(conc):Type 6 374 62 7.412 7.24e-06 ***
factor(conc):Treatment 6 101 17 1.999 0.0811 .
Type:Treatment 1 226 226 26.812 3.15e-06 ***
factor(conc):Type:Treatment 6 112 19 2.216 0.0547 .
Residuals 56 471 8
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 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))
`summarise()` has grouped output by 'conc', 'Type'. You can
override using the `.groups` argument.
> cld <- as.data.frame.list(cld$`factor(conc):Type:Treatment`)
> dt$Tukey <- cld$Letters
> 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)
