ANOVA in R: A Complete Tutorial on Single-Factor and Two-Way Analysis of Variance

The Problem

"How does plant diversity differ between fertilized and control plots?" "Do different ecosystems respond to nitrogen addition in the same way?" These are everyday "compare several groups" questions in ecology. Running pairwise t-tests across many groups is tedious and inflates the probability of a Type I error (a false positive) through multiple comparisons. This is exactly what ANOVA (Analysis of Variance) is built for.

ANOVA uses a single F-test to ask whether all group means are equal, then post-hoc multiple comparisons locate which groups actually differ. Using a "nitrogen addition level × ecosystem type" experiment on plant diversity, this tutorial walks the full workflow: assumption checks, single-factor ANOVA, two-way ANOVA with an interaction effect, and simple effects via emmeans after a significant interaction.

How ANOVA Works

Single-factor ANOVA tests whether one categorical factor (e.g., nitrogen level) affects a continuous response (e.g., the Shannon diversity index). The core idea is splitting the total variation into:

• Between-group variation: how each group mean deviates from the grand mean — the treatment effect
• Within-group variation: how individuals within each group differ — the random error

The F-statistic is between-group mean square / within-group mean square. If the treatment effect far exceeds the error, the F-value is large and p < 0.05 rejects the null that all group means are equal.

Two-way ANOVA examines two factors at once (nitrogen level × ecosystem) and tests their interaction — whether the effect of one factor depends on the level of the other. For example, nitrogen addition might suppress grassland diversity more strongly than cropland; that difference in effect is the interaction.

ANOVA has three assumptions (critical — results are unreliable when they are violated):

Assumption	Test	Criterion
Independence	Guaranteed by design	Samples are independent
Normality	Shapiro-Wilk test (on residuals)	p > 0.05 (do not reject)
Homogeneity of variance	Levene test	p > 0.05 (do not reject)

These assumption tests treat the assumption as the null hypothesis, so p > 0.05 means "no evidence to reject" and the assumption is considered met; p < 0.05 means the assumption is violated and you need another method. How to proceed when assumptions fail depends on the design: for a single-factor design with unequal variances, prefer Welch ANOVA (oneway.test()); for a two-way design, do not simply fall back to Kruskal-Wallis (a rank-based test, not an "ANOVA substitute"), but choose a robust ANOVA, permutation test, or a generalized/mixed-effects model appropriate to the design.

Data

The tutorial uses simulated ecological data representing a "nitrogen addition level × ecosystem" factorial experiment on plant Shannon diversity.

Experimental design (3 × 3 = 9 groups, 20 quadrats each, 180 rows total):

Factor	Levels
Nlevel (nitrogen level)	Control / Low_N / High_N
Ecosystem	Grassland / Forest / Cropland

Response variable: Shannon diversity index (continuous). Simpson, Pielou, and species richness S are also available.

# Read data
dat <- read.csv("anova_diversity.csv", header = TRUE)
dat$Ecosystem <- factor(dat$Ecosystem)
dat$Nlevel <- factor(dat$Nlevel, levels = c("Control", "Low_N", "High_N"))
head(dat)

The data file anova_diversity.csv is provided with the article (or prepare your own with the same columns). The simulated pattern is that nitrogen addition lowers diversity, and grassland is the most sensitive to nitrogen — this is the setup for the interaction we find later.

Core Workflow: From Assumptions to Two-Way ANOVA

Step 1: Assumption checks

Normality — on residuals, not raw data

Many beginners run Shapiro-Wilk on the raw data, which is a common mistake. The ANOVA normality assumption applies to the model residuals — even if observed values per group depart from normality, as long as the residuals are approximately normal the parametric test is still valid.

library(car)
# Fit a single-factor model first, then test its residuals
fit1 <- aov(Shannon ~ Nlevel, data = dat)
shapiro.test(residuals(fit1))

Example output:

    Shapiro-Wilk normality test
W = 0.99373, p-value = 0.6408

p = 0.641 > 0.05: no evidence to reject residual normality, so the residuals are approximately normal (cross-check with the Q-Q plot below).

Homogeneity of variance — Levene test

