LOESS
When \(\hat{x}_{i}\) is unevenly distributed across its range, a fixed bandwidth \(h\) leaves some local fits with many neighbors and others with very few.
LOESS (locally estimated scatterplot smoothing) fits a local linear (or quadratic) regression at each design point using an adaptive bandwidth: instead of fixing \(h\), it includes a fixed share of the data at each \(x\), set by the parameter span.
LOESS is useful for smoothing data with uneven \(X\): because the share of points in each neighborhood is fixed, the local fit has roughly the same effective sample size everywhere, even where the data are sparse. The span parameter plays the role \(h\) does for fixed-bandwidth methods, so smaller spans mean narrower neighborhoods, more flexibility, and noisier curves, while larger spans pool more data and smooth more aggressively.
Code
## Data
library(Ecdat)
xy2_hat <- Wages1[order(Wages1[, 'school']), c('school', 'wage')]
colnames(xy2_hat) <- c('x_hat', 'y_hat') # x_hat=years of schooling, y_hat=hourly wage
# Fit loess with the default span
reg_lo <- loess(y_hat ~ x_hat, data=xy2_hat)
plot(y_hat ~ x_hat, pch=16, col=grey(0, .1), data=xy2_hat,
main=NA, xlab='School', ylab='Wage')
lo_col <- adjustcolor('#20B2AA', alpha.f=.75) # teal green
lines(xy2_hat[, 'x_hat'], predict(reg_lo),
col=lo_col, type='o', pch=2)
legend('topleft',
legend=sprintf('Loess (span = %.2f)', reg_lo$pars$span),
lty=1, col=lo_col, cex=.8)
Bias-Variance Tradeoff
A better fit on the data you have is not automatically better, because a model can hit every point in the sample and still predict new observations poorly.
The bias-variance tradeoff describes two sources of model error: bias (systematic miss from over-smoothing) and variance (sample-to-sample wobble from under-smoothing).
The bias-variance tradeoff is useful for thinking about why a tuning parameter matters and how to set it. For local regressions, the bandwidth \(h\) (equivalently, the LOESS span) is the lever that governs the tradeoff: large \(h\) trades variance down for bias up, small \(h\) does the opposite. A loess fit narrow enough to interpolate every point has \(\hat{MSE}=0\) on the fit data but generalizes poorly because it has just copied the noise.
- Wide window (large span, large \(h\)): each local fit uses many observations, so the fitted curve barely changes from sample to sample (low variance). But a nearly flat fit over a wide region cannot follow a relationship that bends, so it is systematically off (high bias).
- Narrow window (small span, small \(h\)): each local fit uses few observations, so the fitted curve chases noise and swings from sample to sample (high variance). But a narrow window can track curvature closely (low bias).
Code
# Two-span loess with bootstrap confidence bands on one scatterplot
pred_design <- data.frame(x_hat=sort(unique(xy2_hat[, 'x_hat'])))
spans <- c(4.0, 0.4)
n_boot <- 399
span_col <- hcl.colors(3, alpha=.9)[c(1, 3)]
plot(y_hat ~ x_hat, data=xy2_hat, pch=16, col=grey(0, .1),
main=NA, xlab='School', ylab='Wage')
for (s in seq_along(spans)) {
sp <- spans[s]
reg_s <- loess(y_hat ~ x_hat, data=xy2_hat, span=sp)
boot_preds <- matrix(NA, nrow=nrow(pred_design), ncol=n_boot)
for (b in 1:n_boot) {
xy2_boot <- xy2_hat[sample(nrow(xy2_hat), replace=TRUE), ]
reg_s_b <- loess(y_hat ~ x_hat, data=xy2_boot, span=sp)
boot_preds[, b] <- predict(reg_s_b, newdata=pred_design)
}
boot_cb <- apply(boot_preds, 1, quantile,
probs=c(.025, .975), na.rm=TRUE)
polygon(c(pred_design[, 'x_hat'], rev(pred_design[, 'x_hat'])),
c(boot_cb[1, ], rev(boot_cb[2, ])),
col=adjustcolor(span_col[s], alpha.f=0.3), border=NA)
lines(pred_design[, 'x_hat'], predict(reg_s, newdata=pred_design),
col=span_col[s], lwd=2)
}
legend('topleft', title='LOESS span',
legend=paste0(spans),
col=span_col, lty=1, lwd=2, cex=.8)
With span \(4.0\) the fitted curve is very smooth and stable, and the bootstrap band sits tightly around it. With span \(0.4\) the curve tracks more local variation, and the bootstrap band fans wider, directly displaying the variance penalty for chasing local detail. (Schooling x_hat is integer-valued with only \(14\) distinct values, so spans below about \(0.4\) collapse to single-point neighborhoods and loess() cannot fit.)
The bootstrap bands show variance, but not bias, because we never see the true curve in one sample. As in Sampling, pretend the \(3294\) workers in Wages1 are the whole population, so the true curve is known. The population curve is the mean wage at each schooling level over all \(3294\) workers. We then draw many samples of \(n=300\) workers, fit a wide span (\(4.0\)) and a narrow span (\(0.5\)) to each, and see how the fitted curves compare to the population curve. At \(n=300\), span \(0.4\) sometimes fails, so we use \(0.5\) for the narrow window.
Code
# Pretend these 3294 workers are the whole population
data(Wages1, package='Ecdat')
population_xy <- Wages1[, c('school', 'wage')]
colnames(population_xy) <- c('x', 'y') # x=years of schooling, y=hourly wage
population_slope <- coef(lm(y ~ x, data=population_xy))[2]
# Population curve: mean wage at each schooling level
population_curve <- tapply(population_xy[, 'y'], population_xy[, 'x'], mean)
# Evaluate at schooling levels with many workers (8 to 15 years)
# Rare levels are often missing from a sample, and loess does not extrapolate
x_levels <- 8:15
population_levels <- population_curve[as.character(x_levels)]
# Repeated samples, each fit with a wide and a narrow span
spans_sim <- c(4.0, 0.5)
n <- 300
n_rep <- 200
curves_sim <- array(NA, dim=c(length(x_levels), n_rep, length(spans_sim)))
for (r in 1:n_rep) {
s_id <- sample(nrow(population_xy), n, replace=FALSE)
xy_sim <- population_xy[s_id, ]
colnames(xy_sim) <- c('x_hat', 'y_hat')
for (s in seq_along(spans_sim)) {
# Return NA if loess fails on this sample
curves_sim[, r, s] <- tryCatch(
suppressWarnings(predict(
loess(y_hat ~ x_hat, data=xy_sim, span=spans_sim[s]),
newdata=data.frame(x_hat=x_levels))),
error=function(e) rep(NA, length(x_levels)))
}
}
# Plot a few sample curves for each span, plus the population curve
par(mfrow=c(1, 2))
for (s in seq_along(spans_sim)) {
matplot(x_levels, curves_sim[, 1:20, s], type='l', lty=1,
col=adjustcolor(span_col[s], alpha.f=.4),
ylim=range(curves_sim[, 1:20, ], population_levels, na.rm=TRUE),
xlab='School', ylab='Wage', main=NA)
lines(x_levels, population_levels, col=rgb(1, 0, 0, .8), lwd=2, type='o', pch=16)
title(paste0('Span = ', spans_sim[s]), font.main=1)
}
Each thin line is the fit from one sample, and the red line is the population curve. The wide-span curves stay close to one another, but they are too straight to follow the bend in the population curve. The narrow-span curves follow the bend better on average, but they spread out much more from sample to sample.
We can summarize this with two numbers per span. At each schooling level, the squared bias is the squared gap between the average fitted value and the population curve, and the variance is the spread of the fitted values across samples. We then average each over the schooling levels.
Code
# Average squared bias and variance across schooling levels
bias_var <- sapply(seq_along(spans_sim), function(s) {
curves_s <- curves_sim[, , s]
bias2_s <- (rowMeans(curves_s, na.rm=TRUE) - population_levels)^2
var_s <- apply(curves_s, 1, var, na.rm=TRUE)
c(Bias2=mean(bias2_s), Variance=mean(var_s), MSE=mean(bias2_s + var_s))
})
colnames(bias_var) <- paste0('Span ', spans_sim)
round(bias_var, 3)
## Span 4 Span 0.5
## Bias2 0.034 0.009
## Variance 0.128 0.301
## MSE 0.162 0.310
The wide span typically has the larger squared bias and the smaller variance, while the narrow span has the opposite. Here the variance term is the larger of the two, so the wide span often has the lower total error at \(n=300\). With a larger sample, the variance of both fits falls, and the narrow span becomes more attractive. Choosing the span (the bandwidth \(h\)) means trading bias against variance.
Note that when the local fits are linear (degree=1), the simple linear regression is nested as a special case with \(span \to \infty\), because every observation then gets nearly the same weight. By default, loess() fits local quadratics (degree=2), so a very wide default fit approaches a global quadratic instead of a straight line.
Model Fit
Each local model produces fitted values \(\hat{y}_{i}^{\mathrm{fit}}\) and residuals \(\hat{e}_{i} = \hat{y}_{i} - \hat{y}_{i}^{\mathrm{fit}}\). To compare models, we often summarize the residuals with a single number.
The mean squared error averages the squared residuals, \[\begin{aligned}
\hat{MSE} &= \frac{1}{n}\sum_{i=1}^{n} \hat{e}_{i}^2 = \frac{1}{n}\sum_{i=1}^{n} \left(\hat{y}_{i} - \hat{y}_{i}^{\mathrm{fit}}\right)^2.
\end{aligned}\] Squaring penalizes large misses heavily, and the units are the square of \(\hat{y}_{i}\) (squared dollars for wages).
The mean absolute percentage error instead averages each residual relative to its observed value, \[\begin{aligned}
\hat{MAPE} &= \frac{100}{n}\sum_{i=1}^{n} \left| \frac{\hat{e}_{i}}{\hat{y}_{i}} \right|.
\end{aligned}\] MAPE is unit-free, which makes it easy to communicate, but it is undefined when \(\hat{y}_{i}=0\) and over-weights observations with small \(\hat{y}_{i}\). MSE is easier to analyze theoretically (formally deriving the bias-variance tradeoff)
Suppose three workers have wages \(\hat{y}=(8, 10, 15)\) and a model predicts \(\hat{y}^{\mathrm{fit}}=(9, 9, 14)\). The residuals are \(\hat{e}=(-1, 1, 1)\), so \[\begin{aligned}
\hat{MSE} &= \frac{(-1)^2 + 1^2 + 1^2}{3} = \frac{3}{3} = 1, \\
\hat{MAPE} &= \frac{100}{3}\left( \frac{1}{8} + \frac{1}{10} + \frac{1}{15} \right) \approx \frac{100}{3}(0.292) \approx 9.7.
\end{aligned}\] The model misses by \(1\) squared-dollar on average, or about \(10\%\) of each wage.
Code
y0_hat <- c(8, 10, 15)
y0_fit <- c(9, 9, 14)
e0_hat <- y0_hat - y0_fit
mean(e0_hat^2) # MSE
## [1] 1
100*mean(abs(e0_hat/y0_hat)) # MAPE
## [1] 9.722222
We can compute these for the linear, regressogram, and piecewise models on the wage~school data.
Code
# xy2_hat defined above
# Three models for wage on schooling
reg_lin <- lm(y_hat ~ x_hat, data=xy2_hat) # globally linear
xy2_hat[, 'xcf'] <- cut(xy2_hat[, 'x_hat'], 3)
reg_rgram <- lm(y_hat ~ xcf, data=xy2_hat) # regressogram, 3 bins
reg_pw <- lm(y_hat ~ xcf*x_hat, data=xy2_hat) # piecewise, 3 bins
# Mean squared error and mean absolute percentage error
fit_stats <- function(model){
e_hat <- resid(model)
y_hat <- xy2_hat[, 'y_hat']
c(MSE=mean(e_hat^2), MAPE=100*mean(abs(e_hat/y_hat)))
}
round(rbind(
Linear = fit_stats(reg_lin),
Regressogram = fit_stats(reg_rgram),
Piecewise = fit_stats(reg_pw)), 2)
## MSE MAPE
## Linear 9.83 73.27
## Regressogram 10.27 75.88
## Piecewise 9.73 72.82
