k-Nearest Neighbors Regression

Total Phosphorus — a rise that levels into a plateau

1 A purely local rule

Total Phosphorus rises and then settles onto a plateau — not the slope break of piecewise regression, not the level shift of changepoint detection, and no natural single knot an analyst would place by eye. k-nearest-neighbors regression sidesteps the question of a functional form entirely: for a query point \(x_0\) on year_frac, let \(N_k(x_0)\) be the set of \(k\) training observations closest to it (Euclidean distance on year_frac), and predict

\[\hat y(x_0) = \frac{1}{k} \sum_{i \in N_k(x_0)} y_i,\]

the plain average of their responses. There’s no linear predictor at all here, not even the feature-engineered kind covered in The Regression Framework — no coefficients are estimated, just a lookup-and-average rule applied at prediction time.

Concentration data is modeled on the log scale throughout this site for the reasons covered in Transformations — a multiplicative process where variability scales with the mean — and Total Phosphorus is no exception; log(result) is used below.

2 The shape in the record

ggplot(tp_wq, aes(year_frac, result)) +
  geom_point(alpha = 0.5) +
  geom_smooth(se = FALSE) +
  labs(title = "Total Phosphorus — raw scale", y = "Phosphorus (mg/L)") +
  nlt_theme

An early low, roughly flat stretch, a rise through the middle years, and a plateau at a higher level for the remainder of the record — the shape a logistic curve would produce, though nothing here assumes that form.

3 Choosing k: the bias-variance tradeoff

\(k\) is the one tuning parameter this method has, and it controls a direct bias-variance tradeoff. A small \(k\) averages over only a handful of nearby points, so the fitted curve tracks the training data’s noise as much as its signal — low bias, high variance. A large \(k\) averages over a wide neighborhood, smoothing away real curvature along with the noise — low variance, high bias. There’s no fitted coefficient to inspect here, so the tradeoff has to be shown, not asserted:

grid <- data.frame(year_frac = seq(min(tp_wq$year_frac), max(tp_wq$year_frac), length.out = 300))

sweep_ks <- c(3, 20, 60)
sweep_fits <- lapply(sweep_ks, function(k) {
  m <- knnreg(log(result) ~ year_frac, data = tp_wq, k = k)
  grid |> mutate(fit = exp(predict(m, newdata = grid)), k = factor(k))
})
sweep_df <- bind_rows(sweep_fits)

ggplot() +
  geom_point(data = tp_wq, aes(year_frac, result), alpha = 0.35, color = "grey40") +
  geom_line(data = sweep_df, aes(year_frac, fit, color = k), linewidth = 1) +
  labs(title = "Total Phosphorus — k-NN fit at three values of k", y = "Phosphorus (mg/L)", color = "k") +
  nlt_theme +
  theme(legend.position = "right")

At \(k = 3\) the curve is jagged, chasing individual points rather than the underlying rise; at \(k = 60\) it’s smoothed almost flat, blunting the plateau’s onset; \(k = 20\) sits between the two, tracking the rise without visibly reacting to single points.

Rather than pick a value by eye, \(k\) is chosen by 10-fold cross-validation, minimizing mean squared error on the log scale across held-out folds:

set.seed(1)
folds <- createFolds(tp_wq$result, k = 10)
candidate_ks <- c(3, 5, 10, 15, 20, 30, 45, 60)

cv_mse <- sapply(candidate_ks, function(k) {
  fold_errs <- sapply(folds, function(idx) {
    train <- tp_wq[-idx, ]
    test <- tp_wq[idx, ]
    m <- knnreg(log(result) ~ year_frac, data = train, k = k)
    mean((log(test$result) - predict(m, newdata = test))^2)
  })
  mean(fold_errs)
})

data.frame(k = candidate_ks, cv_mse = round(cv_mse, 4))

   k cv_mse
1  3 0.0239
2  5 0.0242
3 10 0.0208
4 15 0.0201
5 20 0.0203
6 30 0.0207
7 45 0.0227
8 60 0.0327

best_k <- candidate_ks[which.min(cv_mse)]
fit_tp <- knnreg(log(result) ~ year_frac, data = tp_wq, k = best_k)
best_k

[1] 15

Cross-validated error is minimized at \(k =\) 15, consistent with the sweep above: small enough to track the rise-and-plateau shape, large enough that no single observation swings the prediction.

4 The fitted curve

grid <- grid |> mutate(fit = exp(predict(fit_tp, newdata = grid)))

ggplot() +
  geom_point(data = tp_wq, aes(year_frac, result), alpha = 0.4, color = "grey40") +
  geom_line(data = grid, aes(year_frac, fit), color = "steelblue", linewidth = 1) +
  labs(title = sprintf("Total Phosphorus — k-NN fit (k = %d)", best_k), y = "Phosphorus (mg/L)") +
  nlt_theme

There’s no confidence band here, and none is available: knnreg() produces a fitted value at each point, not a coefficient with a standard error. Nothing about this fit supports a significance test, a slope estimate, or a p-value — it’s a prediction-oriented baseline, not an inferential model. The cross-validated MSE above (0.0201 on the log scale) is the one number this method offers in place of the usual diagnostics, and it answers “how well does this rule predict held-out points,” not “is there a real trend here.”

5 Interpretation

k-NN recovers the rise-and-plateau shape without being told any functional form, using only a neighborhood-averaging rule tuned by cross-validation. That’s useful for showing what the record looks like, but it stops there: no coefficient says how fast Phosphorus rose, no interval bounds where the plateau sits, and no test says whether the apparent leveling-off is real or just where the data happens to thin out. The next section fits LOESS to the same data — a locally weighted regression in place of k-NN’s local average, still without a global functional form, but a step closer to the kind of fit a formal model can be built on.

6 Extending this analysis

A rise-then-plateau shape like this one is exactly what a logistic growth curve describes, and k-NN’s fitted curve is a reasonable read of that from the raw data alone. If that read holds up, the better next move for this specific shape isn’t necessarily the fully nonparametric route this site takes through LOESS and a GAM — it’s fitting the logistic form directly with nls(), which would recover the inferential machinery k-NN structurally can’t provide (a rate parameter, an inflection point, confidence intervals on both) at a fraction of the parameters a smooth-basis fit needs. That’s a genuine tradeoff against the GAM route taken here, not a strictly better option: a parametric curve fit this way is only as good as the assumption that the shape really is logistic, which the GAM never has to commit to.