Equal variances across groups are required. Prefer the Levene test (in the car package). car::leveneTest() is centered on the median by default (a Brown-Forsythe type test), which is more robust to non-normal data than the classic Bartlett test.

leveneTest(Shannon ~ Nlevel, data = dat)

Example output:

      Df F value  Pr(>F)
group  2  1.9516 0.1451
     177

p = 0.145 > 0.05: no evidence to reject homogeneity; variances are comparable across groups.

Q-Q plot to visualize normality

The points fall roughly along the diagonal with only slight deviation at the ends, visually supporting approximate residual normality.

library(ggplot2)
p_qq <- ggplot(data.frame(res = residuals(fit1)), aes(sample = res)) +
  stat_qq(color = "#5B8CC9", alpha = 0.7) +
  stat_qq_line(color = "#1F3B73", linewidth = 0.5) +
  labs(title = "Q-Q plot of ANOVA residuals", x = "Theoretical quantiles",
       y = "Sample quantiles") +
  theme_classic(base_size = 8, base_family = "Arial") +
  theme(axis.line = element_line(linewidth = 0.4),
        axis.ticks = element_line(linewidth = 0.4),
        panel.grid = element_blank())

Step 2: Single-factor ANOVA — the effect of nitrogen level

With assumptions met, fit a single-factor model with aov() and inspect the F-test with summary():

summary(fit1)

Example output:

             Df Sum Sq Mean Sq F value    Pr(>F)
Nlevel        2 26.860  13.430   125.2 < 2.2e-16 ***
Residuals   177 18.997   0.107

F(2, 177) = 125.2, p < 0.001: nitrogen level significantly affects Shannon diversity. But "overall significant" does not mean "every group differs" — you need post-hoc multiple comparisons to locate the differences.

Why Tukey HSD rather than pairwise t-tests? Three groups require 3 t-tests; judging each at α = 0.05 inflates the family-wise error rate to roughly 1-(1-0.05)³ ≈ 14%. Tukey HSD (Honestly Significant Difference) corrects this inflation using the studentized range distribution, keeping the family-wise error rate within α. Compared with Bonferroni correction (dividing α by the number of comparisons — conservative but low-powered), Tukey is the efficient standard choice when comparing all pairwise combinations. Switch to a Dunnett test only when the goal is "each treatment vs. control" rather than all pairs.

TukeyHSD(fit1)

Example output:

$Nlevel
                     diff       lwr       upr    p adj
Low_N-Control  -0.5491233 -0.690464 -0.407783 0.000000
High_N-Control -0.9419133 -1.083254 -0.800573 0.000000
High_N-Low_N   -0.3927900 -0.534130 -0.251450 0.000000

All three pairwise comparisons are significant (p adj < 0.001). diff is the mean difference; lwr/upr are the 95% confidence intervals. Nitrogen addition clearly reduces diversity: High_N averages 0.94 lower than Control, Low_N 0.55 lower.

In the single-factor boxplot, the three boxes use a blue gradient (light → dark) to convey increasing nitrogen level, with Tukey letters a/b/c above each box — different letters mean significant differences between groups. Median diversity falls from left to right (about 2.86 → 2.23 → 1.89), matching the statistics (means 2.81, 2.26, 1.87) and clearly showing the suppression of diversity by nitrogen.

library(multcompView)   # Tukey letters
# Compute Tukey letters (a/b/c)
tk <- TukeyHSD(fit1)$Nlevel
lv <- levels(dat$Nlevel)
pmat <- matrix(1, 3, 3, dimnames = list(lv, lv))
for (pr in rownames(tk)) {
  # Note: splitting group names by "-" assumes no hyphens in the names;
  # if your factor levels contain "-" (e.g. "Low-N"), use multcompLetters4()
  # or another more robust interface instead
  v <- strsplit(pr, "-")[[1]]
  pmat[v[2], v[1]] <- pmat[v[1], v[2]] <- tk[pr, "p adj"]
}
let1 <- multcompLetters(pmat)$Letters
lab1 <- data.frame(Nlevel = names(let1), letters = let1,
                   y = tapply(dat$Shannon, dat$Nlevel, max) + 0.1)

