Module 8: Comparing Environmental Groups: Two-Sample Tests, One-Way ANOVA & Post-Hoc Analysis
Module 8: Comparing Environmental Groups: Two-Sample Tests, One-Way ANOVA & Post-Hoc Analysis
Executing independent and paired Student's and Welch's t-tests, Mann-Whitney U test, One-Way Analysis of Variance (ANOVA), F-ratio interpretation, Tukey HSD post-hoc testing, and Kruskal-Wallis with Dunn test.
8. Comparing Environmental Groups: Two-Sample Tests, ANOVA & Post-Hoc Analysis
A central objective in environmental science is comparing parameters between experimental treatments or monitoring zones. This module covers executing two-sample tests and multi-group ANOVA in R with post-hoc pairwise comparisons.
Student's $t$-test vs. Welch's $t$-test
In classic Student's $t$-test, equal group variances are assumed ($\sigma_1^2 = \sigma_2^2$). In real environmental data, variances are rarely equal. Welch's $t$-test adjusts degrees of freedom according to the Welch-Satterthwaite equation, protecting against Type I error inflation. In R, t.test() runs Welch's test by default (var.equal = FALSE).
One-Way ANOVA & Post-Hoc Multiple Comparisons
When comparing $\ge 3$ groups, running multiple pairwise $t$-tests causes severe Family-Wise Error Rate (FWER) inflation. For 4 groups (6 pairwise tests), the false positive risk balloons from $5\%$ to $1 - (1 - 0.05)^6 \approx 26.5\%$!
One-Way ANOVA computes an overall omnibus $F$-statistic partitioning total sum of squares into Between-Group Variance vs Within-Group (Residual) Variance: $$F = \frac{ ext{MS}_{ ext{between}}}{ ext{MS}_{ ext{within}}}$$ If the omnibus $F$-test is significant ($p < 0.05$), we proceed to Tukey's Honest Significant Difference (Tukey HSD) to determine which specific pairs differ while controlling FWER at $0.05$.
Hands-On R Script: Running ANOVA and Tukey HSD in R
library(tidyverse)
# Simulated Agricultural Runoff Experiment: Nitrate (mg/L) across 4 Riparian Buffer Treatments
set.seed(77)
buffer_data <- tibble(
treatment = factor(rep(c("Control_NoBuffer", "Grass_Filter", "Forested_Buffer", "Wetland_Pond"), each = 15)),
nitrate = c(rnorm(15, mean = 14.8, sd = 2.1),
rnorm(15, mean = 9.4, sd = 1.9),
rnorm(15, mean = 5.2, sd = 1.6),
rnorm(15, mean = 3.8, sd = 1.4))
)
# 1. Fit One-Way ANOVA Model
anova_model <- aov(nitrate ~ treatment, data = buffer_data)
summary(anova_model)
# Note: F(3, 56) = 114.2, p < 2e-16 -> Extremely significant treatment effect!
# 2. Check ANOVA Residual Normality and Homoscedasticity
shapiro.test(residuals(anova_model)) # p > 0.05 -> Residuals are normal!
# 3. Post-Hoc Pairwise Comparisons with Tukey HSD
tukey_res <- TukeyHSD(anova_model)
print(tukey_res)
# 4. Non-Parametric Fallback: Kruskal-Wallis Test
# If normality fails:
kruskal_res <- kruskal.test(nitrate ~ treatment, data = buffer_data)
print(kruskal_res)Common Pitfalls & Troubleshooting Checklist
- Java Heap Exhaustion (GATK / Picard): If processes terminate with
java.lang.OutOfMemoryError, set your environment variable explicitly:export _JAVA_OPTIONS="-Xmx16g -Djava.io.tmpdir=/tmp". - Index Timestamp Mismatch: If tools complain that a BAM index is missing or out of date, regenerate with
samtools index -@ 4 sample.bam. - Conda Channel Conflicts: Avoid ABI incompatibilities between bioconda and conda-forge by setting strict priority:
conda config --set channel_priority strict. - Chromosome Naming Contamination: Always verify whether reference FASTA headers use
chr1or1; divergent naming causes silent alignment drops.