Three-Way ANOVA in R: Three Factors and Their Interaction

Test three grouping factors at once — their main effects, every two-way interaction, and the three-way interaction — then decompose a significant three-way interaction step by step

Learn to run a three-way ANOVA in R to test three factors and their interactions on a numeric outcome. Use rstatix anova_test() for a tidy table with a ges effect size, and base R aov() + summary() for the classic layout — with the outlier, normality and equal-variance checks, how to read the seven terms, and the full decomposition of a significant three-way interaction into simple two-way interactions, simple simple main effects and simple simple comparisons, with the p-values on the plot.

TipKey takeaways

• A three-way ANOVA tests three grouping factors at once — each factor’s main effect, all three two-way interactions, and the three-way interaction — on a single numeric outcome.
• Use anova_test() from rstatix with an outcome ~ A * B * C formula for a tidy table (F, p, and the ges effect size for each of the seven terms), or base R aov() + summary().
• The three-way interaction (A:B:C) is the headline: a significant one means the two-way interaction between two factors itself changes across the levels of the third — interpret it before anything else.
• Check the same assumptions as the lower-order tests: no extreme outliers, normal residuals (shapiro_test() + Q-Q), and equal variance across cells (levene_test()).
• Decompose a significant three-way interaction in three steps: simple two-way interactions → simple simple main effects → simple simple pairwise comparisons — then plot the result with the p-values on it.

Introduction

A three-way ANOVA (also called a three-factor or factorial ANOVA) extends the two-way version by one more factor: it asks how three categorical variables — and all their combinations — affect a numeric outcome. On top of the three main effects, the model now estimates three two-way interactions (A:B, A:C, B:C) and one three-way interaction (A:B:C).

Use it when each observation is classified by three grouping variables and you suspect the effect of one depends on the other two. The classic scenario: a clinical team tests whether a treatment works the same way for everyone, or whether its effect depends on patients’ sex and their risk level — exactly the kind of question only a three-way design can answer. This lesson runs the whole workflow two ways — the modern, pipe-friendly rstatix anova_test() and the classic base-R aov() — using the headache dataset, and shows the assumption checks plus the full step-by-step breakdown of a significant three-way interaction.

Note

Three-way ANOVA hypotheses. There is a null hypothesis for every term: the three main effects (each factor’s group means are equal), the three two-way interactions (each pair of factors acts independently), and — the one that matters most — the three-way interaction (H₀: the two-way interaction between any two factors is the same at every level of the third). A significant three-way interaction rejects that last null: the factors do not combine additively.

The data: three factors at once

We use the headache dataset from the datarium package — migraine headache pain_score for 72 participants, classified by three factors: gender (male, female), risk of migraine (low, high), and treatment (X, Y, Z). The research question: does the effect of treatment on pain depend on both gender and risk? Here treatment is the focal variable (our primary concern) and gender and risk are the moderators.

library(rstatix)

set.seed(123)
data("headache", package = "datarium")

# one random row per gender × risk × treatment cell
headache %>% sample_n_by(gender, risk, treatment, size = 1)

# A tibble: 12 × 5
      id gender risk  treatment pain_score
   <int> <fct>  <fct> <fct>          <dbl>
 1    21 male   high  X               92.7
 2    30 male   high  Y               80.1
 3    33 male   high  Z               81.3
 4     2 male   low   X               76.8
 5     8 male   low   Y               80.7
 6    18 male   low   Z               74.7
 7    57 female high  X               68.4
 8    65 female high  Y               82.1
 9    70 female high  Z               80.4
10    42 female low   X               78.1
11    48 female low   Y               64.1
12    49 female low   Z               69.0

gender (2) × risk (2) × treatment (3) make a 2 × 2 × 3 factorial design, and the outcome is pain_score. Summarize each cell with get_summary_stats() — the mean and standard deviation per combination, the numbers behind the plot that follows:

library(rstatix)

data("headache", package = "datarium")

headache %>%
  group_by(gender, risk, treatment) %>%
  get_summary_stats(pain_score, type = "mean_sd")