p_oneway <- ggplot(dat, aes(x = Nlevel, y = Shannon, fill = Nlevel)) +
  geom_boxplot(alpha = 0.85, outlier.shape = NA, linewidth = 0.4) +
  geom_jitter(width = 0.12, size = 0.6, alpha = 0.35, color = "grey35") +
  geom_text(data = lab1, aes(x = Nlevel, y = y, label = letters),
            size = 3, fontface = "bold") +
  scale_fill_manual(values = c(Control = "#D6E4F0", Low_N = "#5B8CC9",
                               High_N = "#1F3B73")) +
  labs(x = "Nitrogen addition level", y = "Shannon diversity index") +
  theme_classic(base_size = 8, base_family = "Arial") +
  theme(legend.position = "none",
        axis.line = element_line(linewidth = 0.4),
        axis.ticks = element_line(linewidth = 0.4),
        panel.grid = element_blank())

Step 3: Two-way ANOVA — nitrogen × ecosystem

Now consider both factors. First fit the main-effects model (additive, no interaction):

# Main-effects model: each factor's separate influence
fit2a <- aov(Shannon ~ Nlevel + Ecosystem, data = dat)
summary(fit2a)

Example output:

             Df F value    Pr(>F)
Nlevel        2 154.465 < 2.2e-16 ***
Ecosystem     2  21.693 3.83e-09 ***
Residuals   175

Both factors are highly significant. Next fit the interaction model to see whether the nitrogen effect changes across ecosystems:

# Interaction model: nitrogen and ecosystem interact
fit2b <- aov(Shannon ~ Nlevel * Ecosystem, data = dat)
summary(fit2b)

Example output:

                  Df  F value    Pr(>F)
Nlevel             2 164.3332 < 2.2e-16 ***
Ecosystem          2  23.0791 1.34e-09 ***
Nlevel:Ecosystem   4   3.7951  0.005537 **
Residuals        171

The interaction Nlevel:Ecosystem is significant (F = 3.80, p = 0.0055), meaning the nitrogen effect depends on ecosystem type — different ecosystems respond to nitrogen with different intensities.

Check assumptions on the model you actually use. Step 1 tested the single-factor model fit1 residuals, but the model you actually interpret is the two-way interaction model fit2b. The correct approach is to re-test the assumptions (normality + homogeneity) on the final model's residuals, rather than assuming single-factor assumptions carry over.

# Assumption checks on the final two-way model
shapiro.test(residuals(fit2b))
leveneTest(residuals(fit2b) ~ interaction(dat$Nlevel, dat$Ecosystem))

Example output:

    Shapiro-Wilk normality test
W = 0.99601, p-value = 0.919

Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   8  0.6552 0.7302
      171

fit2b residuals are normal (p = 0.919) and homogeneous (p = 0.730), so the two-way model's assumptions are satisfied. Whether running single- or two-factor, always run residual diagnostics on the final model — that is the habit of rigorous research.

Formula shorthand: A + B is the main-effects model; A * B expands automatically to A + B + A:B (with interaction). Keep the interaction term when it is significant; when it is not, do not delete it mechanically from the p-value alone — decide based on the study design, pre-specified hypotheses, and model-simplification principles, simplifying to an additive model only if appropriate.

When the interaction is significant, do not run an overall comparison of Nlevel (that would pool the three ecosystems together and mask the interaction). Instead, compute simple effects on the full model — fix Ecosystem and compare the three Nlevel groups within each panel. In research, use the emmeans package on the full model fit2b; it uses a single pooled residual variance and degrees of freedom (df = 171), which is more rigorous and accurate than subsetting the data into separate ANOVAs:

library(emmeans)
# Simple effects on the full model (fix Ecosystem, compare Nlevel within panels)
emm <- emmeans(fit2b, ~ Nlevel | Ecosystem)
pairs(emm, adjust = "tukey")   # pairwise Nlevel comparisons within each Ecosystem

# Convert emmeans corrected p-values to per-panel letters (pooled df=171)
ct <- as.data.frame(contrast(emm, method = "pairwise", adjust = "tukey"))
emmeans_letters <- function(ct_df, fac_levels) {
  pmat <- matrix(1, length(fac_levels), length(fac_levels),
                 dimnames = list(fac_levels, fac_levels))
  for (i in seq_len(nrow(ct_df))) {
    nm <- strsplit(as.character(ct_df$contrast[i]), " - ")[[1]]
    pmat[nm[2], nm[1]] <- pmat[nm[1], nm[2]] <- ct_df$p.value[i]
  }
  multcompLetters(pmat)$Letters
}
let2 <- lapply(levels(dat$Ecosystem), function(ec)
  emmeans_letters(ct[ct$Ecosystem == ec, ], levels(dat$Nlevel)))
