Analyzing Interaction Effects and Multi-Factor Grouping
A Two-Way ANOVA examines the effect of two independent categorical variables on a continuous dependent variable. Crucially, it also tests for interaction effects—determining if the effect of one factor depends on the level of the other.
We use the showtext package to bring professional Google Fonts into our R visualizations.
library(tidyverse)
library(ggthemes)
library(multcompView)
library(stats)
library(showtext)
# Enable Poppins font for a modern look
font_add_google("Poppins", "poppins")
showtext_auto()
# Set global theme with custom typography
theme_set(theme_test(base_size = 15) +
theme(text = element_text(family = "poppins")))We use the ToothGrowth dataset. In this analysis, we examine how Supplement Type (supp) and Dose (dose) together affect tooth length.
# Load dataset
df <- ToothGrowth
# Ensure independent variables are factors
df$dose <- as.factor(df$dose)
df$supp <- as.factor(df$supp)We use the asterisk (*) syntax to calculate the main effects of both factors and their interaction.
# len ~ supp * dose means Main Effect 1 + Main Effect 2 + Interaction
anova_result <- aov(len ~ supp * dose, data = df)
summary(anova_result) 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 ' ' 1If the interaction or main effects are significant, we perform a Tukey HSD test to compare all possible group combinations. We then use Compact Letter Display (CLD) to simplify the visual communication of these differences.
# Perform Multiple Mean Comparison
tukey_result <- TukeyHSD(anova_result)
# Extract significance lettering
group_lettering <- multcompLetters4(anova_result, tukey_result)
# Convert to dataframe for merging
group_lettering_df <- data.frame(
Letters = group_lettering$`supp:dose`$Letters
) %>%
rownames_to_column("group")We calculate the mean and standard deviation for each factor combination to prepare our bar plot.
# Calculate summary stats
df_summary <- df %>%
group_by(supp, dose) %>%
summarise(len_mean = mean(len),
sd = sd(len),
.groups = "drop") %>%
mutate(group = paste0(supp, ":", dose)) %>%
left_join(group_lettering_df, by = "group")
# View summary table
knitr::kable(df_summary)| supp | dose | len_mean | sd | group | Letters |
|---|---|---|---|---|---|
| OJ | 0.5 | 13.23 | 4.459708 | OJ:0.5 | b |
| OJ | 1 | 22.70 | 3.910953 | OJ:1 | a |
| OJ | 2 | 26.06 | 2.655058 | OJ:2 | a |
| VC | 0.5 | 7.98 | 2.746634 | VC:0.5 | c |
| VC | 1 | 16.77 | 2.515309 | VC:1 | b |
| VC | 2 | 26.14 | 4.797731 | VC:2 | a |
We create a “Dodged” bar plot where the grouping is determined by the supplement type across different doses.
df_summary %>%
ggplot(aes(x = dose, y = len_mean, fill = supp)) +
geom_bar(position = position_dodge(0.9),
stat = "identity",
alpha = 0.8) +
geom_errorbar(aes(ymin = len_mean - sd,
ymax = len_mean + sd),
width = 0.25,
position = position_dodge(0.9)) +
geom_text(aes(label = Letters,
y = len_mean + sd),
vjust = -0.9,
color = "gray25",
position = position_dodge(0.9),
size = 5,
show.legend = FALSE) +
scale_fill_manual(values = c("OJ" = "#8c52ff",
"VC" = "#cb6ce6")) +
labs(x = "Dose (mg/day)",
y = "Tooth Length (mm)",
title = "Two-Way ANOVA Result Visualization",
subtitle = "Effect of Dose and Supplement Type on Tooth Growth",
caption = "Data Source: ToothGrowth | Analyzed by Abdullah Al Shamim",
fill = "Supplement") +
ylim(0, 40) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(0.15, 0.85))
DV ~ Factor1 * Factor2 (includes interaction).Factor1:Factor2 is \(< 0.05\), the factors depend on each other.supp:dose section of the Tukey output for interaction significance.Excellent Progress! You have now moved from single-factor analysis to complex interaction modeling.