DEV Community: peisen yan

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

peisen yan — Sat, 08 Aug 2026 13:50:31 +0000

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:

Prepare a CSV: factor columns (categorical) first, the continuous response as the last column
Convert factor columns with factor() (e.g., dat$Group <- factor(dat$Group))
Replace the filename in read.csv() and the variable names in the formulas (e.g., aov(value ~ Group))
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.

Cluster Analysis in R: Hierarchical Clustering & K-means

peisen yan — Wed, 10 Jun 2026 11:53:37 +0000

Cluster Analysis in R: Hierarchical Clustering & K-means

Samples with similar species composition often reflect similar environmental conditions. Cluster analysis quantifies this similarity and helps ecologists understand community structure.

This tutorial demonstrates a complete workflow from distance calculation to hierarchical clustering and K-means clustering using simulated community data (30 plots × 20 species).

Data Preparation

The dataset contains 30 plots and 20 species, designed with 3 community types (10 plots each):

Type	Plots	Dominant Species	Characteristic
Type A	Plot01-10	Sp01-Sp08	Higher abundance of species 1-8
Type B	Plot11-20	Sp09-Sp15	Higher abundance of species 9-15
Type C	Plot21-30	Sp16-Sp20	Higher abundance of species 16-20

# Read data
df <- read.csv("community_data.csv", header = TRUE, row.names = 1)
head(df[, 1:6])

Hierarchical Clustering

Step 1: Calculate Distance Matrix

For species abundance data, Jaccard distance or Bray-Curtis distance are most commonly used.

library(vegan)
dist_jac <- vegdist(df, method = "jaccard")

Step 2: UPGMA Clustering

hc_upgma <- hclust(dist_jac, method = "average")

method = "average" is UPGMA (Unweighted Pair Group Method with Arithmetic Mean). Other options:

Method	Description
`"average"`	UPGMA, most common, balanced branches
`"ward.D2"`	Ward's minimum variance, compact clusters
`"complete"`	Complete linkage, based on farthest points
`"single"`	Single linkage, based on nearest points

Step 3: Plot Dendrogram

# Vertical dendrogram
par(mar = c(4, 4, 3, 8))
plot(hc_upgma, hang = -1, main = "Hierarchical Clustering (UPGMA, Jaccard)",
     xlab = "Sample", ylab = "Jaccard Distance", sub = "", cex = 0.7)
rect.hclust(hc_upgma, k = 3, border = c("#E64B35", "#4DBBD5", "#00A087"))

Parameter notes:

hang = -1: Aligns all branch tips to the same horizontal line
rect.hclust(k = 3): Draws rectangles around 3 clusters
cex = 0.7: Adjusts label font size

# Horizontal dendrogram (better for many samples)
plot(as.dendrogram(hc_upgma), horiz = TRUE,
     main = "Hierarchical Clustering (UPGMA, Jaccard)")

Step 4: Cut the Tree

hc_clusters <- cutree(hc_upgma, k = 3)
print(table(hc_clusters))

cutree() cuts the dendrogram at a specified height and returns cluster membership.

K-means Clustering

Hierarchical clustering doesn't require pre-specifying the number of clusters, but when you have a prior estimate, K-means is faster.

Step 1: Dimensionality Reduction

K-means performs poorly on high-dimensional species data. Reduce dimensions first:

dist_bc <- vegdist(df, method = "bray")
pcoa <- cmdscale(dist_bc, k = 4, eig = TRUE)

Step 2: K-means

set.seed(123)
km <- kmeans(pcoa$points, centers = 3, nstart = 50)
km$cluster

Parameters:

centers = 3: Number of clusters
nstart = 50: 50 random initializations, takes the best result
set.seed(123): Ensures reproducibility

Step 3: Visualization

library(factoextra)
fviz_cluster(km, data = pcoa$points,
             ellipse.type = "norm",
             palette = c("#E64B35", "#4DBBD5", "#00A087"),
             ggtheme = theme_bw(),
             main = "K-means Clustering (k = 3)",
             labelsize = 0) +
  labs(color = "Cluster", shape = "Cluster")