# A tibble: 12 × 7
   gender risk  treatment variable       n  mean    sd
   <fct>  <fct> <fct>     <fct>      <dbl> <dbl> <dbl>
 1 male   high  X         pain_score     6  92.7  5.12
 2 male   high  Y         pain_score     6  82.3  5.00
 3 male   high  Z         pain_score     6  79.7  4.05
 4 male   low   X         pain_score     6  76.1  3.86
 5 male   low   Y         pain_score     6  73.1  4.76
 6 male   low   Z         pain_score     6  74.5  4.89
 7 female high  X         pain_score     6  78.9  5.32
 8 female high  Y         pain_score     6  81.2  4.62
 9 female high  Z         pain_score     6  81.0  3.98
10 female low   X         pain_score     6  74.2  3.69
11 female low   Y         pain_score     6  68.4  4.08
12 female low   Z         pain_score     6  69.8  2.72

Look at the data first

Always plot before testing. A boxplot of pain_score by treatment, coloured by risk and faceted by gender, shows all twelve cells at a glance — and already hints that the treatment gap is widest for one corner of the design:

library(ggpubr)

data("headache", package = "datarium")

ggboxplot(
  headache, x = "treatment", y = "pain_score",
  color = "risk", palette = "jco", facet.by = "gender",
  xlab = "Treatment", ylab = "Headache pain score"
)

Run the three-way ANOVA

With rstatix (recommended)

anova_test() takes the data and an outcome ~ A * B * C formula. The * expands to every term — the three main effects, the three two-way interactions, and the three-way interaction. The result is a tidy table — one row per term:

library(rstatix)

data("headache", package = "datarium")

res.aov <- headache %>% anova_test(pain_score ~ gender * risk * treatment)
res.aov

ANOVA Table (type II tests)

                 Effect DFn DFd      F        p p<.05   ges
1                gender   1  60 16.196 1.63e-04     * 0.213
2                  risk   1  60 92.699 8.80e-14     * 0.607
3             treatment   2  60  7.318 1.00e-03     * 0.196
4           gender:risk   1  60  0.141 7.08e-01       0.002
5      gender:treatment   2  60  3.338 4.20e-02     * 0.100
6        risk:treatment   2  60  0.713 4.94e-01       0.023
7 gender:risk:treatment   2  60  7.406 1.00e-03     * 0.198

Read it from the bottom up — the highest-order term first. There is a statistically significant three-way interaction between gender, risk and treatment, F(2, 60) = 7.41, p = 0.001, η²g = 0.20. Because that top term is significant, it governs the whole interpretation: the risk:treatment interaction itself differs between males and females, so you must not stop at the main-effect rows (gender, risk, treatment are all significant, but reporting them alone would hide the real story). Note too that the lower-order gender:risk and risk:treatment interactions are not significant on their own — the structure only emerges once all three factors are crossed.

Tip

Read interactions top-down. In any factorial ANOVA, interpret the highest-order significant interaction first. A significant A:B:C makes the two-way interactions and main effects misleading in isolation, so you decompose A:B:C rather than quoting the lower rows. Only when the three-way term is not significant do you drop down to interpret the two-way interactions.

Tip

Run it on your own data. Load your file — mydata <- read.csv("my-file.csv") — make all three grouping columns factors (mydata$A <- factor(mydata$A)), then pipe it in: mydata %>% anova_test(outcome ~ A * B * C).

With base R

aov() fits the model and summary() prints the classic ANOVA table — one line per term plus the residuals:

data("headache", package = "datarium")

res.aov.base <- aov(pain_score ~ gender * risk * treatment, data = headache)
summary(res.aov.base)

                      Df Sum Sq Mean Sq F value   Pr(>F)    
