library(dplyr)
library(ggplot2)
library(knitr)
library(DT)
source("R/utils.R")
source("R/trend_model.R")
source("R/mann_kendall.R")
source("R/plot_trend.R")

Background

Environmental monitoring programs accumulate years of water quality data across multiple sampling areas. Assessing whether a given analyte is trending up, down, or holding steady over that record — and whether the trend differs by area — informs everything from compliance reporting to the design of the next monitoring plan.

The same underlying question — is this analyte trending, and does the trend hold across every area — comes up repeatedly across different projects and consulting engagements, usually re-implemented from scratch each time. This vignette is also an example of streamlining that common logic into one reusable pipeline rather than a new script per analyte or project. That streamlining isn’t a substitute for looking closely at each analyte, though: outlier removal, assumption checks, and the interaction test below are still run and inspected individually for TSS, Copper, Ammonia, and Iron, not applied blindly. The benefit is in not re-deriving the same statistical logic every time, while still giving each analyte its own answer.

This vignette walks through a reproducible year-over-year trend assessment applied to four case-study analytes, each chosen to exercise a different part of the pipeline:

  1. TSS — a genuine monotonic increase common to every monitoring area (a main effect of year).
  2. Total Copper — flat across the record in every area (nothing to report — the contrast case).
  3. Total Ammonia — a trend confined to one area (Effluent) that other areas don’t share, producing a significant area × year interaction.
  4. Total Iron — flat overall, but with a handful of high-leverage outliers that need to be identified and removed before the model is trustworthy.

For each analyte the pipeline:

  1. Fits a two-way ANOVA (result ~ year * area) and iteratively removes studentized-residual outliers (|residual| > 3).
  2. Checks model assumptions: homoscedasticity (Breusch-Pagan) and normality of residuals (Shapiro-Wilk).
  3. Tests the area × year interaction term. If it’s significant, area-specific trends genuinely differ and the full interaction model is kept, followed by emmeans post-hoc pairwise comparisons between areas. If not, the model is reduced to the additive year + area form.
  4. Runs a Mann-Kendall trend test per area — a nonparametric, assumption-light check on trend direction that complements the ANOVA slope.

A significant year effect only says the slope is distinguishable from zero — it doesn’t say how big the trend is. Because year is fit as a numeric predictor rather than a factor, the model already contains that slope; the pipeline reads it back out and reports it as a %/year rate of change alongside every significance test, and the Mann-Kendall step reports its own nonparametric equivalent (Sen’s slope) alongside tau, for the same reason.

The data below are synthetic; program names and site identifiers have been genericized to three monitoring areas (Near Field, Mid Field, Effluent) so the vignette runs standalone.

This report was developed with the assistance of Claude, refactoring and streamlining analysis code and scripts originally written by hand for individual projects into the shared pipeline described below.

Pipeline design

The four analyte sections below aren’t four separate scripts. Each one is the same handful of function calls, defined once in R/ and reused for every analyte:

  • fit_trend_model() (R/trend_model.R) does the statistical work: removes outliers, fits the two-way ANOVA, checks the model assumptions, tests the area × year interaction, and — depending on that test — either fits the additive model or keeps the interaction model and runs the emmeans post-hoc comparison. It returns a single object holding everything the report needs downstream: the cleaned data, the fitted model, the interaction test, the assumption checks, and the year-effect slope.
  • mann_kendall_summary() (R/mann_kendall.R) runs the nonparametric Mann-Kendall test and Sen’s slope per area.
  • plot_trend(), plot_posthoc(), and plot_overview() (R/plot_trend.R), sharing one theme and color palette from R/utils.R, turn a fit into the figures below.
  • summarize_trend() collapses a fit to one row for the cross-analyte summary table in the final section.

Each section’s code is just fit_trend_model(temporal, "<analyte>") followed by calls to these helpers. The branching — additive vs. interaction, post-hoc or not, how many outliers to drop — is a property of the data, decided once inside fit_trend_model(), not something written per analyte. Of the four analytes below, only Ammonia’s data triggers the interaction branch; Copper, TSS, and Iron all reduce to the additive model. That asymmetry is exactly why the pipeline is centralized rather than copy-pasted: a fixed per-analyte script has to either special-case the interaction/post-hoc logic every time — most of which then goes unused — or risk silently skipping it on the one analyte where a real interaction shows up. One function tests for it on every analyte and only pays for the post-hoc step when the data actually calls for it.

Additive vs. interaction models

Additive model — year + area. Every area gets its own line, but all lines share the same slope — they’re parallel, just shifted up or down by a constant, area-specific offset that’s the same in every year:

\[\log(\text{result}) = \beta_0 + \beta_1 \cdot \text{year} + \beta_2 \cdot \text{area}\]

Interaction model — year * area. Adds a per-area adjustment to the slope itself, so the lines no longer have to be parallel — one area’s trend can rise while another stays flat:

\[\log(\text{result}) = \beta_0 + \beta_1 \cdot \text{year} + \beta_2 \cdot \text{area} + \beta_3 \cdot (\text{year} \times \text{area})\]

Glossary of tests

Glossary of statistical tests used in the pipeline
Term What it does Plain-language meaning
Breusch-Pagan test Tests whether residual variance is constant across fitted values (homoscedasticity). A low p-value flags heteroscedasticity. “Is the model’s scatter roughly the same size everywhere, or does it fan out?”
Shapiro-Wilk test Tests whether model residuals are normally distributed. A low p-value flags non-normal residuals. “Do the leftover errors look bell-shaped, like they should?”
emmeans Computes model-based estimated marginal means and pairwise contrasts between groups (here, areas at the most recent year), with multiple-comparison-adjusted p-values. “After accounting for year and area, which areas actually differ from each other, and by how much?”
Mann-Kendall test A nonparametric test for a monotonic trend over time, run per area. Doesn’t assume normality or a linear relationship — a lightweight cross-check on the ANOVA slope. Paired with Sen’s slope, its nonparametric estimate of trend magnitude. “Ignoring the model entirely, is this area’s data going up, down, or flat over time, and by roughly how much per year?”

1. Data Overview

temporal <- read.csv("data/synthetic/temporal_wq.csv")

Three of the four analytes look broadly similar at a glance — scattered points across ten years, colored by area. The statistical tests below are what actually separate “a real trend” from “no trend” from “a trend that only shows up in one area.”

2. TSS — A Real Trend

TSS rises at a similar rate in every monitoring area, so the area × year interaction should come back non-significant and the pipeline should fall through to the additive model.

tss_fit <- fit_trend_model(temporal, "TSS")
TSS model summary
Analyte Outliers removed Model Interaction p-value Breusch-Pagan p-value Shapiro-Wilk p-value
TSS 1 additive 0.174 0.0613 0.998

TSS Mann-Kendall trend test by area
Area tau p-value Sen’s slope (%/year) 95% CI (%/year) n
Effluent 0.297 7.12e-03 2.99 0.997 to 5.23 40
Near Field 0.146 1.88e-01 1.75 -0.824 to 4.28 40
Mid Field 0.501 7.60e-06 4.79 3.08 to 6.54 39

The interaction term is not significant (p = 0.174), confirming a shared upward trend rather than area-specific behavior. That test only says the shared slope is distinguishable from zero, though — the size of it comes from the year coefficient in the additive model, back-transformed to a %/year rate:

TSS ANOVA slope (shared across areas)
Area Slope (log/year) Slope (%/year) 95% CI (%/year) p-value
All areas 0.0325 3.31 2.13 to 4.5 1e-07

TSS is rising at roughly 3.31% per year. The Mann-Kendall tau values are positive in every area, consistent with the ANOVA result, and its Sen’s slope estimates (also above, in the per-area Mann-Kendall table) land in the same range without assuming a linear model or normal residuals — a nonparametric cross-check on both the direction and the magnitude of the trend.

3. Total Copper — No Trend

Copper concentrations are flat over the record in every area. This is the contrast case: no significant year effect, no interaction, nothing to report.

copper_fit <- fit_trend_model(temporal, "Total Copper")
Total Copper model summary
Analyte Outliers removed Model Interaction p-value Breusch-Pagan p-value Shapiro-Wilk p-value
Total Copper 0 additive 0.752 0.907 0.147

Total Copper Mann-Kendall trend test by area
Area tau p-value Sen’s slope (%/year) 95% CI (%/year) n
Effluent -0.0205 0.861 -0.298 -2.2 to 1.49 40
Near Field 0.0179 0.880 0.112 -1.46 to 1.96 40
Mid Field 0.1150 0.300 0.861 -1.1 to 2.45 40
Total Copper ANOVA slope (shared across areas)
Area Slope (log/year) Slope (%/year) 95% CI (%/year) p-value
All areas 0.00277 0.277 -0.732 to 1.3 0.589

Mann-Kendall p-values are all well above 0.05 in every area — no evidence of a monotonic trend, consistent with the flat ANOVA fit. The year coefficient’s confidence interval also spans zero, so there’s no slope to report here beyond noise — the non-significant p-value and the near-zero estimate agree.

4. Total Ammonia — A Trend Confined to One Area

Ammonia is where the pipeline’s interaction branch actually matters: Effluent concentrations rise steadily while Near Field and Mid Field stay flat. A model that only looked at the overall year effect would miss this; the area × year interaction term is what catches it.

ammonia_fit <- fit_trend_model(temporal, "Total Ammonia")
Total Ammonia model summary
Analyte Outliers removed Model Interaction p-value Breusch-Pagan p-value Shapiro-Wilk p-value
Total Ammonia 0 interaction 0 0.439 0.985