fviz_cluster() projects data onto the first two principal coordinates, with ellipse.type = "norm" adding 95% confidence ellipses.

Results

The dendrogram and K-means plot clearly show:

Three distinct community types — The three main branches in the dendrogram correspond to the three simulated community types.
UPGMA and K-means agree — Both methods produce highly consistent classification results, approaching the problem from different angles.
Applications:
- Classifying community types based on species composition
- Validating predefined habitat/treatment groupings
- Providing categorical variables for downstream ordination (NMDS, RDA)

Complete Runable Script

# ============================================================
# Complete workflow: Hierarchical + K-means clustering
# ============================================================
# First run: install.packages(c("vegan", "factoextra"))

library(vegan)
library(factoextra)

df <- read.csv("community_data.csv", header = TRUE, row.names = 1)

# Hierarchical clustering
dist_jac <- vegdist(df, method = "jaccard")
hc_upgma <- hclust(dist_jac, method = "average")

pdf("dendrogram.pdf", width = 12, height = 7)
par(mar = c(4, 4, 3, 8))
plot(hc_upgma, hang = -1, main = "Hierarchical Clustering (UPGMA, Jaccard)",
     xlab = "Sample", ylab = "Jaccard Distance", sub = "", cex = 0.7)
rect.hclust(hc_upgma, k = 3, border = c("#E64B35", "#4DBBD5", "#00A087"))
dev.off()

pdf("dendrogram_horizontal.pdf", width = 8, height = 10)
par(mar = c(4, 4, 2, 2))
plot(as.dendrogram(hc_upgma), horiz = TRUE,
     main = "Hierarchical Clustering (UPGMA, Jaccard)")
dev.off()

hc_clusters <- cutree(hc_upgma, k = 3)

# K-means
dist_bc <- vegdist(df, method = "bray")
pcoa <- cmdscale(dist_bc, k = 4, eig = TRUE)
set.seed(123)
km <- kmeans(pcoa$points, centers = 3, nstart = 50)

pdf("kmeans_clusters.pdf", width = 8, height = 6)
fviz_cluster(km, data = pcoa$points,
             ellipse.type = "norm",
             palette = c("#E64B35", "#4DBBD5", "#00A087"),
             ggtheme = theme_bw(),
             main = "K-means Clustering (k = 3)",
             labelsize = 0) +
  labs(color = "Cluster", shape = "Cluster")
dev.off()

# Compare results
data.frame(Sample = rownames(df),
           Hierarchical = paste0("G", hc_clusters),
           Kmeans = paste0("G", km$cluster))

To use with your own data:

Prepare a CSV file (rows = samples, columns = species, values = abundance)
Replace community_data.csv with your file name
Adjust k = 3 to your expected number of clusters

Quick Tips

Use Jaccard distance for presence/absence or abundance data with many zeros
Use Bray-Curtis for raw abundance data
Always set nstart >= 50 for K-means to avoid local optima
Cross-validate clustering results with two different methods
Avoid pure numbers as sample names (use "Plot1" instead of "1")

Follow me on Dev.to for more R tutorials.

Correlation Analysis in R: From Pearson to Spearman with corrplot

peisen yan — Sun, 31 May 2026 15:34:17 +0000

When studying how environmental factors shape biodiversity, correlation analysis is usually the first stop. But which method should you use — Pearson, Spearman, or Kendall? Pick wrong and you might draw misleading conclusions.

This tutorial walks through the complete workflow using real-world ecological data (50 sites × 12 variables), covering method selection, correlation computation, significance testing, and corrplot visualization at publication quality.

Quick Reference: Which Method to Use?

