Local regression from scratch

regression

nonparametric

ecology tutorial

Nadaraya-Watson and local linear regression in base R, bandwidth by cross-validation, and why a variance-only confidence band under-covers where the curve bends.

Author

Tidy Ecology

Published

2026-08-10

A Gaussian process fits a curve by putting a prior over functions. Local regression reaches a similar place from the opposite direction: to estimate the response at a point, just fit a simple model to the nearby data, weighting points by how close they are. No global functional form, no basis, only a kernel and a bandwidth. It is the idea behind LOESS, and it is worth building by hand because the machinery makes one honest lesson unavoidable: the usual pointwise confidence band is built from variance alone, and where the curve bends it is centred on a biased estimate, so its real coverage is well short of what the label claims.

Local averaging

The simplest local estimator is a weighted average. To predict the response at a target $x_0$, weight each observation by a kernel of its distance and take the weighted mean. This is the Nadaraya-Watson estimator, a local constant fit:

\[\hat m(x_0) = \frac{\sum_i K\!\big(\tfrac{x_i - x_0}{h}\big)\, y_i}{\sum_i K\!\big(\tfrac{x_i - x_0}{h}\big)}.\]

The bandwidth $h$ is the one knob that matters: small $h$ tracks the data closely but is noisy, large $h$ is smooth but blurs real structure. We use a Gaussian kernel.

set.seed(4219)
f_true <- function(x) 0.50*x +
  0.90*exp(-((x-0.30)^2)/(2*0.008)) +               # sharp bump (sd ~ 0.09)
  0.50*exp(-((x-0.70)^2)/(2*0.040))                 # broad bump (sd ~ 0.20)
sig_true <- 0.12
n  <- 120
x  <- sort(runif(n))
y  <- f_true(x) + rnorm(n, 0, sig_true)
xg <- seq(0, 1, length.out = 400)

gk <- function(u) exp(-0.5*u^2)                        # Gaussian kernel
nw <- function(x0, x, y, h){ w <- gk((x - x0)/h); sum(w*y)/sum(w) }   # local constant

From local constant to local linear

The weighted average has a weakness at the edges. Near a boundary the kernel is one-sided, so if the curve has a slope there, the average is pulled towards the interior and the fit is biased. The fix is to fit a local straight line instead of a local level: run a weighted least squares of $y$ on $x - x_0$ and read off the intercept. Writing the design as $X = [1, x - x_0]$ and the kernel weights as a diagonal $W$, the estimate is a linear combination of the responses, $\hat m(x_0) = e_1^\top (X^\top W X)^{-1} X^\top W\, y$, so we can extract the effective weight each observation carries. Collecting those weight vectors for every fitted point gives the smoother matrix $L$ with $\hat y = L y$, which we will need for both cross-validation and the confidence band.

ll_wts <- function(x0, x, h){                          # local linear effective weights
  w <- gk((x - x0)/h); X <- cbind(1, x - x0)
  XtW <- t(X * w); A <- XtW %*% X
  as.numeric((solve(A) %*% XtW)[1, ])
}
ll_fit <- function(xt, x, y, h){
  L <- t(vapply(xt, ll_wts, numeric(length(x)), x = x, h = h))
  list(mu = as.numeric(L %*% y), L = L)
}

Choosing the bandwidth

The smoother matrix gives leave-one-out cross-validation almost for free: the deleted residual at point $i$ is the ordinary residual divided by $1 - L_{ii}$, so we never refit. We minimise the cross-validation score over a grid of bandwidths.

cv_h <- function(h){
  L <- t(vapply(x, ll_wts, numeric(n), x = x, h = h))
  fit <- as.numeric(L %*% y); dL <- diag(L)
  mean(((y - fit)/(1 - dL))^2)
}
h_grid <- exp(seq(log(0.015), log(0.20), length.out = 40))
cv <- sapply(h_grid, cv_h)
h_cv <- h_grid[which.min(cv)]
cv_thr <- min(cv)*1.05; h_band <- h_grid[cv <= cv_thr]     # within 5% of the best score
h_lo <- min(h_band); h_hi <- max(h_band); h_ratio <- h_hi/h_lo

