The Ultimate Guide to (one-way) ANOVA in R: Code, Graphs & Write-Ups

Introduction

Analysis of variance (ANOVA) is a popular statistical technique for comparing outcomes across multiple groups, but utilizing it effectively in R can be overwhelming with the multitude of functions and options available. This article serves as a personal playbook - showcasing the essential tools and workflows I’ve curated specifically for conducting (one-way) ANOVA in R, from coding to visualization and reporting.

Load required packages

Code

pacman::p_load(tidyverse, ggpubr, rstatix, emmeans, performance)

The data

In this hypothetical scenario¹, we have a clinical trial involving eighteen participants who were randomly assigned to one of three drug groups: placebo, anxifree, or joyzepam. The aim of the trial was to determine if either of the two drugs could improve mood compared to the placebo. Each group consisted of six participants and the primary outcome was the change in mood over the course of the trial:

id	drug	mood_gain
1	placebo	0.5
2	placebo	0.3
3	placebo	0.1
4	anxifree	0.6
5	anxifree	0.4
6	anxifree	0.2
7	joyzepam	1.4
8	joyzepam	1.7
9	joyzepam	1.3
10	placebo	0.6
11	placebo	0.9
12	placebo	0.3
13	anxifree	1.1
14	anxifree	0.8
15	anxifree	1.2
16	joyzepam	1.8
17	joyzepam	1.3
18	joyzepam	1.4

Plotting the data

Let’s plot the data using the ggboxplot() function from the ggpubr package:

Code

data %>%
  ggpubr::ggboxplot(
    y = "mood_gain",
    x = "drug",
    color = "drug",
    add = c("dotplot"),
    add.params = list(
      position = position_jitter(w = 0.09, h = 0.04)
    ),
    palette = c("#00AFBB", "#E7B800", "#FC4E07")
  ) +
  theme(
    legend.position = "none",
    axis.title.x = element_blank()
  )

Straight away, we can see that there appears to be a difference between the drug groups, with joyzepam showing higher mood gain on average than anxifree and placebo. However, while a raw difference in groups provides a sense of the magnitude of an effect, we need statistics to helps us determine if the difference is likely to represent a real effect rather than just chance variation.

To assess whether the observed differences between the three treatment groups are statistically significant, we will perform a one-way ANOVA².

One-way ANOVA in R

My preferred way to compute (one-way) ANOVA in R is with the rstatix package (Kassambara 2023).

The first step is to calculate the ANOVA with the rstatix::anova_test function:

Code

rstatixaov <-
  data %>%
  rstatix::anova_test(dv = mood_gain, between = drug, wid = id)

Then print the output:

Code

rstatixaov %>% 
  tibble::as_tibble() %>% 
  mutate(p = ifelse(p < 0.001, "< .001", signif(p, 2))) %>% 
  mutate(across(where(is.numeric), round,2)) %>% 
  select(-`p<.05`) %>% 
  flextable::flextable() %>% 
  flextable::autofit()

Effect	DFn	DFd	F	p	ges
drug	2	15	18.61	< .001	0.71

The results reveal a statistically significant effect of the drug treatments on changing mood. The small p-value (<0.001) indicates it is highly unlikely (< 0.1% chance) to obtain these results if the drugs had no real effect. Furthermore, the large generalized eta squared value (ges = 0.71) suggests the drug treatments accounted for 71% of the variance in mood scores. This large effect size provides evidence that the drugs substantially impacted mood.

While the overall analysis provides evidence for drug effect, we don’t yet know which specific drug treatments led to superior mood gains compared to others. Our previous visualization gives us some hints, but to explore formally, we need to run post-hoc pairwise comparisons.

Post-hoc comparisons

Post-hoc comparisons are statistical tests performed to determine which specific pairs of drug groups differ from each other significantly, while ensuring the overall risk of false positives remains low³.

I like to run pairwise comparisons with the rstatix::emmeans_test() function:

Code

rstatixcomp <-
  rstatix::emmeans_test(mood_gain ~ drug,
    data = data,
    p.adjust.method = "holm",
    detailed = T
  ) %>%
  mutate(across(where(is.numeric), signif, 2)) %>%
  select(-.y., -term, -se, -null.value, -df, -p) %>% 
  arrange(p.adj)

rstatixcomp %>% 
  flextable::flextable() %>%
  flextable::autofit()

group1	group2	estimate	conf.low	conf.high	statistic	p.adj	p.adj.signif
joyzepam	placebo	1.00	0.66	1.40	5.9	0.000091	****
anxifree	joyzepam	-0.77	-1.10	-0.39	-4.4	0.001100	**
anxifree	placebo	0.27	-0.11	0.64	1.5	0.150000	ns