names(let2) <- levels(dat$Ecosystem)
lab2 <- do.call(rbind, lapply(names(let2), function(ec) {
  data.frame(Ecosystem = ec, Nlevel = names(let2[[ec]]), letters = let2[[ec]],
             y = tapply(dat$Shannon[dat$Ecosystem == ec],
                        dat$Nlevel[dat$Ecosystem == ec], max) + 0.12)
}))

p_twoway <- ggplot(dat, aes(x = Nlevel, y = Shannon, fill = Nlevel)) +
  geom_boxplot(alpha = 0.85, outlier.shape = NA, linewidth = 0.4) +
  geom_jitter(width = 0.12, size = 0.6, alpha = 0.35, color = "grey35") +
  geom_text(data = lab2, aes(x = Nlevel, y = y, label = letters),
            size = 3, fontface = "bold") +
  facet_wrap(~ Ecosystem, ncol = 3) +
  scale_fill_manual(values = c(Control = "#D6E4F0", Low_N = "#5B8CC9",
                               High_N = "#1F3B73")) +
  labs(x = "Nitrogen addition level", y = "Shannon diversity index") +
  theme_classic(base_size = 8, base_family = "Arial") +
  theme(legend.position = "none",
        strip.text = element_text(face = "bold"),
        axis.line = element_line(linewidth = 0.4),
        axis.ticks = element_line(linewidth = 0.4),
        panel.grid = element_blank())

In the faceted two-way boxplot, nitrogen lowers diversity in all three panels (Cropland / Forest / Grassland). Each panel carries Tukey letters based on the emmeans simple effects (a/b/c, pooled df = 171). The Grassland panel starts highest (about 3.2) and drops the steepest, visually exposing the source of the interaction.

Advanced: Visualizing the Interaction Effect

Text alone is not intuitive for an interaction; an interaction plot (connecting group means) is a great aid. Note: non-parallel lines hint at an interaction, but whether one exists is decided by the interaction term's F-test in the ANOVA model (the Nlevel:Ecosystem term p = 0.0055 above); the interaction plot is a visual aid, not the formal test.

# Group means
mu <- aggregate(Shannon ~ Nlevel + Ecosystem, data = dat, FUN = mean)

# Interaction plot: whether lines are parallel reflects the interaction
p_interact <- ggplot(mu, aes(x = Nlevel, y = Shannon,
                             group = Ecosystem, color = Ecosystem)) +
  geom_point(size = 2.2) +
  geom_line(linewidth = 0.7) +
  scale_color_manual(values = c(Grassland = "#2C7A57", Forest = "#7A4D8E",
                                Cropland = "#B85C4A")) +
  labs(x = "Nitrogen addition level", y = "Mean Shannon diversity index") +
  theme_classic(base_size = 8, base_family = "Arial") +
  theme(axis.line = element_line(linewidth = 0.4),
        axis.ticks = element_line(linewidth = 0.4),
        panel.grid = element_blank())

All three lines decline but are not parallel: the Grassland line (green) is steepest (from about 3.10 down to about 1.87), while Forest (purple) and Cropland (red-brown) are nearly parallel and decline more gently (about 2.80→2.01 and 2.53→1.72 respectively). Grassland drops the most, which is why the interaction is significant. In a real study, such a pattern suggests grassland plant communities may be most sensitive to nitrogen addition, and should be interpreted in light of the experimental design.

Reporting simple effects with emmeans after a significant interaction

A significant interaction means the nitrogen effect depends on ecosystem, so reporting the Nlevel main effect alone is of limited value. The correct approach is to fix Ecosystem and compute simple effects with emmeans on the full model (i.e., compare the Nlevel groups within each Ecosystem):

library(emmeans)
fit2b <- aov(Shannon ~ Nlevel * Ecosystem, data = dat)
emm <- emmeans(fit2b, pairwise ~ Nlevel | Ecosystem, adjust = "tukey")
emm$contrasts