Method	Normality Required	Relationship	Outlier Sensitivity	Best For
Pearson	Yes	Linear	Sensitive	Normally distributed continuous data
Spearman	No	Monotonic	Robust	Default for ecological data ✅
Kendall	No	Monotonic	Robust	Small samples (n < 30) with many ties

Rule of thumb: Run Shapiro-Wilk first. If all variables are normal → Pearson. Otherwise → Spearman (most common in ecology).

Data Overview

The dataset contains 50 sampling sites with 12 variables:

Environmental factors (6): pH, TN (total nitrogen), TP (total phosphorus), SOM (soil organic matter), MAP (mean annual precipitation), MAT (mean annual temperature)

Diversity indices (6): PD (phylogenetic diversity), NRI (net relatedness index), NTI (nearest taxon index), Shannon-Wiener index, Simpson index, S (species richness)

# Load data
df <- read.csv("alpha_diversity.csv", header = TRUE)
head(df)

The data file (alpha_diversity.csv) contains 50 rows and 12 numeric columns — one row per sampling site, one column per variable.

Step 1: Normality Test

Before choosing a correlation method, test whether your data follows a normal distribution:

# Shapiro-Wilk normality test for all variables
p_values <- sapply(df, function(x) shapiro.test(x)$p.value)
norm_table <- data.frame(
  Variable = names(df),
  p_value = round(p_values, 4),
  Normal = ifelse(p_values > 0.05, "Yes", "No")
)
print(norm_table)

In our dataset, most environmental factors and diversity indices have p < 0.05 (non-normal). This means Spearman is the appropriate choice.

Step 2: Calculate Correlation Coefficients

R's cor() function computes all three methods in one line:

# All three methods at once
pearson_r  <- cor(df, method = "pearson")
spearman_r <- cor(df, method = "spearman")
kendall_r  <- cor(df, method = "kendall")

The sign and relative magnitude are broadly consistent across methods. When a relationship is significant under Spearman but not Pearson, it often indicates a non-linear monotonic relationship — exactly where Spearman excels.

Step 3: Get Coefficients and p-values Together

cor() only returns the correlation matrix. For significance tests, use rcorr() from the Hmisc package:

library(Hmisc)
res <- rcorr(as.matrix(df), type = "spearman")

res$r  # Correlation coefficients
res$P  # p-values

Example interpretation (pH vs. Shannon):

Correlation coefficient r = 0.54
p < 0.001 (highly significant)
Interpretation: Higher pH is associated with higher Shannon diversity — soil acidity is an important driver of community diversity.

# Export results
write.csv(round(res$r, 3), "spearman_corr_coef.csv")
write.csv(round(res$P, 4), "spearman_p_values.csv")

Step 4: corrplot Visualization

Let's build up from a basic plot to a publication-ready figure. All examples use corrplot.

4.1 Basic corrplot — One Line

library(corrplot)
corrplot(res$r)

This gives you a colored circle matrix out of the box — blue for positive, red for negative, color intensity reflects strength. It's useful for a quick look but not publication-ready (no labels, no sorting).

4.2 Custom Colors + Upper/Lower Triangular Mix

Three improvements:

order = "AOE" — auto-cluster by angular order of eigenvectors
Upper triangle: colored circles; Lower triangle: numbers
Custom blue-white-red color palette

corrplot(res$r, order = "AOE", type = "upper", tl.pos = "d",
         tl.col = "black", tl.cex = 0.8,
         col = colorRampPalette(c("#0072B5", "white", "#BC3C29"))(200))
corrplot(res$r, add = TRUE, type = "lower", method = "number",
         order = "AOE", diag = FALSE, tl.pos = "n", cl.pos = "n",
         col = "black", number.cex = 0.7)

Key parameters:

order = "AOE": Reorders variables to reveal clustered structure
add = TRUE: Overlays the second layer (numbers) on top of the first (circles)
tl.pos = "d": Places variable labels on the diagonal to avoid overlap

4.3 Publication Standard: Significance Star Annotation ⭐

