10. Capstone Project: End-to-End Environmental Case Study & Reproducible Reporting

In this capstone module, you synthesize all foundational statistical concepts into an end-to-end scientific research project. You will analyze a real-world multi-site watershed investigation assessing the impact of urbanization and agricultural runoff on river ecosystem health.

Complete Production R Script: The Watershed Case Study

# Complete Capstone Watershed Analysis Script
library(tidyverse)

# 1. Generate Synthetic Multi-Site Watershed Dataset (n = 80 samples)
set.seed(2026)
watershed <- tibble(
  sample_id  = paste0("SAMP_", sprintf("%03d", 1:80)),
  zone       = factor(rep(c("Pristine_Headwater", "Agricultural_Midstream", "Urban_Downstream", "Restored_Wetland"), each = 20)),
  turbidity  = c(rnorm(20, 2.1, 0.4), rnorm(20, 8.4, 1.2), rnorm(20, 19.5, 3.1), rnorm(20, 4.2, 0.8)), # NTU
  nitrate    = c(rnorm(20, 0.35, 0.08), rnorm(20, 6.2, 1.1), rnorm(20, 8.9, 1.4), rnorm(20, 1.8, 0.4)), # mg/L
  bio_index  = c(rnorm(20, 88, 6), rnorm(20, 52, 8), rnorm(20, 24, 7), rnorm(20, 71, 5)) # EPT Biotic Index
)

# 2. Descriptive Statistics Summary Table
stats_summary <- watershed %>%
  group_by(zone) %>%
  summarise(
    n           = n(),
    mean_bio    = mean(bio_index),
    sd_bio      = sd(bio_index),
    sem_bio     = sd(bio_index) / sqrt(n()),
    mean_nitrate = mean(nitrate),
    sd_nitrate  = sd(nitrate)
  )
print(stats_summary)

# 3. Hypothesis Testing: Does Macroinvertebrate Health Differ Significantly by Zone?
anova_bio <- aov(bio_index ~ zone, data = watershed)
summary(anova_bio)

# Check Normality of Model Residuals
shapiro_res <- shapiro.test(residuals(anova_bio))
cat("Shapiro-Wilk p-value:", shapiro_res$p.value, "
")

# Post-Hoc Comparisons
tukey_bio <- TukeyHSD(anova_bio)
print(tukey_bio)

# 4. Regression Analysis: Predicting Biotic Health from Nutrient Pollution
model_nitrate_bio <- lm(bio_index ~ nitrate, data = watershed)
summary(model_nitrate_bio)

# 5. Publication Visualization with ggplot2
fig_a <- ggplot(watershed, aes(x = zone, y = bio_index, fill = zone)) +
  geom_boxplot(alpha = 0.7, width = 0.5, outlier.shape = NA) +
  geom_jitter(width = 0.15, alpha = 0.8, size = 2, shape = 21, color = "black") +
  theme_bw(base_size = 11) +
  theme(legend.position = "none", axis.text.x = element_text(angle = 20, hjust = 1)) +
  labs(title = "A: Biotic Integrity Across River Zones", y = "Macroinvertebrate Biotic Index", x = "")

fig_b <- ggplot(watershed, aes(x = nitrate, y = bio_index)) +
  geom_point(aes(color = zone), size = 2.5, alpha = 0.85) +
  geom_smooth(method = "lm", color = "black", se = TRUE) +
  theme_bw(base_size = 11) +
  labs(title = "B: Biotic Health vs. Nitrate Pollution", x = "Nitrate Concentration (mg/L)", y = "Biotic Index")

# Print figures
print(fig_a)
print(fig_b)

How to Report Statistical Results in Scientific Journals