Example output (partial):

Ecosystem = Grassland:
 contrast          estimate    SE  df t.ratio p.value
 Control - Low_N      0.732 0.0904 171   8.097 <0.0001
 Control - High_N     1.225 0.0904 171  13.548 <0.0001
 Low_N - High_N       0.493 0.0904 171   5.450 <0.0001

emmeans uses a single pooled residual variance and degrees of freedom (df = 171) across all three ecosystems — the same error estimate everywhere — which is more rigorous than subsetting the data into separate ANOVAs. Subsetting uses each subset's own residual df, does not fully exploit the full-model information, and raises multiple-comparison problems for cross-ecosystem inference. For comparing groups after a significant interaction in a research paper, emmeans-style methods are the recommended approach. Report the pairwise contrast table (estimate, SE, df, corrected p-value) and the interaction plot as primary evidence; the letters on the boxplots are only a visualization aid.

Interpreting the Results

Three conclusions emerge from the analysis and figures:

1. Nitrogen addition broadly suppresses diversity. The single-factor ANOVA (F = 125.2, p < 0.001) shows nitrogen level significantly affects Shannon diversity (the Nlevel main effect in the two-way model is likewise significant, F = 164.3, p < 0.001, but must be read in the interaction context — see point 3). Tukey multiple comparisons show all three groups differ significantly; the letters on the boxplot are all different (a/b/c), and High_N is on average 0.94 lower than Control. The direction is consistent with the classic ecological finding that nitrogen deposition suppresses plant diversity.

2. Ecosystems differ fundamentally. Regardless of nitrogen level, Grassland has the highest diversity, Forest is intermediate, and Cropland the lowest; the Ecosystem main effect is highly significant (F = 21.7, p < 0.001). This gradient matches the general understanding that intensive agricultural management lowers plant diversity — long-term cultivation, weeding, and fertilization suppress the colonization of non-target species, homogenizing the community.

3. The interaction reveals differential responses. The significant interaction (F = 3.80, p = 0.0055) shows the nitrogen effect is not uniform: the interaction plot shows Grassland declining most steeply while Forest and Cropland are nearly parallel and decline more gently — the interaction comes mainly from grassland's steeper response. If a real study showed a similar interaction, it would suggest grassland ecosystems face a greater risk of diversity loss under nitrogen-deposition scenarios.

Reading Tukey letters: the letters above the boxplots come from the post-hoc grouping — groups sharing a letter do not differ significantly; groups with different letters do. All three groups here carry a, b, and c (all different), meaning every pair differs significantly.

Complete Runnable Code

Here is the full R script. Copy it into RStudio to run (the data file anova_diversity.csv is provided with the article):

# ============================================================
# Single-factor & two-way ANOVA: complete workflow
# ============================================================
# Install packages on first run:
# install.packages(c("car", "ggplot2", "multcompView", "emmeans"))

library(car)
library(ggplot2)
library(multcompView)   # Tukey letters
library(emmeans)        # simple effects (group comparison after a significant interaction)

# --- Read data ---
dat <- read.csv("anova_diversity.csv", header = TRUE)
dat$Ecosystem <- factor(dat$Ecosystem)
dat$Nlevel <- factor(dat$Nlevel, levels = c("Control", "Low_N", "High_N"))
head(dat)

# --- Shared theme (Nature-style minimal) ---
theme_pub <- function(base_size = 8, base_family = "Arial") {
  theme_classic(base_size = base_size, base_family = base_family) +
    theme(
      axis.line = element_line(linewidth = 0.4, colour = "black"),
      axis.ticks = element_line(linewidth = 0.4, colour = "black"),
      axis.title = element_text(size = base_size + 1),
      axis.text = element_text(size = base_size - 0.5),
      legend.title = element_text(size = base_size),
      legend.text = element_text(size = base_size - 1),
      strip.text = element_text(size = base_size, face = "bold"),
      legend.key.size = unit(3, "mm"),
      panel.grid = element_blank()
    )
}
# Color palette
cols_nlevel <- c(Control = "#D6E4F0", Low_N = "#5B8CC9", High_N = "#1F3B73")
cols_ecosystem <- c(Grassland = "#2C7A57", Forest = "#7A4D8E", Cropland = "#B85C4A")