Cross-validation picks a bandwidth of 0.029, giving a fit with about 15 effective degrees of freedom. We fit at that bandwidth and attach the textbook variance-only band, using the residual variance and the row norms of the smoother matrix.

onx <- ll_fit(x, x, y, h_cv); L <- onx$L
nu  <- sum(diag(L))
sig2 <- sum((y - onx$mu)^2)/(n - nu)
fg  <- ll_fit(xg, x, y, h_cv); Lg <- fg$L
sd_g <- sqrt(sig2 * rowSums(Lg^2))                     # variance only, ignores bias
hw_g <- 1.96*sd_g; lo_g <- fg$mu - hw_g; hi_g <- fg$mu + hw_g
tru_g <- f_true(xg); ftx <- f_true(x)
inr <- xg >= min(x) & xg <= max(x)

Left: scattered points with a smooth green fit and a grey band; a dashed curve rises above the fit at a sharp peak near x = 0.3. Right: a U-shaped cross-validation curve with a wide flat trough shaded. — Figure 1: The fit and the bandwidth. (A) Local linear fit (green) with the variance-only 95 per cent band, the data, and the true curve (dashed). The fit sits below the true peak: the band is centred a little low there. (B) The cross-validation score against bandwidth; every bandwidth in the shaded band fits within five per cent of the best.

Why the band under-covers

The band above is honest about noise and silent about bias. To see the bias without any noise getting in the way, apply the smoother to the true mean function: the result is the fit you would get on average, and it differs from the truth by exactly the smoothing bias.

efit_g <- as.numeric(Lg %*% ftx)                       # expected fit = smoother on the truth
bias_g <- efit_g - tru_g                               # exact smoothing bias
ipk <- which.max(tru_g); xpk <- xg[ipk]                # the sharp optimum
bias_pk <- bias_g[ipk]; se_pk <- sd_g[ipk]; bias_in_se_pk <- abs(bias_pk)/se_pk
## true coverage of the NOMINAL 95% band: the fit is centred at truth + bias
cover_g  <- pnorm(1.96 - bias_g/sd_g) - pnorm(-1.96 - bias_g/sd_g)
cover_pk <- cover_g[ipk]; min_cov <- min(cover_g[inr])
frac_undercover <- mean(cover_g[inr] < 0.90)

At the sharp peak the fit sits 0.046 below the truth, which is 1.2 standard errors. The band is not too narrow; it is centred in the wrong place. A pointwise interval labelled 95 per cent therefore contains the true value at the peak only about 79 per cent of the time, and on the steep flank its true coverage falls to 62 per cent. Across the interior, 16 per cent of the grid has real coverage below 90 per cent. More data at this same cross-validated bandwidth would not fix it, because the missing piece is bias, not variance.

The bandwidth itself is barely pinned down. Bandwidths from 0.018 to 0.041, a 2.2-fold range, all fit within five per cent of the best cross-validation score, and across that range the fitted peak height swings from 1.04 to 1.13. How smooth the curve should be is a judgement the data leaves largely open.

Finally, the boundary is where the local constant estimator earns its bad name.

xb <- xg[xg <= 0.08]                                   # left edge, where the curve has a slope
nw_eb <- sapply(xb, function(x0){ w <- gk((x - x0)/h_cv); sum(w*ftx)/sum(w) })   # NW on the truth
Lb <- t(vapply(xb, ll_wts, numeric(n), x = x, h = h_cv)); ll_eb <- as.numeric(Lb %*% ftx)
bias_nw_0 <- abs(nw_eb[1] - f_true(xb[1])); bias_ll_0 <- abs(ll_eb[1] - f_true(xb[1]))

Right at the left edge, where the response is rising, the local constant fit is off by 0.013 while the local linear fit is off by only 0.003. A one-sided average cannot follow a trend; a local line can, which is why local linear is the sensible default.