This is the format most commonly seen in peer-reviewed papers — all correlations shown, significant ones marked with stars (* p < 0.05, ** p < 0.01, *** p < 0.001):

p_mat <- res$P
diag(p_mat) <- 1
order_aoe <- corrMatOrder(res$r, order = "AOE")

corrplot(res$r[order_aoe, order_aoe], type = "upper", tl.pos = "d",
         tl.col = "black", tl.cex = 0.8,
         col = colorRampPalette(c("#0072B5", "white", "#BC3C29"))(200),
         p.mat = p_mat[order_aoe, order_aoe],
         sig.level = c(0.001, 0.01, 0.05),
         insig = "label_sig",
         pch.cex = 1.0, pch.col = "black")
corrplot(res$r[order_aoe, order_aoe], add = TRUE, type = "lower",
         method = "number", diag = FALSE, tl.pos = "n", cl.pos = "n",
         col = "black", number.cex = 0.7)

Critical details:

sig.level must be in ascending order (0.001, 0.01, 0.05)
insig = "label_sig" annotates significant correlations with stars
p.mat must use the same AOE ordering as the correlation matrix to keep upper and lower triangles aligned
Pre-compute order_aoe with corrMatOrder() and apply it to both res$r and p_mat

Step 5: Reading the Results

The corrplot reveals clear patterns in environment–diversity relationships:

Environment vs. Diversity:

pH positively correlated with Shannon (r = 0.54, p < 0.001), Simpson (r = 0.46, p < 0.001), and S (r = 0.44, p < 0.01) — near-neutral to alkaline soils support higher diversity
TN shows consistent positive correlations with diversity indices (r ≈ 0.43–0.53) — nitrogen availability is another key driver
SOM only weakly correlates with PD (r = 0.32, p < 0.05)

Among diversity indices:

Shannon and Simpson are strongly correlated (r = 0.85) — they capture overlapping information
Shannon and S also highly correlated (r = 0.74) — Shannon is inherently influenced by richness
NRI and NTI moderately correlated (r = 0.53) — both reflect phylogenetic structure

Method note: Comparing Pearson and Spearman results, most significant relationships agree in direction, but Spearman's rank correlation is more robust when normality assumptions are violated. Spearman is recommended as the default in ecological studies.

Complete Runnable Script

# ============================================================
# Complete Correlation Analysis: Pearson / Spearman / Kendall
# ============================================================
# First run: install.packages(c("Hmisc", "corrplot"))

library(Hmisc)
library(corrplot)

# --- Load data ---
df <- read.csv("alpha_diversity.csv", header = TRUE)

# --- Normality test ---
p_vals <- sapply(df, function(x) shapiro.test(x)$p.value)
data.frame(Variable = names(df), p_value = round(p_vals, 4),
           Normal = ifelse(p_vals > 0.05, "Yes", "No"))

# --- Three correlation methods ---
cor(df, method = "pearson")
cor(df, method = "spearman")
cor(df, method = "kendall")

# --- Spearman with p-values ---
res <- rcorr(as.matrix(df), type = "spearman")
res$r
res$P

# --- Visualization ---

# 1) Basic
pdf("corrplot_basic.pdf", width = 10, height = 8)
corrplot(res$r)
dev.off()

# 2) Custom colors + upper/lower mix
pdf("corrplot_color.pdf", width = 10, height = 8)
corrplot(res$r, order = "AOE", type = "upper", tl.pos = "d",
         tl.col = "black", tl.cex = 0.8,
         col = colorRampPalette(c("#0072B5", "white", "#BC3C29"))(200))
corrplot(res$r, add = TRUE, type = "lower", method = "number",
         order = "AOE", diag = FALSE, tl.pos = "n", cl.pos = "n",
         col = "black", number.cex = 0.7)
dev.off()

