DEV Community: peisen yan

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