Left: a coverage curve mostly near 0.95 but dipping to about 0.6 around x = 0.3. Right: a zoom at the left edge where a dotted curve separates upward from a dashed curve, and a green curve tracks the dashed one. — Figure 2: The two failures. (A) True coverage of the nominal 95 per cent band across the gradient; it drops well below the dashed 0.95 line wherever the curve bends. (B) Near the left boundary the local constant fit (dotted) peels away from the truth (dashed) while the local linear fit (green) stays with it.

Where this leaves you

Local regression is an easy, transparent smoother, and local linear is the version to reach for because it behaves at the edges. But read its pointwise band with the same caution as any plug-in interval. It accounts for variance and ignores the smoothing bias, so where the response curves sharply the interval is centred low and its true coverage falls well below the label, here to around 62 per cent. Undersmoothing on purpose shrinks the bias but widens the band and defeats the point; the honest response is to stop trusting a variance-only interval as if it were calibrated.

That is exactly what conformal prediction offers: a band around any of these fits, this one or the Gaussian process, with a finite-sample coverage guarantee that does not depend on getting the smoothness or the bias right. And whichever interval you use, the honest final step is to check its coverage on held-out data rather than trust the formula.

References

Cleveland WS 1979. Journal of the American Statistical Association 74(368):829-836 (10.1080/01621459.1979.10481038).

Cleveland WS, Devlin SJ 1988. Journal of the American Statistical Association 83(403):596-610 (10.1080/01621459.1988.10478639).

Fan J, Gijbels I 1996. Local Polynomial Modelling and Its Applications. Chapman and Hall/CRC. ISBN 978-0-412-98321-4.

Wasserman L 2006. All of Nonparametric Statistics. Springer, New York.

