ggplot(cu_wq, aes(year_frac, result)) +
geom_point(alpha = 0.5) +
geom_smooth(se = FALSE) +
labs(title = "Total Copper — raw scale", y = "Copper (µg/L)") +
nlt_theme
Generalized Additive Models
Total Copper — step-like, non-monotonic change with no fixed pattern
1 What a GAM is
Total Copper doesn’t settle into flat-then-sloped segments (piecewise regression), a single level jump (changepoint detection), or a repeating annual shape with a correlated-noise problem (GAMM, next chapter, on Dissolved Oxygen). It moves up, up again into a brief spike, back down, further down below where it started, then partway back up — a shape with no periodic structure and no small set of knots an analyst would confidently place by eye. That’s the case a GAM is actually for.
1.1 Still linear, in the sense that matters
A GAM’s linear predictor is \(\eta = \beta_0 + f(x)\), where \(f\) is a smooth function built from a spline basis rather than a single coefficient times \(x\) — a weighted sum of many small, fixed basis columns, \(f(x) = \sum_{m=1}^{k} \beta_m b_m(x)\), the same “linear in the parameters” idea covered in The Regression Framework and used there for the piecewise chapter’s hinge term. What’s different is how much basis-building mgcv::gam() does automatically: instead of one or two hand-picked hinge terms, it builds a rich basis (up to k candidate knots) and a fitted smoothing penalty decides how much of that flexibility survives, shrinking toward a straight line wherever the data doesn’t support the extra wiggle.
1.2 How the penalty works, concretely
Fitting isn’t just least squares on the basis columns — the coefficients \(\beta_1, \dots, \beta_k\) are chosen to minimize
\[\sum_{i=1}^n \left(y_i - f(x_i)\right)^2 + \lambda \int f''(x)^2\,dx,\]
a squared-error term traded off against a roughness penalty: \(\int f''(x)^2\,dx\) is large when \(f\) bends sharply and zero when \(f\) is a straight line, since a line has zero second derivative everywhere. \(\lambda\) sets the exchange rate between the two — \(\lambda \to 0\) lets the fit use as much of the basis’s flexibility as it wants, \(\lambda \to \infty\) forces \(f\) back to a straight line regardless of what the basis could represent. method = "REML" below estimates \(\lambda\) from the data itself, rather than fixing it or choosing it by eye.
1.3 What edf means
The effective degrees of freedom (edf) of a fitted smooth is \(\mathrm{edf} = \mathrm{tr}(S_\lambda)\), the trace of the smoother matrix \(S_\lambda\) that maps the observed \(y\)’s onto the fitted \(\hat f\)’s at the fitted \(\lambda\). It generalizes “number of parameters” from an ordinary linear term — a straight line has exactly one slope, edf = 1 — to a penalized smooth, where \(\lambda\) can land anywhere between edf ≈ 1 (penalized all the way back to a straight line) and edf ≈ k (using nearly all of the basis’s allowed flexibility). It’s the one number, reported alongside every smooth term on this site, that says how much curvature the fit actually kept.
1.4 Semi-parametric, not fully either
The basis dimension k caps how wiggly the fit is allowed to get; the fitted penalty decides how much of that allowance actually gets used, based on the data. That’s closer to “as flexible as the data supports” than either a fixed parametric form (piecewise regression) or no form at all (k-NN, LOESS).
1.5 Inferential, not just a curve
The fitted smooth carries an edf and a significance test — is this smooth distinguishable from flat at all — the same kind of reportable number a slope or a breakpoint gives elsewhere in this vignette. That’s what makes “is this concentration currently rising, falling, or holding” an answerable, testable question below, rather than a read off a noisy plot.
2 The shape in the record
2.1 Raw series
An early low, stable period; a rise to a higher plateau; a brief spike well above that; a drop back down; a further drop below the original level; and a partial recovery — none of it periodic, none of it a straight line or two.
2.2 Checking the scale: log vs. raw
fit_check_log <- gam(log(result) ~ s(year_frac, bs = "tp", k = 20), data = cu_wq, method = "REML")
fit_check_raw <- gam(result ~ s(year_frac, bs = "tp", k = 20), data = cu_wq, method = "REML")
bptest(lm(resid(fit_check_log) ~ fitted(fit_check_log)))
studentized Breusch-Pagan test
data: lm(resid(fit_check_log) ~ fitted(fit_check_log))
BP = 1.426, df = 1, p-value = 0.2324bptest(lm(resid(fit_check_raw) ~ fitted(fit_check_raw)))
studentized Breusch-Pagan test
data: lm(resid(fit_check_raw) ~ fitted(fit_check_raw))
BP = 50.418, df = 1, p-value = 1.242e-12Total Copper is a concentration, so the usual dilution/loading argument for logging applies, and the check confirms it: the raw-scale fit fails Breusch-Pagan decisively (\(p = 1.2\times10^{-12}\)), the log-scale fit doesn’t (\(p = 0.232\)). log(result) is used from here on.
3 Fitting the shape
fit_cu <- gam(log(result) ~ s(year_frac, bs = "tp", k = 20), data = cu_wq, method = "REML")
summary(fit_cu)
Family: gaussian
Link function: identity
Formula:
log(result) ~ s(year_frac, bs = "tp", k = 20)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.01387 0.01998 100.8 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(year_frac) 15.91 17.87 37.56 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.824 Deviance explained = 84.4%
-REML = 32.922 Scale est. = 0.057484 n = 144k = 20 gives the basis enough candidate knots to represent several sharp transitions if the data supports them — a small k would cap the fit’s flexibility below what this record’s shape actually needs, forcing the penalty to smooth over real transitions it doesn’t have the basis capacity to represent. The fitted edf of 15.9 shows the penalty used most of that allowance rather than shrinking toward a straight line, and the smooth term is significant by a wide margin (\(p < 2\times10^{-16}\)).
4 The fitted curve, with a confidence band
grid <- data.frame(year_frac = seq(min(cu_wq$year_frac), max(cu_wq$year_frac), length.out = 300))
pred <- predict(fit_cu, newdata = grid, se.fit = TRUE)
grid <- grid |>
mutate(
fit = pred$fit,
lower = exp(fit - 1.96 * pred$se.fit),
upper = exp(fit + 1.96 * pred$se.fit),
fit = exp(fit)
)
ggplot() +
geom_point(data = cu_wq, aes(year_frac, result), alpha = 0.4, color = "grey40") +
geom_ribbon(data = grid, aes(year_frac, ymin = lower, ymax = upper), fill = "steelblue", alpha = 0.25) +
geom_line(data = grid, aes(year_frac, fit), color = "steelblue", linewidth = 1) +
labs(title = "Total Copper — fitted GAM with 95% confidence band", y = "Copper (µg/L)") +
nlt_theme
predict(fit_cu, se.fit = TRUE) returns a standard error alongside the fitted value at every point, the same as predict.lm() would — a GAM’s confidence band isn’t a special case, it’s the ordinary consequence of the smooth term being a fitted linear combination of coefficients with a covariance matrix, propagated through exp() for the raw-scale band shown here. The band widens over the brief spike around year 4.5–5, where only a few months of data inform that stretch of the curve, and the fitted peak (about 18.7 µg/L) sits below the highest individual raw readings there — a short, few-month excursion is exactly what a smoothing penalty will partly discount, since a spike consistent with only a handful of points is harder to distinguish from noise than a shift that persists across many.
5 Where the record is rising, falling, or flat
A single edf and p-value confirm the smooth term overall isn’t flat, but say nothing about where the record is moving in a given direction, which is the question a regulator actually asks. That question is answered by the derivative of the fitted curve — its local slope — and a confidence interval on that derivative, estimated by the same finite-difference approach used in the GAMM chapter: nudge the fitted curve forward and backward by a small step, take \((\hat y(t+h) - \hat y(t-h))/2h\) as the local slope, and get its standard error from predict(..., type = "lpmatrix") propagated through the model’s coefficient covariance.
5.1 Computing the local slope
h <- 1e-3
Xp_plus <- predict(fit_cu, newdata = data.frame(year_frac = grid$year_frac + h), type = "lpmatrix")
Xp_minus <- predict(fit_cu, newdata = data.frame(year_frac = grid$year_frac - h), type = "lpmatrix")
D <- (Xp_plus - Xp_minus) / (2 * h)
beta <- coef(fit_cu)
Vp <- vcov(fit_cu)
deriv <- as.numeric(D %*% beta)
se <- sqrt(rowSums((D %*% Vp) * D))
grid <- grid |>
mutate(
deriv_lower = deriv - 1.96 * se,
deriv_upper = deriv + 1.96 * se,
trend = case_when(
deriv_lower > 0 ~ "rising",
deriv_upper < 0 ~ "falling",
TRUE ~ "flat"
)
)
table(grid$trend)
falling flat rising
45 179 76 5.2 Visualizing significant stretches
grid_seg <- grid |>
mutate(year_frac_end = lead(year_frac), fit_end = lead(fit)) |>
filter(!is.na(year_frac_end))
trend_colors <- c(rising = "#2E7D32", falling = "#C62828", flat = "grey55")
ggplot() +
geom_point(data = cu_wq, aes(year_frac, result), alpha = 0.35, color = "grey40") +
geom_segment(
data = grid_seg,
aes(x = year_frac, y = fit, xend = year_frac_end, yend = fit_end, color = trend),
linewidth = 1.1
) +
scale_color_manual(values = trend_colors, name = "local trend") +
labs(title = "Total Copper — fitted trajectory by significance of local slope", y = "Copper (µg/L)") +
nlt_theme +
theme(legend.position = "right")
The coloring tracks the shape described from the raw plot, but now with each stretch labeled by whether the change is statistically distinguishable from flat: a significant rise onto the first plateau, a significant further rise into the spike, a significant fall back down, a significant fall below the original baseline, a significant partial recovery, and flat stretches — not significantly different from zero slope — everywhere the record is actually holding at a level rather than moving. That’s the concrete answer to “is it currently getting better or worse”: at any point on this record, the color at that point is a tested answer, not a visual impression.
6 Model diagnostics
6.1 Homoscedasticity and normality
bptest(lm(resid(fit_cu) ~ fitted(fit_cu)))
studentized Breusch-Pagan test
data: lm(resid(fit_cu) ~ fitted(fit_cu))
BP = 1.426, df = 1, p-value = 0.2324shapiro.test(resid(fit_cu))
Shapiro-Wilk normality test
data: resid(fit_cu)
W = 0.98827, p-value = 0.2661Both tests fail to reject their null: the Breusch-Pagan test (\(p = 0.232\)) finds no evidence of heteroscedasticity, and the Shapiro-Wilk test (\(p = 0.266\)) finds no evidence of non-normal residuals, on the log scale settled on above.
The GAM’s own F-test on the smooth term rests on the same independence assumption as any OLS-family test, so it’s worth checking directly rather than assuming it away because this fit isn’t the seasonal-cycle case GAMM was built to catch:
6.2 Independence
acf(resid(fit_cu), main = "Total Copper GAM — residual autocorrelation")
Box.test(resid(fit_cu), lag = 1, type = "Ljung-Box")
Box-Ljung test
data: resid(fit_cu)
X-squared = 13.558, df = 1, p-value = 0.0002313Unlike the piecewise and changepoint chapters’ clean results, this one isn’t clean: lag-1 autocorrelation is a little over 0.3, and the Ljung-Box test rejects independence decisively (\(p = 2.3\times10^{-4}\)).
6.3 Residuals over time
cu_wq |>
mutate(resid = resid(fit_cu)) |>
ggplot(aes(year_frac, resid)) +
geom_point(alpha = 0.5) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "Total Copper GAM — residuals over time", y = "residual (log scale)") +
nlt_theme
Plotted against time, the largest residuals aren’t spread evenly across the record — they cluster around the same stretches already flagged as sharp transitions in “The shape in the record” above, while the flatter stretches on either side of them show tight, patternless scatter. That’s a different mechanism than GAMM’s case: a penalized spline is constrained to be smooth, so wherever the true mean genuinely bends fast — a step, a spike — the fitted curve necessarily lags or overshoots for a few points to either side, producing runs of same-signed residuals right at the transition without any correlated noise process behind them. It’s a real limitation worth naming rather than a false alarm, but it’s concentrated at the transitions rather than spread through the whole record, which is why the overall significance conclusion (\(p < 2\times10^{-16}\), an effect far too large for a handful of locally correlated points to manufacture) still stands. It does mean the confidence band and the derivative-significance call right at each transition edge should be read with a bit more caution than the flat, well-fit stretches between them.
7 Why not just engineer the right OLS features
The GAM’s spline basis is still feature engineering handed to a linear fit, the same category of move as Piecewise Regression’s hand-built hinge term. Unlike Dissolved Oxygen’s clean annual cycle, though, there’s no sin/cos shortcut here — nothing about a step-and-spike pattern is periodic. The honest OLS alternative is a plain polynomial in year_frac, tried at increasing degree:
poly_fits <- lapply(c(4, 8, 12, 16), function(deg) lm(log(result) ~ poly(year_frac, deg), data = cu_wq))
data.frame(
model = c(sprintf("OLS, degree-%d polynomial", c(4, 8, 12, 16)), "GAM"),
AIC = c(sapply(poly_fits, AIC), AIC(fit_cu)),
R2_or_dev_explained = c(sapply(poly_fits, function(f) summary(f)$r.squared), summary(fit_cu)$dev.expl)
) model AIC R2_or_dev_explained
1 OLS, degree-4 polynomial 173.57620 0.4453331
2 OLS, degree-8 polynomial 89.32471 0.7077141
3 OLS, degree-12 polynomial 65.88346 0.7650459
4 OLS, degree-16 polynomial 56.34146 0.7919932
5 GAM 16.53372 0.8435067Even at degree 16 — well past where a polynomial starts ringing at the record’s edges (Runge’s phenomenon) rather than fitting anything real — the polynomial’s \(R^2\) (0.79) and AIC (56.3) both trail the GAM’s (0.84 deviance explained, AIC 16.5) by a wide margin, and the gap isn’t closing quickly as degree increases. The reason is structural, not just a matter of trying harder: a polynomial term is global — every coefficient shapes the curve’s behavior everywhere, so fitting a sharp local step or spike well tends to distort the fit far away from it, or costs an escalating number of parameters trying not to. A spline basis is local — each basis function’s influence is confined to a neighborhood of its knot, so representing one sharp transition doesn’t fight with representing another one elsewhere in the record. That locality, not just a higher effective degree, is what a hand-built polynomial feature can’t reproduce no matter how high its degree goes.
8 Interpretation
Total Copper moved through five distinct, individually significant phases over the twelve-year record: a rise onto an elevated plateau, a further rise into a brief spike, a fall back down, a fall below the original baseline, and a partial recovery — each stretch tested, not eyeballed, via the derivative-significance plot above, and none of it describable by a straight line, a single knot, or a periodic function. The GAM’s smoothing penalty found this shape without being told any of it in advance, using only a smoothness assumption, and a degree-16 polynomial handed every opportunity to match it still fell well short — the concrete case, in this vignette, for reaching past a fixed parametric form when a series’ shape isn’t known ahead of time. For a regulator asking “is this currently a problem,” the record closes on the flat, partially-recovered plateau rather than the elevated spike, which is a materially different answer than either the raw high-water mark or a single straight-line trend fit to the whole record would give.
9 Extending this analysis
s(year_frac) alone is the single-predictor case; a GAM’s real strength shows up once more than time is believed to drive the analyte. If a plausible confound exists — flow, a co-occurring pollutant, an upstream operational variable — the natural extension is a multi-predictor GAM, s(year_frac) + s(flow) or similar, which lets each smooth absorb its own share of the variation instead of forcing everything through a single time trend the way this chapter does. The same machinery generalizes to multiple monitoring sites too, via a by-factor smooth (s(year_frac, by = site)) that lets the fitted shape differ across sites while still sharing one model, rather than fitting each site’s record separately the way this single-composite-site case study does.