gender                 1  313.4   313.4  16.196 0.000163 ***
risk                   1 1793.6  1793.6  92.699  8.8e-14 ***
treatment              2  283.2   141.6   7.318 0.001433 ** 
gender:risk            1    2.7     2.7   0.141 0.708487    
gender:treatment       2  129.2    64.6   3.338 0.042200 *  
risk:treatment         2   27.6    13.8   0.713 0.494221    
gender:risk:treatment  2  286.6   143.3   7.406 0.001334 ** 
Residuals             60 1160.9    19.3                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The rows are the three main effects, the three two-way interactions, the gender:risk:treatment line and the Residuals; F value and Pr(>F) (with significance stars) carry the same conclusion as the rstatix call. aov() keeps the fitted model, which the diagnostics and the post-hoc decomposition below reuse.

Effect size: how big is each term?

A significant p-value says a term matters; the effect size says how much. anova_test() already prints ges (generalized eta-squared) per term. As a rough benchmark, η² ≈ 0.01 is small, ≈ 0.06 medium and ≥ 0.14 large. Here risk dominates (ges ≈ 0.61, very large), the three-way interaction is a real, large effect (ges ≈ 0.20), and the non-significant two-way interactions sit near zero — report the effect size alongside each p-value, never just “p < 0.05”.

Check the assumptions

A three-way ANOVA trusts the same assumptions as the lower-order tests (beyond independent observations): no extreme outliers, the residuals are normally distributed, and the cells have equal variance.

Outliers

identify_outliers() flags extreme values in each cell:

library(rstatix)

data("headache", package = "datarium")

headache %>%
  group_by(gender, risk, treatment) %>%
  identify_outliers(pain_score)

# A tibble: 4 × 7
  gender risk  treatment    id pain_score is.outlier is.extreme
  <fct>  <fct> <fct>     <int>      <dbl> <lgl>      <lgl>     
1 female high  X            57       68.4 TRUE       TRUE      
2 female high  Y            62       73.1 TRUE       FALSE     
3 female high  Z            67       75.0 TRUE       FALSE     
4 female high  Z            71       87.1 TRUE       FALSE     

There is one extreme outlier (id 57 — a female at high risk taking drug X). With many small cells a single extreme value is common; you can keep it and check whether it changes the conclusion (re-run the ANOVA without it and compare), or use a robust ANOVA from the WRS2 package. Here it does not alter the verdict, so we keep it.

Normality of residuals — Shapiro-Wilk + Q-Q

ANOVA assumes the residuals (not each raw cell) are normal. Fit the model, test the residuals with Shapiro-Wilk (p > 0.05 → no departure from normal), and back it with a Q-Q plot:

library(rstatix)

data("headache", package = "datarium")

model <- lm(pain_score ~ gender * risk * treatment, data = headache)
shapiro_test(residuals(model))

# A tibble: 1 × 3
  variable         statistic p.value
  <chr>                <dbl>   <dbl>
1 residuals(model)     0.982   0.398

library(ggpubr)

data("headache", package = "datarium")
model <- lm(pain_score ~ gender * risk * treatment, data = headache)

ggqqplot(residuals(model))

The points fall along the reference line and Shapiro-Wilk is non-significant (p = 0.4) — the residuals are normal. You can also check normality within each cell and a faceted Q-Q plot:

library(rstatix)

data("headache", package = "datarium")

headache %>%
  group_by(gender, risk, treatment) %>%
  shapiro_test(pain_score)

# A tibble: 12 × 6
   gender risk  treatment variable   statistic       p
   <fct>  <fct> <fct>     <chr>          <dbl>   <dbl>
 1 male   high  X         pain_score     0.958 0.808  
 2 male   high  Y         pain_score     0.902 0.384  
 3 male   high  Z         pain_score     0.955 0.784  
 4 male   low   X         pain_score     0.982 0.962  
 5 male   low   Y         pain_score     0.920 0.507  
 6 male   low   Z         pain_score     0.924 0.535  
 7 female high  X         pain_score     0.714 0.00869
 8 female high  Y         pain_score     0.939 0.654  
 9 female high  Z         pain_score     0.971 0.901  
10 female low   X         pain_score     0.933 0.600  
11 female low   Y         pain_score     0.927 0.555  
12 female low   Z         pain_score     0.958 0.801  