# 3) Publication standard with significance stars
p_mat <- res$P
diag(p_mat) <- 1
order_aoe <- corrMatOrder(res$r, order = "AOE")
pdf("corrplot_significance.pdf", width = 10, height = 8)
corrplot(res$r[order_aoe, order_aoe], type = "upper", tl.pos = "d",
         tl.col = "black", tl.cex = 0.8,
         col = colorRampPalette(c("#0072B5", "white", "#BC3C29"))(200),
         p.mat = p_mat[order_aoe, order_aoe],
         sig.level = c(0.001, 0.01, 0.05),
         insig = "label_sig",
         pch.cex = 1.0, pch.col = "black")
corrplot(res$r[order_aoe, order_aoe], add = TRUE, type = "lower",
         method = "number", diag = FALSE, tl.pos = "n", cl.pos = "n",
         col = "black", number.cex = 0.7)
dev.off()

To use with your own data: Prepare a CSV file where each column is a numeric variable and each row is an observation. Replace "alpha_diversity.csv" in the read.csv() call with your file name.

Key Takeaways

Check normality first — Shapiro-Wilk test tells you which correlation method to use
Spearman is the safe default for ecological and environmental data
corrplot with order = "AOE" reveals clustered structure automatically
Significance star annotations make your figure publication-ready
Correlation ≠ causation — even r = 0.9 doesn't mean A causes B
Multiple testing — with 12 variables (66 pairs), some significant results may be random. Consider Bonferroni or FDR correction for strict inference

Need help with R data analysis? I provide data analysis, visualization, and publication-ready figure services.

Follow me on Dev.to for more R tutorials. GitHub and project links available on request.

Let me know in the comments if you have questions or want the example dataset!

Scatter Plots with Linear Regression in R (with Multi-Panel Layouts)

peisen yan — Wed, 20 May 2026 09:44:22 +0000

The Problem

You have environmental data (pH, organic matter, nutrients) and biodiversity data (Shannon index). You want to know: which environmental factors significantly affect diversity?

The quickest way to answer this is scatter plots + linear regression.

Here's a complete, runnable R workflow using real soil microbiome data — from a simple lm() call to a multi-panel publication-ready figure. Everything uses base R graphics, no ggplot2 required.

Data

This dataset has 24 soil samples with 6 environmental factors and 2 diversity indices:

Variable	Description	Unit
pH	Soil acidity	—
AK	Available potassium	mg/kg
AP	Available phosphorus	mg/kg
AN	Available nitrogen	mg/kg
OM	Organic matter	g/kg
W	Water content	%
Shannon	Species diversity	—
PD	Phylogenetic diversity	—

# Read and inspect
df <- read.csv("sandian.csv", header = TRUE)
str(df)
head(df)

Single Factor: pH → Shannon

x <- df$pH
y <- df$Shannon
fit <- lm(y ~ x)
summary(fit)

Multiple R-squared:  0.0004, p-value: 0.927

R² = 0.0004, p = 0.927 — pH has no significant linear relationship with Shannon diversity in this dataset. That's OK. Non-significant results are normal in real research.

Plot it:

par(mar = c(5, 5, 3, 1))
plot(x, y,
     xlab = "pH", ylab = "Shannon Index",
     pch = 16, cex = 1.5, col = "#00AFBB",
     cex.lab = 1.8, cex.axis = 1.3, las = 1)
abline(fit, lwd = 3, col = "#00AFBB")

Multi-Panel: All 6 Factors at Once

Doing one plot at a time is inefficient. Use layout() + for loop:

x_vars <- c("pH", "AK", "AP", "AN", "OM", "W")
x_labels <- c("pH", "AK (mg/kg)", "AP (mg/kg)",
              "AN (mg/kg)", "SOM (g/kg)", "W (%)")

layout(matrix(1:6, 2, 3, byrow = TRUE))
par(mar = c(4.5, 4.5, 3, 1))