# ============================================================
# Step 1: assumption checks
# ============================================================
fit1 <- aov(Shannon ~ Nlevel, data = dat)
shapiro.test(residuals(fit1))   # normality (on residuals)
leveneTest(Shannon ~ Nlevel, data = dat)   # homogeneity of variance

# ============================================================
# Step 2: single-factor ANOVA + Tukey multiple comparison
# ============================================================
summary(fit1)
TukeyHSD(fit1)

# ============================================================
# Step 3: two-way ANOVA (main-effects vs interaction model)
# ============================================================
fit2a <- aov(Shannon ~ Nlevel + Ecosystem, data = dat)
summary(fit2a)          # main-effects model

fit2b <- aov(Shannon ~ Nlevel * Ecosystem, data = dat)
summary(fit2b)          # interaction model

# Note: when the interaction is significant, do not run an overall Tukey
# comparison on Nlevel — instead use emmeans to look at simple effects
# within each Ecosystem (see the two-way visualization section)

# ============================================================
# Step 4: visualization
# ============================================================
# --- Tukey letters helper ---
tukey_letters <- function(model, fac, data) {
  tk <- TukeyHSD(model)[[fac]]
  lv <- levels(data[[fac]])
  pmat <- matrix(1, length(lv), length(lv), dimnames = list(lv, lv))
  for (pr in rownames(tk)) {
    v <- strsplit(pr, "-")[[1]]
    pmat[v[2], v[1]] <- pmat[v[1], v[2]] <- tk[pr, "p adj"]
  }
  multcompLetters(pmat)$Letters
}

# Single-factor boxplot (with Tukey letters)
let1 <- tukey_letters(fit1, "Nlevel", dat)
lab1 <- data.frame(Nlevel = names(let1), letters = let1,
                   y = tapply(dat$Shannon, dat$Nlevel, max) + 0.1)
p_oneway <- ggplot(dat, aes(x = Nlevel, y = Shannon, fill = Nlevel)) +
  geom_boxplot(alpha = 0.85, outlier.shape = NA, linewidth = 0.4) +
  geom_jitter(width = 0.12, size = 0.6, alpha = 0.35, color = "grey35") +
  geom_text(data = lab1, aes(x = Nlevel, y = y, label = letters),
            size = 3, fontface = "bold") +
  scale_fill_manual(values = cols_nlevel) +
  labs(x = "Nitrogen addition level", y = "Shannon diversity index") +
  theme_pub() +
  theme(legend.position = "none")

# Two-way faceted boxplot (per-panel letters, emmeans simple effects, pooled df=171)
emm <- emmeans(fit2b, ~ Nlevel | Ecosystem)
ct <- as.data.frame(contrast(emm, method = "pairwise", adjust = "tukey"))
emmeans_letters <- function(ct_df, fac_levels) {
  pmat <- matrix(1, length(fac_levels), length(fac_levels),
                 dimnames = list(fac_levels, fac_levels))
  for (i in seq_len(nrow(ct_df))) {
    nm <- strsplit(as.character(ct_df$contrast[i]), " - ")[[1]]
    pmat[nm[2], nm[1]] <- pmat[nm[1], nm[2]] <- ct_df$p.value[i]
  }
  multcompLetters(pmat)$Letters
}
let2 <- lapply(levels(dat$Ecosystem), function(ec)
  emmeans_letters(ct[ct$Ecosystem == ec, ], levels(dat$Nlevel)))
names(let2) <- levels(dat$Ecosystem)
lab2 <- do.call(rbind, lapply(names(let2), function(ec) {
  data.frame(Ecosystem = ec, Nlevel = names(let2[[ec]]), letters = let2[[ec]],
             y = tapply(dat$Shannon[dat$Ecosystem == ec],
                        dat$Nlevel[dat$Ecosystem == ec], max) + 0.12)
}))
p_twoway <- ggplot(dat, aes(x = Nlevel, y = Shannon, fill = Nlevel)) +
  geom_boxplot(alpha = 0.85, outlier.shape = NA, linewidth = 0.4) +
  geom_jitter(width = 0.12, size = 0.6, alpha = 0.35, color = "grey35") +
  geom_text(data = lab2, aes(x = Nlevel, y = y, label = letters),
            size = 3, fontface = "bold") +
  facet_wrap(~ Ecosystem, ncol = 3) +
  scale_fill_manual(values = cols_nlevel) +
  labs(x = "Nitrogen addition level", y = "Shannon diversity index") +
  theme_pub() +
  theme(legend.position = "none")