library(ggpubr)

data("headache", package = "datarium")

ggqqplot(headache, "pain_score", ggtheme = theme_bw()) +
  facet_grid(gender + risk ~ treatment, labeller = "label_both")

Scores are normal in every cell (p > 0.05) except one (female, high risk, drug X, p = 0.0086) — the same cell that held the extreme outlier.

Homogeneity of variance — Levene’s test

Levene’s test asks whether all twelve cells share a common variance. Cross all three factors so it tests every cell; p > 0.05 means no evidence the variances differ — what we want:

library(rstatix)

data("headache", package = "datarium")

headache %>% levene_test(pain_score ~ gender * risk * treatment)

# A tibble: 1 × 4
    df1   df2 statistic     p
  <int> <int>     <dbl> <dbl>
1    11    60     0.179 0.998

Levene’s test is non-significant, so equal-variance holds across the cells.

Note

If the residuals stray far from the line, or Levene’s test is significant, the data don’t fit the classic three-way ANOVA. See the dedicated normality test lesson for the full assumption workflow and the robust/non-parametric routes.

Decompose the significant three-way interaction

Because the gender:risk:treatment interaction is significant, the right follow-up is a structured, three-step decomposition (Keppel & Wickens). At each step we pass the overall model as error so the analyses borrow the pooled error term from the full three-way model — more powerful when, as here, the equal-variance assumption holds.

Step 1 — simple two-way interactions

First, run the risk:treatment two-way interaction at each level of gender. Group by gender, fit pain_score ~ risk * treatment, and pass error = model:

library(rstatix)

data("headache", package = "datarium")

model <- lm(pain_score ~ gender * risk * treatment, data = headache)
headache %>%
  group_by(gender) %>%
  anova_test(pain_score ~ risk * treatment, error = model)

# A tibble: 6 × 8
  gender Effect           DFn   DFd      F             p `p<.05`   ges
* <fct>  <chr>          <dbl> <dbl>  <dbl>         <dbl> <chr>   <dbl>
1 male   risk               1    60 50.0   0.00000000187 "*"     0.455
2 male   treatment          2    60 10.2   0.000157      "*"     0.253
3 male   risk:treatment     2    60  5.25  0.008         "*"     0.149
4 female risk               1    60 42.8   0.000000015   "*"     0.416
5 female treatment          2    60  0.482 0.62          ""      0.016
6 female risk:treatment     2    60  2.87  0.065         ""      0.087

There is a statistically significant simple two-way interaction between risk and treatment for males, F(2, 60) = 5.25, p = 0.008, but not for females, F(2, 60) = 2.87, p = 0.065. For males, the effect of treatment on pain depends on their migraine risk.

Note

Adjust your alpha. Statistical significance of a simple two-way interaction is judged at a Bonferroni-adjusted level of 0.05 / 2 = 0.025 (two simple two-way interactions, one per gender) — so the male interaction (p = 0.008) clears the bar and the female one (p = 0.065) does not.

Step 2 — simple simple main effects

A significant simple two-way interaction is followed up with simple simple main effects. We only need this for males (the only significant simple two-way interaction): group by gender and risk, and test the effect of treatment — again with error = model. Use base R subset() to read off the males:

library(rstatix)

data("headache", package = "datarium")

model <- lm(pain_score ~ gender * risk * treatment, data = headache)
treatment.effect <- headache %>%
  group_by(gender, risk) %>%
  anova_test(pain_score ~ treatment, error = model)

subset(treatment.effect, gender == "male")   # base-R filter

# A tibble: 2 × 9
  gender risk  Effect      DFn   DFd     F         p `p<.05`   ges
  <fct>  <fct> <chr>     <dbl> <dbl> <dbl>     <dbl> <chr>   <dbl>
1 male   high  treatment     2    60 14.8  0.0000061 "*"     0.33 
2 male   low   treatment     2    60  0.66 0.521     ""      0.022