for (i in seq_along(x_vars)) {
  fit <- lm(df$Shannon ~ df[[x_vars[i]]])
  s <- summary(fit)

  plot(df[[x_vars[i]]], df$Shannon,
       xlab = x_labels[i], ylab = "Shannon",
       pch = 16, cex = 1.3, col = "#00AFBB",
       cex.lab = 1.5, cex.axis = 1.2, las = 1)
  abline(fit, lwd = 3, col = "#00AFBB")

  # Annotate R² and p-value
  r2 <- round(s$r.squared, 3)
  p  <- s$coefficients[2, 4]
  p_label <- ifelse(p < 0.001, "p<0.001", paste0("p=", round(p, 3)))

  x_pos <- par("usr")[2] - diff(par("usr")[1:2]) * 0.25
  y_pos <- par("usr")[4] - diff(par("usr")[3:4]) * 0.1
  text(x_pos, y_pos, labels = bquote(R^2 == .(r2) ~ " " ~ .(p_label)),
       cex = 1.1, adj = c(1, 1))
}

Save at 300 dpi for publication:

png("scatter_6panel.png", width = 3000, height = 2000, res = 300)
# ... plotting code ...
dev.off()

Results at a Glance

Factor	R²	p-value	Significant?
pH	0.000	0.927	❌
AK	0.088	0.160	❌
AP	0.012	0.615	❌
AN	0.220	0.021	✅
OM	0.540	<0.001	✅✅
W	0.000	0.956	❌

Organic matter (OM) is the strongest predictor — it explains 54% of the variance in Shannon diversity.

Complete Script (Copy & Run)

# ================================================================
# Scatter Plots + Linear Regression: Environmental Factors vs. Diversity
# ================================================================

plot_one <- function(x, y, xlab, ylab) {
  fit <- lm(y ~ x)
  s <- summary(fit)
  r2 <- round(s$r.squared, 3)
  p  <- s$coefficients[2, 4]
  p_label <- ifelse(p < 0.001, "p<0.001", paste0("p=", round(p, 3)))

  plot(x, y, xlab = xlab, ylab = ylab,
       pch = 16, cex = 1.3, col = "#00AFBB",
       cex.lab = 1.5, cex.axis = 1.2, las = 1)
  abline(fit, lwd = 3, col = "#00AFBB")

  x_pos <- par("usr")[2] - diff(par("usr")[1:2]) * 0.25
  y_pos <- par("usr")[4] - diff(par("usr")[3:4]) * 0.1
  text(x_pos, y_pos, labels = bquote(R^2 == .(r2) ~ " " ~ .(p_label)),
       cex = 1.1, adj = c(1, 1))
}

df <- read.csv("sandian.csv", header = TRUE)

x_vars   <- c("pH", "AK", "AP", "AN", "OM", "W")
x_labels <- c("pH", "AK (mg/kg)", "AP (mg/kg)",
              "AN (mg/kg)", "SOM (g/kg)", "W (%)")

png("scatter_6panel.png", width = 3000, height = 2000, res = 300)
layout(matrix(1:6, 2, 3, byrow = TRUE))
par(mar = c(4.5, 4.5, 3, 1))

for (i in seq_along(x_vars)) {
  plot_one(df[[x_vars[i]]], df$Shannon, x_labels[i], "Shannon")
}
dev.off()

To adapt for your data: change x_vars to your column names, x_labels to your axis labels, and df$Shannon to your diversity index.

Parameter Cheat Sheet

param	what it does	common values
`pch`	point shape	16 (filled), 1 (empty), 17 (triangle)
`cex`	point/text size	1.3–1.5 (moderate)
`col`	color	`#00AFBB` (teal), `#E69F00` (orange)
`lty`	line type	1 (solid), 2 (dashed)
`lwd`	line width	2–3 (publication)
`res`	export dpi	300 (journal)

When R² Is Not Significant

Normal. Don't only show significant results.
Try non-linear: poly(x, 2) for quadratic terms.
Small n (< 30) : may lack statistical power, not ecological meaning.

*Follow me on Dev.to for more R tutorials.

Tags: r tutorial datascience ecology