Before deciding which model to report, the pipeline tests the null hypothesis that the area × year interaction term is zero — i.e., that all three areas share the same trend — against the alternative that at least one area’s trend differs, using a Type III F-test (car::Anova) at α = 0.05. If the interaction p-value is at or above 0.05, that term is dropped and the simpler additive model (year + area) is reported instead, since there’s no statistical evidence the trends differ; because the additive model’s area offset is constant across years, a single ANOVA main-effect test for area already tells you whether areas differ overall, so there’s no need for a post-hoc, year-specific comparison. If it’s below 0.05, as is the case here, the full interaction model is kept and emmeans post-hoc pairwise comparisons are run between areas at the most recent year, to identify specifically which areas differ once the shared trend is allowed to vary by area.

The Type III ANOVA table below has one row per term: the year and area main effects, and the year:area interaction. When the interaction is significant, it’s the row that matters most — a significant year:area means the year and area main-effect rows are no longer straightforwardly interpretable on their own (the “effect of year” isn’t a single number once it’s allowed to differ by area), so the interaction result is what actually drives what gets reported below:

Total Ammonia Type III ANOVA
Sum Sq Df F value Pr(>F)
(Intercept) 2.1272 1 71.5583 0
year 2.2171 1 74.5823 0
area 1.4922 2 25.0985 0
year:area 1.5106 2 25.4089 0
Residuals 3.3888 114 NA NA

The interaction is significant (p = 7.49^{-10}), so the pipeline keeps the full interaction model and runs emmeans post-hoc comparisons between areas at the most recent year, rather than reducing to the additive model:

Total Ammonia emmeans post-hoc pairwise comparisons, most recent year
contrast estimate SE df t.ratio p.value
Effluent - Mid Field 1.6163 0.0717 114 22.5561 0
Effluent - Near Field 1.2886 0.0717 114 17.9827 0
Mid Field - Near Field -0.3277 0.0717 114 -4.5735 0

Because the interaction is significant, there’s no single shared slope to report — each area gets its own, estimated from the fitted model via emmeans::emtrends():

Total Ammonia ANOVA slope by area
Area Slope (log/year) Slope (%/year) 95% CI (%/year) p-value
Effluent 0.08200 8.540 6.52 to 10.6 0.000
Near Field -0.00676 -0.673 -2.52 to 1.21 0.478
Mid Field 0.00658 0.660 -1.21 to 2.57 0.490
Total Ammonia Mann-Kendall trend test by area
Area tau p-value Sen’s slope (%/year) 95% CI (%/year) n
Effluent 0.6380 0.000 8.440 6.52 to 10.4 40
Near Field -0.0513 0.650 -0.565 -2.89 to 1.48 40
Mid Field 0.0667 0.552 0.703 -1.17 to 2.68 40

The Mann-Kendall results confirm the same asymmetry at the individual-area level: a strong, significant upward trend in Effluent, and no trend in Near Field or Mid Field. Its Sen’s slope for Effluent lands close to the ANOVA estimate above, giving the same magnitude by an assumption-light route.

5. Total Iron — Outlier Removal in Action

Iron is flat like copper, but the record includes a handful of erroneously high results — the kind of data-entry or instrument spike that shows up in any long-running monitoring dataset. fit_trend_model() flags these via studentized residuals and refits until none remain.

iron_fit <- fit_trend_model(temporal, "Total Iron")
Total Iron model summary
Analyte Outliers removed Model Interaction p-value Breusch-Pagan p-value Shapiro-Wilk p-value
Total Iron 5 additive 0.64 0.608 0.226

5 observation(s) were removed as outliers before the reported model was fit.

Total Iron Mann-Kendall trend test by area
Area tau p-value Sen’s slope (%/year) 95% CI (%/year) n
Effluent -0.0901 0.440 -0.961 -2.59 to 1.14 37
Near Field 0.0850 0.453 0.703 -1.29 to 2.7 39
Mid Field 0.0472 0.681 0.439 -1.93 to 2.87 39
Total Iron ANOVA slope (shared across areas)
Area Slope (log/year) Slope (%/year) 95% CI (%/year) p-value
All areas -0.00156 -0.155 -1.27 to 0.969 0.784

As with Copper, the slope estimate is small and its confidence interval spans zero — flat, once the outliers are out of the way.

6. Summary Across Analytes

Model summary across all four analytes
Analyte Outliers removed Model Interaction p-value Breusch-Pagan p-value Shapiro-Wilk p-value
TSS 1 additive 0.174 0.0613 0.998
Total Copper 0 additive 0.752 0.9070 0.147
Total Ammonia 0 interaction 0.000 0.4390 0.985
Total Iron 5 additive 0.640 0.6080 0.226

TSS shows a real, shared upward trend; Copper shows none; Ammonia’s trend is confined to Effluent; Iron is flat once its outliers are removed. Four different outcomes, produced by running the same pipeline described at the outset against each analyte in turn — not four different analyses.