There is a significant simple simple main effect of treatment for males at high risk, F(2, 60) = 14.8, p < 0.0001, but not for males at low risk, F(2, 60) = 0.66, p = 0.521 (judged at the Bonferroni level 0.05 / 2 = 0.025). So the type of treatment matters for males who are at high risk — and not otherwise.

Step 3 — simple simple pairwise comparisons

Finally, find which treatments differ for males at high risk with pairwise comparisons. emmeans_test() (estimated marginal means, a wrapper around the emmeans package) runs all pairwise comparisons with a Bonferroni adjustment; index the rows we need with base-R [ ]:

library(rstatix)
library(emmeans)

data("headache", package = "datarium")

pwc <- headache %>%
  group_by(gender, risk) %>%
  emmeans_test(pain_score ~ treatment, p.adjust.method = "bonferroni")

# males at high risk only
pwc[pwc$gender == "male" & pwc$risk == "high", c("group1", "group2", "p.adj", "p.adj.signif")]

# A tibble: 3 × 4
  group1 group2      p.adj p.adj.signif
  <chr>  <chr>       <dbl> <chr>       
1 X      Y      0.000386   ***         
2 X      Z      0.00000942 ****        
3 Y      Z      0.897      ns          

For males at high risk, treatment X differs significantly from Y — a mean difference of 10.4 pain points (p.adj < 0.001) — and X from Z by 13.1 points (p.adj < 0.0001), while Y and Z do not differ (difference 2.66, p.adj = 0.897). In words: drugs Y and Z both lower pain substantially relative to X for high-risk males, and are statistically indistinguishable from each other.

Plot it with the p-values on it

The publication-ready figure puts the result on the plot — the overall three-way test in the subtitle (get_test_label()) and the significant pairwise brackets (males at high risk) via stat_pvalue_manual(), after add_xy_position() computes where the brackets sit:

library(rstatix)
library(ggpubr)
library(emmeans)

data("headache", package = "datarium")

res.aov <- headache %>% anova_test(pain_score ~ gender * risk * treatment)
pwc <- headache %>%
  group_by(gender, risk) %>%
  emmeans_test(pain_score ~ treatment, p.adjust.method = "bonferroni") %>%
  add_xy_position(x = "treatment")
pwc.male.high <- pwc[pwc$gender == "male" & pwc$risk == "high", ]

ggboxplot(
  headache, x = "treatment", y = "pain_score",
  color = "risk", palette = "jco", facet.by = "gender",
  xlab = "Treatment", ylab = "Headache pain score"
) +
  stat_pvalue_manual(
    pwc.male.high, color = "risk", linetype = "risk", hide.ns = TRUE,
    tip.length = 0, step.increase = 0.1
  ) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc)
  )

This is the figure to report: the three-way ANOVA’s F, degrees of freedom, p-value and effect size up top, and the treatment differences that survive the decomposition marked right on the high-risk males.

Report the result

A publication-style write-up leads with the three-way interaction, then walks down the decomposition, each step with its F, degrees of freedom, p-value — for example:

A three-way ANOVA was conducted to determine the effects of gender, risk and treatment on migraine headache pain score. Residuals were normally distributed (Shapiro-Wilk, p > 0.05) and there was homogeneity of variances (Levene’s test, p > 0.05). There was a statistically significant three-way interaction between gender, risk and treatment, F(2, 60) = 7.41, p = 0.001, generalized η² = 0.20. Statistical significance was accepted at p < 0.025 for simple two-way interactions and simple simple main effects. There was a significant simple two-way interaction between risk and treatment for males, F(2, 60) = 5.25, p = 0.008, but not for females, F(2, 60) = 2.87, p = 0.065. For males, there was a significant simple simple main effect of treatment at high risk, F(2, 60) = 14.8, p < 0.0001, but not at low risk, F(2, 60) = 0.66, p = 0.521. All pairwise comparisons used a Bonferroni adjustment; for high-risk males, treatment X differed significantly from both Y and Z, while Y and Z did not differ.

NoteThe math (optional)