The results provide evidence of a difference in mood gain between joyzepam and anxifree (diff = 0.77), as well as between joyzepam and placebo (diff = 1), but not between anxifree and placebo (diff = 0.27).

Testing assumptions

To ensure the validity of the results from a one-way ANOVA, several key assumptions must be tested. Firstly, the residuals are assumed to be normally distributed, often assessed visually through residual plots. Secondly, the homogeneity of variances assumption requires equal population variances across groups, typically checked with Levene’s test and residual plots. Thirdly, observations must be independent, unrelated, and uninfluenced by each other, usually satisfied through proper experimental design and randomization (Navarro 2019).

The trial participants were not related (i think) and randomly allocated to each group so it is safe to assume each observation is independent.

However, let’s check the normality and equal variance assumptions using the performance::check_models function:

Code

# homogenity of variance
assumptions <-
  performance::check_model(lm(mood_gain ~ drug, data = data),
    check = c("homogeneity", "normality")
  )

assumptions

Here we see no major violations to the equal variance and normality assumptions, giving us confidence that our results are valid.

Putting it all together

My favourite aspect of rstatix is the ability to incorporate the statistical results into plots produced by the ggpubr package⁴.

For example, here is how to put all the most relevant information into one plot:

Code

# add xy position to rstatix dataframe
stat_test <-
  rstatixcomp %>%
  rstatix::add_xy_position(x = "drug")

# boxplot using ggpubr package add means, ANOVA, and pairwise results
data %>%
  ggpubr::ggboxplot(
    y = "mood_gain",
    x = "drug",
    color = "drug",
    add = c("dotplot"),
    add.params = list(
      position = position_jitter(w = 0.09, h = 0.04),
      size = .8
    ),
    palette = c("#00AFBB", "#E7B800", "#FC4E07")
  ) +
  scale_y_continuous(breaks = seq(.5, 2, .5)) +
  coord_cartesian(ylim = c(0.1, 2.3)) +
  theme(
    legend.position = "none",
    axis.title.x = element_blank()
  ) +
  labs(y = "Mood Gain") +

    stat_summary(
    geom = "point",
    fun.y = "mean",
    col = "black",
    size = 2.5,
    shape = 24,
  ) +
  ggpubr::stat_pvalue_manual(stat_test,
    label = "p.adj.signif",
    tip.length = 0.01,
    hide.ns = T,
    bracket.nudge.y = -.03
  ) +
  ggpubr::stat_anova_test(label = "as_detailed_italic", label.y.npc = 1)

Writing up results

Here is an example for writing up ANOVA and pairwise comparison results for publication in a scientific journal:

“A one-way ANOVA revealed a significant difference in mood gain between the three groups, F(2,15) = 18.611, p < 0.01, η² = 0.713. Post hoc pairwise comparisons were conducted using the Holm method to control for the family-wise error rate across the multiple comparisons. The Holm-adjusted p-values indicated that the joyzepan group had significantly higher mood gain (M = 1.483, SD = 0.214) compared to the placebo (M = 0.45, SD = 0.281) and anxifree (M = 0.717, SD = 0.392) group. There was no significant difference between the anxifree group and the placebo group.”

Final thoughts

With so many options available running ANOVA in R can seem daunting at first, but breaking it down into manageable steps makes the process much more approachable. Leveraging the rstatix package is my preferred approach as it not only simplifies the calculations required for ANOVA but also provides a range of functions to help communicate the results effectively.

References

Kassambara, Alboukadel. 2023. “Rstatix: Pipe-Friendly Framework for Basic Statistical Tests.” https://CRAN.R-project.org/package=rstatix.

Navarro, Danielle. 2019. “Learning Statistics with r: A Tutorial for Psychology Students and Other Beginners (Version 0.6.1).” University of New South Wales.

Footnotes

This data set is taken from an article explaining how to calculate ANOVA in jamovi. This allowed me to directly compare the results generated by R, with those that can be recreated using the jamovi software. The data set can be accessed here.↩︎
Read more about ANOVA here ↩︎
Controlling the overall Type I error rate across the multiple comparisons.↩︎
A good article on how to incorporate pairwise comparisons into plots can be found here ↩︎

