ggplot(do_wq, aes(year_frac, result)) +
geom_point(alpha = 0.5) +
geom_smooth(se = FALSE) +
labs(title = "Dissolved Oxygen — raw scale", y = "DO (mg/L)") +
nlt_theme
Generalized Additive Mixed Models
Dissolved Oxygen — checking whether the GAM’s significance test can be trusted
1 What changes between a GAM and a GAMM
1.1 Why the GAM’s F-test assumes independence
A GAM’s significance test on a smooth term — the F-test behind an edf and a p-value — is built on the same assumption ordinary least squares makes: the residuals are independent draws. Nothing about fitting s(year_frac) and s(month) checks that assumption, and monthly water-quality data routinely violates it — this month’s reading tends to sit closer to last month’s than an independent draw would, for reasons the model’s mean structure hasn’t captured (slow-moving hydrology, persistence in whatever’s driving oxygen demand). If that correlation is there and unmodeled, the standard errors on a smooth term come out too small, edf can come out inflated, and “the trend is significant” can mean “the trend, plus some of the correlated noise the model mistook for it.”
1.2 What a GAMM adds: an AR(1) correlation structure
A generalized additive mixed model keeps the exact same additive-smooth structure — the same s(year_frac) + s(month), the same basis functions, the same idea of a fitted penalty deciding how much wiggle survives — and adds one thing: an explicit model for how the residuals are correlated, rather than assuming they aren’t. mgcv::gamm() fits this by handing the smooth-term machinery to nlme::lme() underneath, which accepts a correlation argument — here, corAR1(), an AR(1) structure where the correlation between any two residuals decays geometrically with how far apart in time they are:
\[\mathrm{Corr}(\varepsilon_i, \varepsilon_{i-h}) = \phi^{|h|}, \qquad h = 0, 1, 2, \dots\]
A single parameter \(\phi\) (with \(|\phi| < 1\) for the process to be stationary) sets that whole decay: \(h = 1\) apart (adjacent months) gives correlation \(\phi\) itself, \(h = 2\) apart gives \(\phi^2\), and so on, fading toward zero for observations far apart rather than dropping to zero immediately after one lag. Nothing about the mean structure changes; what changes is that the smooth terms’ significance tests now account for observations not being independent, which is the actual point of this chapter.
2 The shape in the record
The same shape as before: an annual saw-tooth with a slight downward drift, hard to judge by eye against the size of the seasonal swing.
2.1 Does log(result) still make sense here
Dissolved Oxygen isn’t a concentration governed by dilution or loading the way Total Nitrogen or Specific Conductance are — it’s set by gas solubility, reaeration, and biological demand — so the usual multiplicative argument for logging concentration data doesn’t obviously carry over. Checked empirically instead: fitting the same additive-smooth structure on both scales and comparing residual homoscedasticity, the raw-scale fit fails the Breusch-Pagan test (\(p = 0.0024\)) and the log-scale fit doesn’t (\(p = 0.133\)). log(result) is used from here on, on the strength of that check rather than the weaker a priori argument.
3 Fitting trend and season, assuming independence first
fit_gam <- gam(log(result) ~ s(year_frac, bs = "tp", k = 10) + s(month, bs = "cc", k = 12),
data = do_wq, knots = list(month = c(0.5, 12.5)), method = "REML")
summary(fit_gam)
Family: gaussian
Link function: identity
Formula:
log(result) ~ s(year_frac, bs = "tp", k = 10) + s(month, bs = "cc",
k = 12)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.136958 0.004466 478.5 <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) 3.846 4.756 16.66 <2e-16 ***
s(month) 6.626 10.000 155.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.923 Deviance explained = 92.9%
-REML = -195.28 Scale est. = 0.0028726 n = 144Read at face value, s(year_frac) looks like real nonlinear structure — edf 3.85, not a straight line — on top of a strongly significant seasonal term. Before trusting that, it’s worth checking the assumption the F-test above depends on.
4 Checking the independence assumption
acf(resid(fit_gam), main = "Residual autocorrelation, naive GAM")
Lag-1 autocorrelation in the residuals is a little over 0.5 — the residual left over after month and year_frac’s smooths are accounted for is still substantially predictable from the previous month’s residual, which an independent-errors model has no way to represent. That’s the assumption s(year_frac)’s edf-of-3.85 and its F-test were resting on, and it doesn’t hold here. The concern isn’t hypothetical: a smoothing penalty free to explain slow-moving, positively-autocorrelated noise as if it were curvature will do exactly that, because from the fitting procedure’s point of view, a real multi-year wiggle and a run of correlated residuals both look like “structure the flat line doesn’t capture.”
5 Refitting with an explicit correlation structure
5.1 Fitting the AR(1)-corrected model
fit_gamm <- gamm(log(result) ~ s(year_frac, bs = "tp", k = 10) + s(month, bs = "cc", k = 12),
data = do_wq, knots = list(month = c(0.5, 12.5)),
correlation = corAR1(form = ~1))
summary(fit_gamm$gam)
Family: gaussian
Link function: identity
Formula:
log(result) ~ s(year_frac, bs = "tp", k = 10) + s(month, bs = "cc",
k = 12)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.136947 0.008704 245.5 <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) 1.000 1 21.63 8.1e-06 ***
s(month) 7.494 10 72.87 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.916
Scale est. = 0.0030753 n = 1445.2 Confirming the correlation structure is earning its keep
phi <- coef(fit_gamm$lme$modelStruct$corStruct, unconstrained = FALSE)
phi Phi
0.5637352 fit_gamm0 <- gamm(log(result) ~ s(year_frac, bs = "tp", k = 10) + s(month, bs = "cc", k = 12),
data = do_wq, knots = list(month = c(0.5, 12.5)))
AIC(fit_gamm0$lme, fit_gamm$lme) df AIC
fit_gamm0$lme 5 -394.6952
fit_gamm$lme 6 -445.2578The estimated correlation is \(\hat\phi \approx\) 0.56, and the AR(1) model’s AIC is decisively lower than the no-correlation version’s — the correlation structure is earning its keep, not just adding a parameter for its own sake.
Once it’s in the model, s(year_frac)’s edf drops back to 1.00: a straight line. The wiggle the naive fit found on the trend term evaporates once the residual correlation is given somewhere else to go. The smooth term is still significant (\(p = 8.1\times10^{-6}\)) — there’s a real long-term decline here — but it’s a simple, constant-rate decline, not the more complex-looking trajectory the naive GAM implied.
5.3 Corrected rate per year
year_grid <- data.frame(year_frac = range(do_wq$year_frac), month = c(1, 1))
pred <- predict(fit_gamm$gam, newdata = year_grid, type = "terms", terms = "s(year_frac)")
slope_per_year <- diff(pred[, 1]) / diff(year_grid$year_frac)
slope_per_year 2
-0.01150398 (exp(slope_per_year) - 1) * 100 2
-1.143806 The corrected rate is about -1.14% per year.
6 The trend, with seasonality removed
A single slope-per-year number is a summary; seeing the trend itself is a better check on it. predict(..., type = "terms") isolates s(year_frac)’s contribution on its own, separate from s(month), which is exactly the “trend component” a classical decomposition (e.g. STL) would plot — recovered here as a term of the fitted GAMM rather than by a separate decomposition step.
trend_grid <- data.frame(
year_frac = seq(min(do_wq$year_frac), max(do_wq$year_frac), length.out = 200),
month = 6
)
trend_term <- predict(fit_gamm$gam, newdata = trend_grid, type = "terms",
terms = "s(year_frac)", se.fit = TRUE)
intercept <- coef(fit_gamm$gam)[["(Intercept)"]]
trend_grid <- trend_grid |>
mutate(
fit = intercept + trend_term$fit[, "s(year_frac)"],
lower = fit - 1.96 * trend_term$se.fit[, "s(year_frac)"],
upper = fit + 1.96 * trend_term$se.fit[, "s(year_frac)"]
)
ggplot(trend_grid, aes(year_frac, exp(fit))) +
geom_ribbon(aes(ymin = exp(lower), ymax = exp(upper)), alpha = 0.2) +
geom_line(linewidth = 1) +
labs(title = "Dissolved Oxygen — trend component, seasonality removed",
subtitle = "s(year_frac) term from the AR(1)-corrected GAMM, 95% CI",
y = "DO (mg/L)", x = "Year") +
nlt_theme
The band widens gradually toward both ends of the record, as fitted-curve uncertainty usually does away from the bulk of the data, but never crosses back to flat or rising — consistent with the significant, monotonic decline reported above, and a direct visual check on the edf = 1.00 finding: the line drawn here really is straight, not a curve that happens to average out to a negative slope.
7 Diagnostics
resid_norm <- as.numeric(resid(fit_gamm$lme, type = "normalized"))
shapiro.test(resid_norm)
Shapiro-Wilk normality test
data: resid_norm
W = 0.99221, p-value = 0.6189bptest(lm(resid_norm ~ fitted(fit_gamm$gam)))
studentized Breusch-Pagan test
data: lm(resid_norm ~ fitted(fit_gamm$gam))
BP = 0.8266, df = 1, p-value = 0.3633acf(resid_norm, main = "Residual autocorrelation, AR(1)-corrected")
The normalized residuals — the model’s own accounting for the AR(1) structure, standardized so they should look like independent noise if the correction worked — pass both tests (Shapiro-Wilk \(p = 0.619\), Breusch-Pagan \(p = 0.363\)), and the lag-1 autocorrelation that was above 0.5 in the naive fit is now essentially zero. The correction did what it was supposed to.
8 Interpretation
The naive GAM and the GAMM agree that Dissolved Oxygen is declining, but they disagree on how: ignoring residual autocorrelation made the long-term trend look like a genuinely wiggly, complex trajectory (edf 3.85); accounting for it collapsed the trend back down to a simple straight line (edf 1.00) declining at roughly 1.1% per year, still clearly significant once the correlation structure is given the noise it was actually explaining. That distinction isn’t cosmetic — a regulator or client reading “Dissolved Oxygen shows a complex, accelerating decline” would reasonably ask different questions than one reading “Dissolved Oxygen is declining steadily and slowly,” and the first version of that claim was an artifact of an assumption the naive fit never checked. The general lesson carries beyond this one series: fitting a GAM to a regularly sampled monitoring record and reading its smooth-term significance at face value skips a step: checking acf(resid(fit)) first, and reaching for gamm() with a correlation structure when it shows real autocorrelation, the way it did here.
9 Extending this analysis
corAR1() was the right correction here because a single lag-1 dependency was what the ACF plot actually showed, but it’s one member of a broader family nlme::lme() accepts — corARMA() for a richer autoregressive-moving-average structure, or a seasonal variant, would be worth checking if a residual ACF plot on a different record showed a shape AR(1) doesn’t match (a slower decay, or a spike at the seasonal lag instead of lag 1). And this remains a single-site model: extending it to a monitoring network of several sites would call for a hierarchical structure — site as a random effect on top of the same additive-smooth mean, so the trend and seasonal terms are estimated once across sites while still letting each site’s baseline differ, rather than refitting this same GAMM independently at every location.