A three-way ANOVA partitions the total variation in the outcome into the parts each term explains plus what’s left over — the total (T), the three main effects (A, B, C), the three two-way interactions (AB, AC, BC), the three-way interaction (ABC) and the residual (R):

\[SS_T = SS_A + SS_B + SS_C + SS_{AB} + SS_{AC} + SS_{BC} + SS_{ABC} + SS_R\]

Each sum of squares becomes a mean square by dividing by its degrees of freedom, \(MS = SS / df\), and each term’s F-statistic is its mean square over the residual mean square:

\[F_{ABC} = \frac{MS_{ABC}}{MS_R}\]

So there is one F per term, each comparing the variation that term explains against the unexplained (within-cell) variation. The three-way term \(F_{ABC}\) asks whether the two-way interaction between two factors changes across the third — the question a three-way design exists to answer.

Try it live

Run a full three-way ANOVA on a different dataset — npk (a classic agricultural trial: crop yield under nitrogen N, phosphate P and potassium K, each present or absent). Edit the factors or the formula and re-run; the sandbox boots on first Run.

🟢 With an AI agent

Ask Prova “I have a numeric outcome and three grouping columns — run a three-way ANOVA on my data, check the assumptions, and if the three-way interaction is significant, walk me through decomposing it” — it answers with rstatix code you can run on your own data, then helps you read the seven terms and the simple two-way / simple simple follow-ups. The runtime is the judge. Ask Prova →

Common issues

Empty or tiny cells. Three factors multiply quickly — a 2 × 2 × 3 design already has 12 cells, and a missing factor combination makes the three-way interaction inestimable (or wildly unstable). Run table(mydata$A, mydata$B, mydata$C) first; every cell needs at least a few observations.

Reporting main effects when the three-way interaction is significant. When A:B:C is significant, the main-effect and two-way rows are misleading on their own — the whole point is that the lower-order pattern changes across the third factor. Interpret the three-way term first and decompose it (simple two-way → simple simple main → pairwise); don’t quote the top-level main effects in isolation.

Forgetting the pooled error term. In the decomposition, pass error = model (the full three-way lm()) to anova_test(). This borrows the three-way model’s residual error and degrees of freedom, which is more powerful than letting each subset estimate its own error — recommended when the equal-variance assumption holds.

Frequently asked questions

NoteWhat is a three-way ANOVA?

A three-way ANOVA tests the effect of three grouping factors — and all their interactions — on a single numeric outcome. It estimates three main effects, three two-way interactions and one three-way interaction. In rstatix it’s mydata %>% anova_test(outcome ~ A * B * C); in base R, summary(aov(outcome ~ A * B * C, data = mydata)).

NoteHow do I interpret a significant three-way interaction?

Interpret it before any lower-order term. A significant A:B:C means the two-way interaction between two factors itself changes across the levels of the third, so the main effects and two-way interactions are misleading alone. Decompose it in three steps: simple two-way interactions (the B:C interaction at each level of A), then simple simple main effects (the effect of one factor within each combination of the other two), then simple simple pairwise comparisons — each using the pooled error term and a Bonferroni-adjusted alpha.

NoteWhat is the difference between a two-way and a three-way ANOVA?

A two-way ANOVA crosses two factors (two main effects + one interaction). A three-way ANOVA crosses three factors, adding two more two-way interactions and a three-way interaction — seven terms in all. Use a three-way design when each observation is classified by three grouping variables and you suspect the effect of one depends on the combination of the other two.

NoteWhat is a simple two-way interaction and a simple simple main effect?

They are the follow-up analyses for a significant three-way interaction. A simple two-way interaction is the interaction between two factors computed at one level of the third (e.g. risk:treatment for males only). A simple simple main effect drills down once more: the effect of a single factor within each combination of the other two (e.g. the effect of treatment for males at high risk). Both are run with the pooled error = model term and a Bonferroni-adjusted significance level.

NoteWhat are the assumptions of a three-way ANOVA?