--- title: "The Ultimate Guide to (one-way) ANOVA in R: Code, Graphs & Write-Ups" execute: echo: true warning: false code-fold: show bibliography: references.bib --- # Introduction Analysis of variance (ANOVA) is a popular statistical technique for comparing outcomes across multiple groups, but utilizing it effectively in R can be overwhelming with the multitude of functions and options available. This article serves as a personal playbook - showcasing the essential tools and workflows I've curated specifically for conducting (one-way) ANOVA in R, from coding to visualization and reporting. # Load required packages ```{r} pacman::p_load(tidyverse, ggpubr, rstatix, emmeans, performance) ``` # The data ```{r echo=FALSE} data <- tibble::tribble( ~id, ~drug, ~cbt, ~mood_gain, 1L, "placebo", "no.therapy", 0.5, 2L, "placebo", "no.therapy", 0.3, 3L, "placebo", "no.therapy", 0.1, 4L, "anxifree", "no.therapy", 0.6, 5L, "anxifree", "no.therapy", 0.4, 6L, "anxifree", "no.therapy", 0.2, 7L, "joyzepam", "no.therapy", 1.4, 8L, "joyzepam", "no.therapy", 1.7, 9L, "joyzepam", "no.therapy", 1.3, 10L, "placebo", "CBT", 0.6, 11L, "placebo", "CBT", 0.9, 12L, "placebo", "CBT", 0.3, 13L, "anxifree", "CBT", 1.1, 14L, "anxifree", "CBT", 0.8, 15L, "anxifree", "CBT", 1.2, 16L, "joyzepam", "CBT", 1.8, 17L, "joyzepam", "CBT", 1.3, 18L, "joyzepam", "CBT", 1.4 ) %>% mutate(across(c(id, drug, cbt), as.factor)) %>% select(-cbt) ``` In this hypothetical scenario[^1], we have a clinical trial involving eighteen participants who were randomly assigned to one of three drug groups: placebo, anxifree, or joyzepam. The aim of the trial was to determine if either of the two drugs could improve mood compared to the placebo. Each group consisted of six participants and the primary outcome was the change in mood over the course of the trial: [^1]: This data set is taken from an [article](https://danawanzer.github.io/stats-with-jamovi/one-way-anova.html) explaining how to calculate ANOVA in jamovi. This allowed me to directly compare the results generated by R, with those that can be recreated using the jamovi software. The data set can be accessed [here](https://lsj.readthedocs.io/en/latest/_static/data/clinicaltrial.omv). ```{r warning=FALSE, message=FALSE, echo=FALSE} data %>% flextable::flextable() %>% flextable::autofit() ``` # Plotting the data Let's plot the data using the `ggboxplot()` function from the `ggpubr` package: ```{r} data %>% ggpubr::ggboxplot( y = "mood_gain", x = "drug", color = "drug", add = c("dotplot"), add.params = list( position = position_jitter(w = 0.09, h = 0.04) ), palette = c("#00AFBB", "#E7B800", "#FC4E07") ) + theme( legend.position = "none", axis.title.x = element_blank() ) ``` Straight away, we can see that there appears to be a difference between the drug groups, with joyzepam showing higher mood gain on average than anxifree and placebo. However, while a raw difference in groups provides a sense of the magnitude of an effect, we need statistics to helps us determine if the difference is likely to represent a real effect rather than just chance variation. To assess whether the observed differences between the three treatment groups are statistically significant, we will perform a one-way ANOVA[^2]. [^2]: Read more about ANOVA [here](https://stats.libretexts.org/Bookshelves/Introductory_Statistics/Mostly_Harmless_Statistics_(Webb)/11%3A_Analysis_of_Variance/11.01%3A_One-Way_ANOVA) # One-way ANOVA in R My preferred way to compute (one-way) ANOVA in R is with the `rstatix` package [@rstatix]. The first step is to calculate the ANOVA with the `rstatix::anova_test` function: ```{r} rstatixaov <- data %>% rstatix::anova_test(dv = mood_gain, between = drug, wid = id) ``` Then print the output: ```{r} rstatixaov %>% tibble::as_tibble() %>% mutate(p = ifelse(p < 0.001, "< .001", signif(p, 2))) %>% mutate(across(where(is.numeric), round,2)) %>% select(-`p<.05`) %>% flextable::flextable() %>% flextable::autofit() ``` The results reveal a statistically significant effect of the drug treatments on changing mood. The small p-value (\<0.001) indicates it is highly unlikely (\< 0.1% chance) to obtain these results if the drugs had no real effect. Furthermore, the large generalized eta squared value (ges = 0.71) suggests the drug treatments accounted for 71% of the variance in mood scores. This large effect size provides evidence that the drugs substantially impacted mood. While the overall analysis provides evidence for drug effect, we don't yet know which specific drug treatments led to superior mood gains compared to others. Our previous visualization gives us some hints, but to explore formally, we need to run post-hoc pairwise comparisons. # Post-hoc comparisons Post-hoc comparisons are statistical tests performed to determine which specific pairs of drug groups differ from each other significantly, while ensuring the overall risk of false positives remains low[^3]. [^3]: Controlling the overall Type I error rate across the multiple comparisons. I like to run pairwise comparisons with the `rstatix::emmeans_test()` function: ```{r} rstatixcomp <- rstatix::emmeans_test(mood_gain ~ drug, data = data, p.adjust.method = "holm", detailed = T ) %>% mutate(across(where(is.numeric), signif, 2)) %>% select(-.y., -term, -se, -null.value, -df, -p) %>% arrange(p.adj) rstatixcomp %>% flextable::flextable() %>% flextable::autofit() ``` The results provide evidence of a difference in mood gain between joyzepam and anxifree (diff = 0.77), as well as between joyzepam and placebo (diff = 1), but not between anxifree and placebo (diff = 0.27). # Testing assumptions To ensure the validity of the results from a one-way ANOVA, several key assumptions must be tested. Firstly, the residuals are assumed to be normally distributed, often assessed visually through residual plots. Secondly, the homogeneity of variances assumption requires equal population variances across groups, typically checked with Levene's test and residual plots. Thirdly, observations must be independent, unrelated, and uninfluenced by each other, usually satisfied through proper experimental design and randomization [@navarro2015learning]. The trial participants were not related (i think) and randomly allocated to each group so it is safe to assume each observation is independent. However, let's check the normality and equal variance assumptions using the `performance::check_models` function: ```{r} # homogenity of variance assumptions <- performance::check_model(lm(mood_gain ~ drug, data = data), check = c("homogeneity", "normality") ) assumptions ``` Here we see no major violations to the equal variance and normality assumptions, giving us confidence that our results are valid. # Putting it all together My favourite aspect of `rstatix` is the ability to incorporate the statistical results into plots produced by the `ggpubr` package[^4]. [^4]: A good article on how to incorporate pairwise comparisons into plots can be found [here](https://www.datanovia.com/en/blog/how-to-add-p-values-onto-basic-ggplots/) For example, here is how to put all the most relevant information into one plot: ```{r warning=FALSE} # add xy position to rstatix dataframe stat_test <- rstatixcomp %>% rstatix::add_xy_position(x = "drug") # boxplot using ggpubr package add means, ANOVA, and pairwise results data %>% ggpubr::ggboxplot( y = "mood_gain", x = "drug", color = "drug", add = c("dotplot"), add.params = list( position = position_jitter(w = 0.09, h = 0.04), size = .8 ), palette = c("#00AFBB", "#E7B800", "#FC4E07") ) + scale_y_continuous(breaks = seq(.5, 2, .5)) + coord_cartesian(ylim = c(0.1, 2.3)) + theme( legend.position = "none", axis.title.x = element_blank() ) + labs(y = "Mood Gain") + stat_summary( geom = "point", fun.y = "mean", col = "black", size = 2.5, shape = 24, ) + ggpubr::stat_pvalue_manual(stat_test, label = "p.adj.signif", tip.length = 0.01, hide.ns = T, bracket.nudge.y = -.03 ) + ggpubr::stat_anova_test(label = "as_detailed_italic", label.y.npc = 1) ``` # Writing up results Here is an example for writing up ANOVA and pairwise comparison results for publication in a scientific journal: ```{r echo=FALSE, include=FALSE} ngroups <- as.character(english::as.english(rstatixaov$DFn + 1)) df1 <- rstatixaov$DFn df2 <- rstatixaov$DFd fval <- rstatixaov$F ges <- rstatixaov$ges pval <- ifelse(rstatixaov$p < 0.01, "< 0.01", paste0("= ", as.character(rstatixaov$p))) sig_write_up <- ifelse(rstatixaov$p < 0.05, "a significant difference", "no significant difference") f_write_up <- paste0("F(", df1, ",", df2, ") = ", fval, ", p ", pval, ", η² = ", ges, ".") rstatixsummary <- data %>% group_by(drug) %>% rstatix::get_summary_stats( type = c("common") ) %>% select(drug, mean, sd) ``` "A one-way ANOVA revealed `r sig_write_up` in mood gain between the `r ngroups` groups, `r f_write_up` Post hoc pairwise comparisons were conducted using the Holm method to control for the family-wise error rate across the multiple comparisons. The Holm-adjusted p-values indicated that the joyzepan group had significantly higher mood gain (M = `r rstatixsummary[2,2]`, SD = `r rstatixsummary[2,3]`) compared to the placebo (M = `r rstatixsummary[3,2]`, SD = `r rstatixsummary[3,3]`) and anxifree (M = `r rstatixsummary[1,2]`, SD = `r rstatixsummary[1,3]`) group. There was no significant difference between the anxifree group and the placebo group." # Final thoughts With so many options available running ANOVA in R can seem daunting at first, but breaking it down into manageable steps makes the process much more approachable. Leveraging the `rstatix` package is my preferred approach as it not only simplifies the calculations required for ANOVA but also provides a range of functions to help communicate the results effectively.