# Interaction plot
mu <- aggregate(Shannon ~ Nlevel + Ecosystem, data = dat, FUN = mean)
p_interact <- ggplot(mu, aes(x = Nlevel, y = Shannon,
                             group = Ecosystem, color = Ecosystem)) +
  geom_point(size = 2.2) +
  geom_line(linewidth = 0.7) +
  scale_color_manual(values = cols_ecosystem) +
  labs(x = "Nitrogen addition level", y = "Mean Shannon diversity index") +
  theme_pub()

# Save figures
ggsave("boxplot_oneway.png", p_oneway, width = 110, height = 100, units = "mm", dpi = 300)
ggsave("boxplot_twoway.png", p_twoway, width = 183, height = 90, units = "mm", dpi = 300)
ggsave("interaction_plot.png", p_interact, width = 120, height = 100, units = "mm", dpi = 300)

Adapting to your own data:

1. Prepare a CSV: factor columns (categorical) first, the continuous response as the last column
2. Convert factor columns with factor() (e.g., dat$Group <- factor(dat$Group))
3. Replace the filename in read.csv() and the variable names in the formulas (e.g., aov(value ~ Group))
4. For a single-factor design with unequal variances, switch to oneway.test(Shannon ~ Nlevel, data = dat) (Welch ANOVA); for a two-way design with unequal variances, choose a robust method appropriate to the design

Practical Tips

aov formula cheat sheet:

Formula	Meaning
aov(y ~ A)	Single-factor ANOVA
aov(y ~ A + B)	Two-way main-effects model
aov(y ~ A * B)	Two-way interaction model (equivalent to A + B + A:B)

Common gotchas:

• Test residuals, not raw data. Running Shapiro-Wilk on raw data often falsely reports non-normality; the ANOVA normality assumption applies to the model residuals, not the observations themselves
• Test assumptions on the final model. Single-factor and two-factor are different models; each needs its own residual diagnostics — there is no "test once, done forever"
• Use Levene, not Bartlett, for homogeneity. car::leveneTest() is median-centered by default (Brown-Forsythe type), more robust to non-normal data than Bartlett
• "p > 0.05" is "no evidence to reject", not "proven". The null is that the assumption holds; with large samples Shapiro can be sensitive to tiny deviations, so combine it with a Q-Q plot
• Locate differences after a significant result. ANOVA significant only says "there is a difference"; which pairs differ requires post-hoc tests such as TukeyHSD
• Tukey letters are one common annotation. Use multcompView::multcompLetters() to turn post-hoc results into a/b/c letters on boxplots; different letters mean significant differences. Common in ecology/agriculture papers (not a journal-mandated standard; some journals prefer reporting adjusted p-values directly)
• Use a gradient for dose-response factors. When treatment levels have an increasing meaning (e.g., 0/low/high nitrogen), shades of one color from light to dark convey "increasing dose" better than different hues
• Be cautious about main effects when the interaction is significant. A significant interaction means one factor's effect depends on the other, so read main effects in the interaction context; after a significant interaction, use emmeans for simple effects
• Watch for pseudoreplication. If several quadrats come from the same site/plot (nested or repeated-measures structure), they are not fully independent replicates and ordinary ANOVA overestimates sample size. Use a mixed-effects model instead, e.g., lmer(Shannon ~ Nlevel * Ecosystem + (1 | Site), data = dat)
• A non-significant interaction does not force deletion. Whether to simplify to an additive model should be decided from the study design and pre-specified hypotheses, not mechanically from the p-value
• Handle unequal variance by design. Single-factor: Welch ANOVA (oneway.test()); two-factor: do not simply use Kruskal-Wallis, but choose a robust/permutation/generalized/mixed model appropriate to the design

---

Follow me on Dev.to for more R tutorials on ecological data analysis.