The same as the one- and two-way tests, applied to the cells: the observations are independent, there are no extreme outliers (identify_outliers()), the residuals are normally distributed (shapiro_test() on the model residuals + a Q-Q plot), and the cells have equal variance (levene_test(y ~ A * B * C)). Balanced cell sizes also keep the test well behaved; for an unbalanced design rstatix’s anova_test() already defaults to Type-II sums of squares.

Test your understanding

ImportantPractice

1. Run it. In the live cell below, fit a three-way ANOVA of migraine pain_score by gender, risk and treatment with all their interactions. Fill the blanks. Is the gender:risk:treatment term significant?
2. Conceptual. Your three-way interaction is significant. Why is it wrong to report only “the main effect of treatment”, and what is the correct sequence of follow-up analyses?

NoteHint

Fill the blanks so the formula is pain_score ~ gender * risk * treatment — the * between all three factors requests every main effect, every two-way interaction and the three-way interaction. Look at the gender:risk:treatment row’s p and p<.05 columns. For question 2, recall that a significant three-way interaction means a two-way pattern changes across the third factor.

NoteSolution

library(rstatix)
library(datarium)
data("headache", package = "datarium")
headache %>% anova_test(pain_score ~ gender * risk * treatment)
#> gender:risk:treatment is significant: F(2, 60) = 7.41, p = 0.001, ges = 0.20.

The gender:risk:treatment interaction is significant (p = 0.001).

For question 2: a significant three-way interaction means the effect of treatment depends on the combination of gender and risk, so a single “main effect of treatment” averages over a relationship that genuinely changes. Decompose it instead: simple two-way interactions (risk:treatment within each gender) → simple simple main effects (treatment within each gender × risk cell) → simple simple pairwise comparisons, each with the pooled error term and a Bonferroni adjustment.

TipWhich test when?

• One factor, between-subjects → one-way ANOVA.
• Two factors, between-subjects → two-way ANOVA.
• Three factors, between-subjects → three-way ANOVA (this lesson).
• Within-subjects (the same subjects measured under every condition) → repeated-measures ANOVA.
• Assumptions violated (non-normal residuals, unequal variances) → a robust ANOVA (WRS2) — see the normality test lesson for the full workflow.

Conclusion

You can now run and interpret a three-way ANOVA in R two ways: the tidy rstatix anova_test() (with a ges effect size per term) and the classic base-R aov() + summary(). The headline is the three-way interaction: when it’s significant — as it was for gender, risk and treatment on migraine pain — the factors don’t combine additively, so you interpret it first and decompose it through simple two-way interactions, simple simple main effects and simple simple pairwise comparisons. Check the assumptions — no extreme outliers, normal residuals (shapiro_test) and equal variance (levene_test) — and put the surviving differences on the plot with stat_pvalue_manual() so the figure reports itself.

Related lessons

• Two-way ANOVA in R — the two-factor version this extends, where the decomposition is just simple main effects. · One-way ANOVA in R — the single-factor starting point. · Repeated-measures ANOVA in R — the within-subjects design, when the same subjects are measured under every condition. · Normality test in R — the Shapiro-Wilk + Q-Q residual check the assumptions step relies on.
• Go further: The ANOVA family · Biostatistics

Was this page helpful?

Prove you can do it. Master the whole ANOVA in R series — track your path, build projects, and earn a certificate.

Start free →

Go Pro — unlimited Prova on your own data and a verifiable certificate that proves the skill.

from $15/mo billed yearly

Go Pro →

✓ You're Pro — keep going. The runtime is the judge.

Ready to level up?

Prove it — earn the ANOVA & Multi-Group Comparisons certificate.Take the assessment →Try this on your own data with Prova — the runtime is the judge.Open Prova →

Get new R & Python lessons by email

Practical, reproducible, no spam. Unsubscribe anytime.

Double opt-in. We never share your email.

Share this pageXLinkedInRedditHN

Note

Every result on this page was produced by the code shown, run at build time against a pinned R environment — edit any block and Run to reproduce it yourself.