--- title: "Local regression from scratch" description: "Nadaraya-Watson and local linear regression in base R, bandwidth by cross-validation, and why a variance-only confidence band under-covers where the curve bends." date: "2026-08-10 10:00" categories: [R, regression, nonparametric, ecology tutorial] image: thumbnail.png image-alt: "A scatter of points with a smooth fitted curve, a shaded confidence band and a dashed true curve running above the fit at a sharp peak." --- A [Gaussian process](../gaussian-process-regression/) fits a curve by putting a prior over functions. Local regression reaches a similar place from the opposite direction: to estimate the response at a point, just fit a simple model to the nearby data, weighting points by how close they are. No global functional form, no basis, only a kernel and a bandwidth. It is the idea behind LOESS, and it is worth building by hand because the machinery makes one honest lesson unavoidable: the usual pointwise confidence band is built from variance alone, and where the curve bends it is centred on a biased estimate, so its real coverage is well short of what the label claims. ## Local averaging The simplest local estimator is a weighted average. To predict the response at a target $x_0$, weight each observation by a kernel of its distance and take the weighted mean. This is the Nadaraya-Watson estimator, a local constant fit: $$\hat m(x_0) = \frac{\sum_i K\!\big(\tfrac{x_i - x_0}{h}\big)\, y_i}{\sum_i K\!\big(\tfrac{x_i - x_0}{h}\big)}.$$ The bandwidth $h$ is the one knob that matters: small $h$ tracks the data closely but is noisy, large $h$ is smooth but blurs real structure. We use a Gaussian kernel. ```{r} #| label: dgp set.seed(4219) f_true <- function(x) 0.50*x + 0.90*exp(-((x-0.30)^2)/(2*0.008)) + # sharp bump (sd ~ 0.09) 0.50*exp(-((x-0.70)^2)/(2*0.040)) # broad bump (sd ~ 0.20) sig_true <- 0.12 n <- 120 x <- sort(runif(n)) y <- f_true(x) + rnorm(n, 0, sig_true) xg <- seq(0, 1, length.out = 400) gk <- function(u) exp(-0.5*u^2) # Gaussian kernel nw <- function(x0, x, y, h){ w <- gk((x - x0)/h); sum(w*y)/sum(w) } # local constant ``` ## From local constant to local linear The weighted average has a weakness at the edges. Near a boundary the kernel is one-sided, so if the curve has a slope there, the average is pulled towards the interior and the fit is biased. The fix is to fit a local straight line instead of a local level: run a weighted least squares of $y$ on $x - x_0$ and read off the intercept. Writing the design as $X = [1, x - x_0]$ and the kernel weights as a diagonal $W$, the estimate is a linear combination of the responses, $\hat m(x_0) = e_1^\top (X^\top W X)^{-1} X^\top W\, y$, so we can extract the effective weight each observation carries. Collecting those weight vectors for every fitted point gives the smoother matrix $L$ with $\hat y = L y$, which we will need for both cross-validation and the confidence band. ```{r} #| label: local-linear ll_wts <- function(x0, x, h){ # local linear effective weights w <- gk((x - x0)/h); X <- cbind(1, x - x0) XtW <- t(X * w); A <- XtW %*% X as.numeric((solve(A) %*% XtW)[1, ]) } ll_fit <- function(xt, x, y, h){ L <- t(vapply(xt, ll_wts, numeric(length(x)), x = x, h = h)) list(mu = as.numeric(L %*% y), L = L) } ``` ## Choosing the bandwidth The smoother matrix gives leave-one-out cross-validation almost for free: the deleted residual at point $i$ is the ordinary residual divided by $1 - L_{ii}$, so we never refit. We minimise the cross-validation score over a grid of bandwidths. ```{r} #| label: bandwidth cv_h <- function(h){ L <- t(vapply(x, ll_wts, numeric(n), x = x, h = h)) fit <- as.numeric(L %*% y); dL <- diag(L) mean(((y - fit)/(1 - dL))^2) } h_grid <- exp(seq(log(0.015), log(0.20), length.out = 40)) cv <- sapply(h_grid, cv_h) h_cv <- h_grid[which.min(cv)] cv_thr <- min(cv)*1.05; h_band <- h_grid[cv <= cv_thr] # within 5% of the best score h_lo <- min(h_band); h_hi <- max(h_band); h_ratio <- h_hi/h_lo ``` Cross-validation picks a bandwidth of `r sprintf("%.3f", h_cv)`, giving a fit with about `r sprintf("%.0f", sum(diag(ll_fit(x, x, y, h_cv)$L)))` effective degrees of freedom. We fit at that bandwidth and attach the textbook variance-only band, using the residual variance and the row norms of the smoother matrix. ```{r} #| label: fit-band onx <- ll_fit(x, x, y, h_cv); L <- onx$L nu <- sum(diag(L)) sig2 <- sum((y - onx$mu)^2)/(n - nu) fg <- ll_fit(xg, x, y, h_cv); Lg <- fg$L sd_g <- sqrt(sig2 * rowSums(Lg^2)) # variance only, ignores bias hw_g <- 1.96*sd_g; lo_g <- fg$mu - hw_g; hi_g <- fg$mu + hw_g tru_g <- f_true(xg); ftx <- f_true(x) inr <- xg >= min(x) & xg <= max(x) ``` ```{r} #| label: fig-fit #| echo: false #| fig-cap: "The fit and the bandwidth. (A) Local linear fit (green) with the variance-only 95 per cent band, the data, and the true curve (dashed). The fit sits below the true peak: the band is centred a little low there. (B) The cross-validation score against bandwidth; every bandwidth in the shaded band fits within five per cent of the best." #| fig-alt: "Left: scattered points with a smooth green fit and a grey band; a dashed curve rises above the fit at a sharp peak near x = 0.3. Right: a U-shaped cross-validation curve with a wide flat trough shaded." suppressPackageStartupMessages(library(ggplot2)) theme_te <- function() theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#dad9ca", linewidth = 0.3), axis.title = element_text(colour = "#2c3a31"), axis.text = element_text(colour = "#5d6b61"), plot.subtitle = element_text(colour = "#46604a", size = 10), plot.background = element_rect(fill = "white", colour = NA), panel.background = element_rect(fill = "white", colour = NA)) forest <- "#275139"; terra <- "#b5534e" dg <- data.frame(x = xg, mu = fg$mu, lo = lo_g, hi = hi_g, tru = tru_g) pA <- ggplot() + geom_point(data = data.frame(x, y), aes(x, y), colour = "#93a87f", alpha = 0.55, size = 1.4) + geom_ribbon(data = dg, aes(x, ymin = lo, ymax = hi), fill = "#93a87f", alpha = 0.30) + geom_line(data = dg, aes(x, tru), colour = "#16241d", linetype = "dashed", linewidth = 0.8) + geom_line(data = dg, aes(x, mu), colour = forest, linewidth = 1.1) + labs(x = "environmental gradient", y = "response", subtitle = "(A) local linear fit and 95% band") + theme_te() dcv <- data.frame(h = h_grid, cv = cv) pB <- ggplot(dcv, aes(h, cv)) + annotate("rect", xmin = h_lo, xmax = h_hi, ymin = -Inf, ymax = Inf, fill = "#93a87f", alpha = 0.30) + geom_line(colour = forest, linewidth = 1) + geom_vline(xintercept = h_cv, colour = "#5d6b61", linetype = "dotted") + labs(x = "bandwidth", y = "leave-one-out CV score", subtitle = "(B) the bandwidth is weakly pinned") + theme_te() grid::pushViewport(grid::viewport(layout = grid::grid.layout(1, 2))) print(pA, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1)) print(pB, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` ## Why the band under-covers The band above is honest about noise and silent about bias. To see the bias without any noise getting in the way, apply the smoother to the true mean function: the result is the fit you would get on average, and it differs from the truth by exactly the smoothing bias. ```{r} #| label: coverage efit_g <- as.numeric(Lg %*% ftx) # expected fit = smoother on the truth bias_g <- efit_g - tru_g # exact smoothing bias ipk <- which.max(tru_g); xpk <- xg[ipk] # the sharp optimum bias_pk <- bias_g[ipk]; se_pk <- sd_g[ipk]; bias_in_se_pk <- abs(bias_pk)/se_pk ## true coverage of the NOMINAL 95% band: the fit is centred at truth + bias cover_g <- pnorm(1.96 - bias_g/sd_g) - pnorm(-1.96 - bias_g/sd_g) cover_pk <- cover_g[ipk]; min_cov <- min(cover_g[inr]) frac_undercover <- mean(cover_g[inr] < 0.90) ``` At the sharp peak the fit sits `r sprintf("%.3f", abs(bias_pk))` below the truth, which is `r sprintf("%.1f", bias_in_se_pk)` standard errors. The band is not too narrow; it is centred in the wrong place. A pointwise interval labelled 95 per cent therefore contains the true value at the peak only about `r sprintf("%.0f", 100*cover_pk)` per cent of the time, and on the steep flank its true coverage falls to `r sprintf("%.0f", 100*min_cov)` per cent. Across the interior, `r sprintf("%.0f", 100*frac_undercover)` per cent of the grid has real coverage below 90 per cent. More data at this same cross-validated bandwidth would not fix it, because the missing piece is bias, not variance. The bandwidth itself is barely pinned down. Bandwidths from `r sprintf("%.3f", h_lo)` to `r sprintf("%.3f", h_hi)`, a `r sprintf("%.1f", h_ratio)`-fold range, all fit within five per cent of the best cross-validation score, and across that range the fitted peak height swings from `r sprintf("%.2f", ll_fit(xpk, x, y, h_hi)$mu)` to `r sprintf("%.2f", ll_fit(xpk, x, y, h_lo)$mu)`. How smooth the curve should be is a judgement the data leaves largely open. Finally, the boundary is where the local constant estimator earns its bad name. ```{r} #| label: boundary xb <- xg[xg <= 0.08] # left edge, where the curve has a slope nw_eb <- sapply(xb, function(x0){ w <- gk((x - x0)/h_cv); sum(w*ftx)/sum(w) }) # NW on the truth Lb <- t(vapply(xb, ll_wts, numeric(n), x = x, h = h_cv)); ll_eb <- as.numeric(Lb %*% ftx) bias_nw_0 <- abs(nw_eb[1] - f_true(xb[1])); bias_ll_0 <- abs(ll_eb[1] - f_true(xb[1])) ``` Right at the left edge, where the response is rising, the local constant fit is off by `r sprintf("%.3f", bias_nw_0)` while the local linear fit is off by only `r sprintf("%.3f", bias_ll_0)`. A one-sided average cannot follow a trend; a local line can, which is why local linear is the sensible default. ```{r} #| label: fig-coverage #| echo: false #| fig-cap: "The two failures. (A) True coverage of the nominal 95 per cent band across the gradient; it drops well below the dashed 0.95 line wherever the curve bends. (B) Near the left boundary the local constant fit (dotted) peels away from the truth (dashed) while the local linear fit (green) stays with it." #| fig-alt: "Left: a coverage curve mostly near 0.95 but dipping to about 0.6 around x = 0.3. Right: a zoom at the left edge where a dotted curve separates upward from a dashed curve, and a green curve tracks the dashed one." dcov <- data.frame(x = xg, cover = cover_g)[inr[order(xg)], ] qA <- ggplot(data.frame(x = xg[inr], cover = cover_g[inr]), aes(x, cover)) + geom_hline(yintercept = 0.95, colour = "#5d6b61", linetype = "dashed") + geom_line(colour = forest, linewidth = 1) + ylim(0.5, 1) + labs(x = "environmental gradient", y = "true coverage of the 95% band", subtitle = "(A) real coverage dips where the curve bends") + theme_te() xz <- seq(0, 0.15, length.out = 120) nwz <- sapply(xz, nw, x = x, y = y, h = h_cv); llz <- ll_fit(xz, x, y, h_cv)$mu dz <- rbind(data.frame(x = xz, val = f_true(xz), w = "truth"), data.frame(x = xz, val = nwz, w = "local constant"), data.frame(x = xz, val = llz, w = "local linear")) qB <- ggplot(dz, aes(x, val, colour = w, linetype = w)) + geom_line(linewidth = 0.9) + scale_colour_manual(values = c(truth = "#16241d", `local constant` = terra, `local linear` = forest)) + scale_linetype_manual(values = c(truth = "dashed", `local constant` = "dotted", `local linear` = "solid")) + labs(x = "environmental gradient", y = "response", colour = NULL, linetype = NULL, subtitle = "(B) local constant is biased at the edge") + theme_te() + theme(legend.position = c(0.68, 0.22), legend.background = element_rect(fill = "white", colour = "#dad9ca")) grid::pushViewport(grid::viewport(layout = grid::grid.layout(1, 2))) print(qA, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1)) print(qB, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` ## Where this leaves you Local regression is an easy, transparent smoother, and local linear is the version to reach for because it behaves at the edges. But read its pointwise band with the same caution as any plug-in interval. It accounts for variance and ignores the smoothing bias, so where the response curves sharply the interval is centred low and its true coverage falls well below the label, here to around `r sprintf("%.0f", 100*min_cov)` per cent. Undersmoothing on purpose shrinks the bias but widens the band and defeats the point; the honest response is to stop trusting a variance-only interval as if it were calibrated. That is exactly what [conformal prediction](../conformal-prediction-intervals/) offers: a band around any of these fits, this one or the [Gaussian process](../gaussian-process-regression/), with a finite-sample coverage guarantee that does not depend on getting the smoothness or the bias right. And whichever interval you use, the honest final step is to [check its coverage on held-out data](../checking-predictive-calibration/) rather than trust the formula. ## References Cleveland WS 1979. Journal of the American Statistical Association 74(368):829-836 (10.1080/01621459.1979.10481038). Cleveland WS, Devlin SJ 1988. Journal of the American Statistical Association 83(403):596-610 (10.1080/01621459.1988.10478639). Fan J, Gijbels I 1996. Local Polynomial Modelling and Its Applications. Chapman and Hall/CRC. ISBN 978-0-412-98321-4. Wasserman L 2006. All of Nonparametric Statistics. Springer, New York. ## Related tutorials - [Gaussian process regression from scratch](../gaussian-process-regression/) - [Conformal prediction intervals from scratch](../conformal-prediction-intervals/) - [Checking a predictive model: calibration and coverage](../checking-predictive-calibration/) - [Penalised regression splines from scratch](../penalised-regression-splines/)