This chapter asks three questions.
- How does the average of \(\hat{y}_{i}\) change with \(\hat{x}_{i}\) without assuming a single straight line?
- Regressogram and piecewise regression: these split \(\hat{x}\) into fixed bins and fit a constant or a line in each bin. Both are new, and each bin runs a simple regression on its own observations.
- How do we avoid jumps at the bin edges?
- Locally constant and local linear regression: these fit a constant or a line in a window that slides along \(x\). Both are new, and each window uses the uniform kernel from the kernel density.
- How local should each estimate be?
- Bandwidth: this half-width \(h\) is the same one that set the bins of a histogram. Larger \(h\) pools more observations, and smaller \(h\) tracks the data more closely.
Local Relationships
Scatterplots are a great and simple plot for bivariate data that simply plots each observation. There are many extensions and similar tools. Many of them summarize \(\hat{y}_{i}\) using only the observations whose \(\hat{x}_{i}\) is near some value \(x\).
An observation is “local” to \(x\) when \(|\hat{x}_{i}-x|\leq h\), for some bandwidth \(h\).
The bandwidth \(h\) is the same \(h\) that set the bin half-width of a histogram and the spread of a kernel density. It is useful for controlling how smooth a local estimate is.
The example below uses a small bandwidth for local means and medians to help understand how the central tendency changes. Each value of \(\hat{x}\) is a group, so the red triangles are the conditional means \(\hat{m}_{Y}(g)\) from Comparing Multiple Groups, now ordered along \(\hat{x}\). This chapter writes means \(\hat{m}_{Y}(g)\) as \(\hat{m}_{Y}(x,h)\), because each group is the set of observations with \(\hat{x}_{i}~\text{in}~[x-h,x+h]\). The bandwidth \(h\) is small enough here that each group holds only that value, but we could expand it to pool more data. The histogram bandwidth \(h\) controlled how much data went into a density estimate at \(x\) whereas here \(h\) controls how much data goes into a conditional estimate at \(x\).
Overlayed on a scatterplot, we can also get a sense of how dispersion changes too. The blue squares show the grouped medians, the conditional quantile \(\hat{F}_{Y}^{-1}(0.5 \mid x)\) at each value of \(\hat{x}\). A median below its mean is a sign that \(\hat{y}\) is skewed right within that group, so we can also get a sense of whether skew is changing with \(x\).
Code
## MAIN EXAMPLE
# Running example for the whole chapter: wages and years of schooling
data(Wages1, package='Ecdat')
xy_hat <- Wages1[order(Wages1[, 'school']), c('school', 'wage')]
colnames(xy_hat) <- c('x_hat', 'y_hat') # x_hat=years of schooling, y_hat=hourly wage
# Scatter of every observation
plot(y_hat ~ x_hat, data=xy_hat, pch=16, col=grey(0, .1),
main=NA, xlab='School', ylab='Wage')
# Conditional mean at each integer schooling level
# Each group is the window of half-width 1/2 around one schooling level
school_means <- aggregate(y_hat ~ x_hat, data=xy_hat, mean)
points(y_hat ~ x_hat, data=school_means, # connect dots so the trend pops
pch=17, col=rgb(1, 0, 0, .8), type='b')
# Conditional median at each integer schooling level
school_medians <- aggregate(y_hat ~ x_hat, data=xy_hat, median)
points(y_hat ~ x_hat, data=school_medians,
pch=15, col=rgb(0, 0, 1, .8), type='b')
legend('topleft', legend=c('Mean', 'Median'),
pch=c(17, 15), lty=1, col=c(rgb(1, 0, 0, .8), rgb(0, 0, 1, .8)), bty='n')
title('Grouped Means, Medians, and Scatterplot', font.main=1)
Code
# Nearly all workers have 7 to 15 years of schooling.
# At each of those levels, the median sits below the mean.
# That is a sign that wages within each level are skewed right.
Regression models provide a general way to describe how the mean \(\hat{y}_{i}\) changes with \(\hat{x}_{i}\). A simple linear regression fits one straight line to all of the data, which can miss a relationship that bends. The methods in this chapter instead fit a simple model to the observations within a window of \(\hat{x}_{i}\) values. Each method splits the outcome into two parts, \[
\hat{y}_{i} = \hat{y}_{i}^{\mathrm{fit}} + \hat{e}_{i},
\] where the fitted value \(\hat{y}_{i}^{\mathrm{fit}}\) can bend with \(\hat{x}_{i}\) instead of following a single straight line, and \(\hat{e}_{i}\) is the residual. After minimizing the sum of squared errors in the regressions below, our estimated models will also predict local averages.
The methods differ in two choices: whether the windows are fixed bins or slide along with \(x\), and whether the model within each window is a constant or a line. All four methods in the table use a bandwidth \(h\), which controls how “local” a local estimate actually is.
Local Regression Methods
| Fixed bins |
Regressogram |
Piecewise regression |
| Window slides with \(x\) |
Locally constant regression |
Local linear regression |
Each step fixes one problem with the step before it. Adding a slope captures the trend within a bin. Sliding the window replaces the jump at the bin edge with smaller steps. Adding a slope to the sliding window reduces bias near the edges of the data. A window near an edge can hold many observations from one side of the design point, which pulls the window average away from the trend. The local linear fit follows the trend. The code below checks this at one design point near an edge.
Code
# The figure above uses half-width h=3.5 and two fixed bins.
# The bins hold schooling levels 3 to 9 and 10 to 16, with midpoints 6 and 13.
# The sliding window around 15 years of schooling, near the top of the range
xy_window <- xy_hat[abs(xy_hat[, 'x_hat']-15) <= 3.5, ]
table(xy_window[, 'x_hat'])
##
## 12 13 14 15 16
## 1188 339 264 134 16
# Of its 1941 workers, 1188 have 12 years of schooling.
# Grouped mean at 15 years, then the locally constant and local linear fits
school_means[school_means[, 'x_hat']==15, 'y_hat']
## [1] 8.375111
mean(xy_window[, 'y_hat'])
## [1] 6.291988
predict(lm(y_hat ~ x_hat, data=xy_window), newdata=data.frame(x_hat=15))
## 1
## 8.244139
# Near the top, the locally constant fit falls below the grouped means.
# At 15 years it is about 6.29, and the grouped mean is about 8.38.
# The many workers with 12 years of schooling pull the window average down.
# The local linear fit, about 8.24, follows the trend.
The figure below fits all four methods to the running example, each with the same half-width \(h\). The left panel uses two fixed bins, and the right panel slides a window of the same width along \(x\). The grey triangles are the grouped means from the figure above. The dotted lines in the left panel mark the bin midpoints, and the circles mark the predictions there. At these midpoints, each window holds exactly one bin, so the sliding fits equal the fixed-bin fits there.
Fixed Bins
The methods in the first row of the table split \(\hat{x}\) into fixed bins and fit a separate model within each bin.
Regressograms
Just as a histogram summarizes how observations are distributed across bins of \(\hat{x}\), a regressogram summarizes how the conditional mean of \(\hat{y}_{i}\) changes across bins of \(\hat{x}\).
A regressogram splits \(\hat{x}\) into exclusive intervals (bins) of half-width \(h\) and reports the average of \(\hat{y}_{i}\) within each bin. With bin midpoints \(x_{1}, ..., x_{L}\), the model is \[\begin{aligned}
\hat{y}_{i} &= \sum_{x~\text{in}~\{x_{1}, ..., x_{L} \}} b_{0}(x,h) \hat{d}_{i}(x,h) + \hat{e}_{i},
\end{aligned}\] where \(\hat{d}_{i}(x,h) = \mathbf{1}\left(\hat{x}_{i}~\text{in}~\left(x-h,x+h\right] \right)\) is a dummy variable for the bin with midpoint \(x\), and each bin has a coefficient \(b_{0}(x,h)\). The model is fit using ordinary least squares (OLS).
The regressogram is useful when we expect the relationship between \(\hat{x}\) and \(\hat{y}\) to bend in ways a single line cannot capture, because the bin averages can rise, fall, or stay flat from one bin to the next without imposing any global shape.
Code
## MAIN EXAMPLE
# Data
plot(y_hat ~ x_hat, data=xy_hat, pch=16, col=grey(0, .1),
main=NA, xlab='School', ylab='Wage')
# xy_hat defined above
## Simple Regression
reg <- lm(y_hat ~ x_hat, data=xy_hat) ## OLS
# Regressogram: Two Schooling Bins, levels 3-9 and 10-16
xy_hat[, 'xcc'] <- cut(xy_hat[, 'x_hat'], breaks=c(2.5, 9.5, 16.5))
rgram_c <- lm(y_hat ~ xcc, data=xy_hat)
# Regressogram: Three Schooling Bins
xy_hat[, 'xcf'] <- cut(xy_hat[, 'x_hat'], 3)
rgram_f <- lm(y_hat ~ xcf, data=xy_hat)
## Compare Models (Only 2 for simplicity)
lines( xy_hat[, 'x_hat'], predict(reg), lwd=2, col=rgb(1, 0, 0, .8))
lines( xy_hat[, 'x_hat'], predict(rgram_f), lwd=2, col=rgb(0, 0, 1, .8), type='s')
legend('topleft',
legend=c('Linear Regression', 'Regressogram (3)'),
col=c(rgb(1, 0, 0, .8), rgb(0, 0, 1, .8)),
lty=1, lwd=2, cex=.8)
In the regressogram model, \(\hat{e}_{i}\) depends on the trial bin coefficients \(b_{0}(x,h)\); the final residuals use the fitted coefficients \(\hat{b}_{0}(x,h)\).
When minimizing the sum of squared errors, the optimal coefficients are denoted as \(\hat{b}_{0}(x,h)\) and we can find them analytically. To do so, notice that each bin has \(n(x,h) = \sum_{i}^{n}\hat{d}_{i}(x,h)\) observations. This means we can split the dataset into parts associated with each bin \[\begin{aligned}
\sum_{i}^{n}\left[\hat{e}_{i}\right]^2
&= \sum_{i}^{n}\left[\hat{y}_{i}- \sum_{x~\text{in}~\{x_{1}, ..., x_{L} \}} b_{0}(x,h) \hat{d}_{i}(x,h) \right]^2 \\
&= \sum_{i}^{n(x_{1},h)}\left[\hat{y}_{i}- \sum_{x~\text{in}~\{x_{1}, ..., x_{L} \}} b_{0}(x,h) \hat{d}_{i}(x,h) \right]^2
+ ... \sum_{i}^{n(x_{L},h)}\left[\hat{y}_{i}- \sum_{x~\text{in}~\{x_{1}, ..., x_{L} \}} b_{0}(x,h) \hat{d}_{i}(x,h) \right]^2 \\
&= \sum_{i}^{n(x_{1},h)}\left[\hat{y}_{i}- b_{0}\left(x_{1},h\right) \right]^2 + ... \sum_{i}^{n(x_{L},h)}\left[\hat{y}_{i}-b_{0}\left(x_{L},h\right) \right]^2 % +~ (N-1)\sum_{i}\hat{y}_{i}.
\end{aligned}\] This separation allows us to analytically optimize for each bin separately \[\begin{aligned}
\min_{ \left\{ b_{0}(x,h) \right\} } \sum_{i}^{n}\left[\hat{e}_{i}\right]^2
&= \sum_{x~\text{in}~\{x_{1}, ..., x_{L} \}} \min_{ b_{0}(x,h) } \sum_{i}^{n(x,h)}\left[\hat{y}_{i}- b_{0}\left(x,h\right) \right]^2.
\end{aligned}\] Within each bin, this is the simple constant regression applied to the observations in that bin, so minimizing yields the optimal coefficient as follows \[\begin{aligned}
0 &= -2 \sum_{i}^{n(x,h)}\left[ \hat{y}_{i} - b_{0}(x,h) \right] \\
\hat{b}_{0}(x,h) &= \frac{\sum_{i}^{n(x,h)} \hat{y}_{i}}{ n(x,h) } = \hat{m}_{Y}(x,h).
\end{aligned}\] As such, the OLS regression yields coefficients that are interpreted as the conditional mean: \(\hat{m}_{Y}(x,h)\). We can directly compute the same statistic by simply taking the average value of \(\hat{y}_{i}\) for all \(i\) observations in a particular bin.
In any case, the values predicted by the model are found as \[
\hat{y}_{i}^{\mathrm{fit}} = \sum_{x} \hat{b}_{0}(x,h) \hat{d}_{i}(x,h) = \sum_{x} \hat{m}_{Y}(x,h) \hat{d}_{i}(x,h).
\] I.e., the predicted value of observation \(i\) is the average value for its bin.
Consider this two-bin regressogram example of how \(\hat{y}_{i}\) differs between a lower and an upper bin of \(\hat{x}_{i}\). \[\begin{aligned}
\hat{y}_{i} &= b_{0}(x=6, h=7/2) \mathbf{1}\left(\hat{x}_{i}~\text{in}~(2.5,9.5]\right) + b_{0}(x=13, h=7/2) \mathbf{1}\left(\hat{x}_{i}~\text{in}~(9.5,16.5] \right) + \hat{e}_{i}.
\end{aligned}\]
Here is a simple example with three data points \((\hat{x}_{i}, \hat{y}_{i})~\text{in}~\{ (3,2), (4,4), (12,3) \}\), for which we can easily compute \[\begin{aligned}
n(x=6,h=7/2) &= 2 \\
\hat{b}_{0}(x=6, h=7/2) &= [2 + 4] / 2 \\
n(x=13,h=7/2) &= 1 \\
\hat{b}_{0}(x=13, h=7/2) &= 3 / 1
\end{aligned}\]
Here is a simple example with data
Code
# Two bins: workers with at most 9 years of schooling and workers with more
y_fit_c <- predict(rgram_c)
fit_dat <- data.frame(xcc=xy_hat[, 'xcc'], y_fit=round(y_fit_c, 6))
table(fit_dat)
## y_fit
## xcc 4.115873 5.917872
## (2.5,9.5] 293 0
## (9.5,16.5] 0 3001
## Compare to simple aggregation
aggregate(y_hat ~ xcc, data=xy_hat, mean)
## xcc y_hat
## 1 (2.5,9.5] 4.115873
## 2 (9.5,16.5] 5.917872
# The mean wage is about 4.12 dollars in the first bin and 5.92 in the second.
Piecewise Regression
To let the slope vary across the range of \(\hat{x}\) instead of just the level, we add a slope term within each bin.
A piecewise regression (also segmented regression) divides \(\hat{x}\) into bins and runs a separate simple linear regression on each bin. The model is \[\begin{aligned}
\hat{y}_{i} &= \sum_{x} \left[b_{0}(x,h) + b_{1}(x,h)\hat{x}_{i} \right] \hat{d}_{i}(x,h) + \hat{e}_{i},
\end{aligned}\] where \(\hat{d}_{i}(x,h)\) is the dummy variable for the bin with midpoint \(x\), and each bin has its own intercept \(b_{0}(x,h)\) and slope \(b_{1}(x,h)\). The model is fit using OLS. Within a bin the relationship is a straight line; across bins the slope and intercept can change.
Piecewise regression is useful when the relationship is approximately linear in pieces but the slope itself changes. Returns to schooling might flatten after a degree, for example, or marginal cost curves might kink at capacity. A regressogram only captures the average within each bin, and a single line has the same slope everywhere. Piecewise regression captures the trend within each bin and lets the slope change at bin boundaries. A bin boundary can also mark a real change in the data, such as a new political regime.
Code
## MAIN EXAMPLE
# Data
plot(y_hat ~ x_hat, data=xy_hat, pch=16, col=grey(0, .1),
main=NA, xlab='School', ylab='Wage')
# xy_hat, with the bins xcc and xcf, defined above
# Piecewise: Two Schooling Bins
preg_c <- lm(y_hat ~ xcc*x_hat, data=xy_hat)
# Piecewise: Three Schooling Bins
preg_f <- lm(y_hat ~ xcf*x_hat, data=xy_hat)
## Compare Models
lines( xy_hat[, 'x_hat'], predict(preg_c), lwd=2, col=rgb(1, 0, 0, .8))
lines( xy_hat[, 'x_hat'], predict(preg_f), lwd=2, col=rgb(1, 0, 0, .8), lty=2)
legend('topleft',
legend=c('2 bins', '3 bins'),
lty=c(1, 2), lwd=2, col=rgb(1, 0, 0, .8), cex=.8)
title('Piecewise Regressions', font.main=1)
This same separation as above allows us to analytically optimize for each bin separately. I.e. we run separate regressions on the split samples. From simple linear regression, we know the solutions are \[\begin{aligned}
\hat{b}_{0}(x,h) &= \hat{m}_{Y}(x,h)-\hat{b}_{1}(x,h)\hat{m}_{X}(x,h) \\
\hat{b}_{1}(x,h) &= \frac{\sum_{i}^{n(x,h)}(\hat{x}_{i}-\hat{m}_{X}(x,h))(\hat{y}_{i}-\hat{m}_{Y}(x,h))}{\sum_{i}^{n(x,h)}(\hat{x}_{i}-\hat{m}_{X}(x,h))^2} = \frac{\hat{c}_{XY}(x,h)}{\hat{v}_{X}(x,h)}.
\end{aligned}\]
Code
# See that two methods give the same predictions
# Piecewise: Two Schooling Bins
xy_hat[, 'xcc'] <- cut(xy_hat[, 'x_hat'], breaks=c(2.5, 9.5, 16.5))
preg_c <- lm(y_hat ~ xcc*x_hat, data=xy_hat)
y_fit_piecewise <- predict(preg_c)
## Split Sample Regressions
xy_split_hat <- split( xy_hat, xy_hat[, 'xcc'])
y_fit_split <- vector('list', length(xy_split_hat))
names(y_fit_split) <- names(xy_split_hat)
for (i in seq_along(xy_split_hat)) {
reg2 <- lm(y_hat ~ x_hat, xy_split_hat[[i]])
y_fit_split[[i]] <- predict(reg2)
}
# Any differences?
all( abs(y_fit_piecewise - unlist(y_fit_split)) < 1e-10 )
## [1] TRUE
Compare a simple regression to a regressogram and a piecewise regression. Examine the relationship between ‘Murder’ and ‘Urbanization’ in the USArrests dataset.
Code
xy2_hat <- USArrests[, c('UrbanPop', 'Murder')]
colnames(xy2_hat) <- c('x_hat', 'y_hat') # x_hat=urban population share, y_hat=murder arrests per 100K
# Globally Linear
reg <- lm(y_hat ~ x_hat, data=xy2_hat)
Now examine the local relationship between ‘Murder’ and ‘Urbanization’ using custom bins. Try to do so on your own before looking at the code below, which compares two piecewise regressions to a simple linear regression.
Code
# Diagnose Fit
#plot( fitted(reg), resid(reg), pch=16, col=grey(0, .5))
#plot( xy2_hat[, 'x_hat'], resid(reg), pch=16, col=grey(0, .5))
# Linear in 2 Pieces (subsets)
xcut2 <- cut(xy2_hat[, 'x_hat'], 2)
xy2_list2 <- split(xy2_hat, xcut2)
regs2 <- vector('list', length(xy2_list2))
names(regs2) <- names(xy2_list2)
for (i in seq_along(xy2_list2)) {
regs2[[i]] <- lm(y_hat ~ x_hat, data=xy2_list2[[i]])
}
# Linear in 3 Pieces (subsets or bins)
xcut3 <- cut(xy2_hat[, 'x_hat'], seq(32, 92, by=20), include.lowest=TRUE) # Finer Bins
xy2_list3 <- split(xy2_hat, xcut3)
regs3 <- vector('list', length(xy2_list3))
names(regs3) <- names(xy2_list3)
for (i in seq_along(xy2_list3)) {
regs3[[i]] <- lm(y_hat ~ x_hat, data=xy2_list3[[i]])
}
## Make Predictions
fit1 <- data.frame(y_fit=predict(reg),
x_hat=reg[['model']][, 'x_hat'])
fit1 <- fit1[order(fit1[, 'x_hat']), ]
fit2 <- vector('list', length(regs2))
for (i in seq_along(regs2)) {
fit2[[i]] <- data.frame(y_fit=predict(regs2[[i]]),
x_hat=regs2[[i]][['model']][, 'x_hat'])
}
fit2 <- do.call(rbind, fit2)
fit2 <- fit2[order(fit2[, 'x_hat']), ]
fit3 <- vector('list', length(regs3))
for (i in seq_along(regs3)) {
fit3[[i]] <- data.frame(y_fit=predict(regs3[[i]]),
x_hat=regs3[[i]][['model']][, 'x_hat'])
}
fit3 <- do.call(rbind, fit3)
fit3 <- fit3[order(fit3[, 'x_hat']), ]
# Compare Predictions
plot(y_hat ~ x_hat, pch=16, col=grey(0, .5), data=xy2_hat)
lines(y_fit ~ x_hat, fit1, lwd=2, col=rgb(1, 0, 0, .8))
lines(y_fit ~ x_hat, fit2, lwd=2, col=rgb(1, 0, 0, .8), lty=2)
lines(y_fit ~ x_hat, fit3, lwd=2, col=rgb(1, 0, 0, .8), lty=3)
legend('topleft',
legend=c('Globally Linear',
'Piecewise Linear (2)',
'Piecewise Linear (3)'),
lty=1:3, lwd=2, col=rgb(1, 0, 0, .8), cex=.8)
For many things, a simple linear regression, regressograms, or piecewise regression is “good enough”. Simple linear regressions struggle with nonlinear relationships but are very easy to run with a computer. Regressograms and piecewise regressions are intuitive ways to capture nonlinear relationships that are computationally efficient but have obvious problems where the bins change. Sometimes we want smoother predictions or to estimate derivatives (gradients).
Sliding Windows
To avoid the jumps at bin edges, we let the window slide along with \(x\). We first write a bin as a set of regression weights, using the uniform kernel from kernel density estimation, and then move those weights from one design point to the next.
Locally Constant Regression
For a bin of half-width \(h\) around the design point \(x\), the “uniform weights” are \[\begin{aligned}
k_{U}\left( \hat{x}_{i}, x, h \right)
&= \frac{\mathbf{1}\left(\frac{|\hat{x}_{i}-x|}{h} \leq 1\right)}{2}
= \frac{\mathbf{1}\left( \hat{x}_{i}~\text{in}~\left[ x-h, x + h\right] \right) }{2}
\end{aligned}\] which is nearly identical to \(\hat{d}_{i}(x,h)/ 2\), but uses intervals \([]\) instead of \((]\).
Here is an example going into the details of the weights.
Code
## Generate Sample Data
x0_hat <- 1:5
y0_hat <- round( rnorm(length(x0_hat)) , 2)
y0_hat
## [1] -0.63 0.18 -0.84 1.60 0.33
## plot(x0_hat, y0_hat)
## Manually Compute Estimate at x=3
h <- 1
k3_weights <- dunif( abs(x0_hat-3)/h, -1, 1) #(x0_hat >= 2)*(x0_hat <= 4)
k3_weights
## [1] 0.0 0.5 0.5 0.5 0.0
# Divide by the total so the weights sum to one
w3 <- k3_weights/sum(k3_weights)
y_pred_3 <- sum(w3*y0_hat)
y_pred_3
## [1] 0.3133333
# Equals the simple average of the y's within h of x=3
mean( y0_hat[ abs(x0_hat-3) <= h ] )
## [1] 0.3133333
Once a bin is written as weights around a design point \(x\), nothing requires \(x\) to be a bin midpoint, so we can slide the window and fit a local average at every design point.
A locally constant (kernel) regression fits a single constant at each design point \(x\). For data near \(x\), the model is \[\begin{aligned}
\hat{y}_{i} &= b_{0}(x,h) + \hat{e}_{i},
\end{aligned}\] and the coefficient minimizes a kernel-weighted sum of squared residuals: \[\min_{b_{0}(x, h)} \sum_{i=1}^{n}\left[\hat{y}_{i} - b_{0}(x, h) \right]^{2} k\left(\hat{x}_{i}, x, h\right),\] where \(k(\hat{x}_{i}, x, h)\) is a kernel weight that down-weights observations far from \(x\).
The locally constant regression is useful as the simplest smoother that avoids jumps at fixed bin edges: at each \(x\) it returns a kernel-weighted average of nearby \(\hat{y}_{i}\). With a uniform kernel and exclusive bins it reduces to a regressogram, so the regressogram is just its simplest special case.
For some intuition, note that this model is a weighted version of the simple constant regression. The weighted OLS estimator with uniform kernel weights \(k_{U}\) yields \[\begin{aligned}
& \min_{b_{0}(x,h)}~ \sum_{i}^{n}\left[\hat{e}_{i} \right]^2 k_{U}\left( \hat{x}_{i}, x, h \right) \\
\Rightarrow & -2 \sum_{i}^{n}\left[\hat{y}_{i}- b_{0}(x, h) \right] k_{U}\left(\hat{x}_{i}, x, h\right) = 0 \\
\Rightarrow & \hat{b}_{0}(x,h)
=\frac{\sum_{i} \hat{y}_{i} k_{U} \left( \hat{x}_{i}, x, h \right) }{ \sum_{i} k_{U}\left( \hat{x}_{i}, x, h \right) }
= \sum_{i} \hat{y}_{i} \left[ \frac{ k_{U} \left( \hat{x}_{i}, x, h \right) }{ \sum_{i} k_{U}\left( \hat{x}_{i}, x, h \right)} \right]
= \sum_{i} \hat{y}_{i} w_{i}(x, h),
\end{aligned}\] where weights are determined from \[\begin{aligned}
\sum_{i}^{n} k_{U} \left( \hat{x}_{i}, x, h \right) &= \sum_{i}^{n(x,h)} k_{U} \left( \hat{x}_{i}, x, h \right) + \sum_{i}^{n - n(x,h)} 0 = \frac{n(x, h)}{2} \\
w_{i}(x, h) &= \left[ \frac{ k_{U} \left( \hat{x}_{i}, x, h \right) }{ \sum_{i} k_{U}\left( \hat{x}_{i}, x, h \right)} \right]
= \frac{\mathbf{1}\left( |\hat{x}_{i} - x| \leq h \right)}{n(x, h)}
\end{aligned}\] So locally constant kernel regression recovers the weighted mean of \(\hat{y}_{i}\) around design point \(x\). If we use exclusive bins, then we are essentially running a regressogram. As such, the regressogram is more crude but can be estimated more easily.
Here is the locally constant regression on the running example. At each design point, we take the kernel-weighted average of \(\hat{y}_{i}\) for observations with nearby values of \(\hat{x}_{i}\).
Code
## MAIN EXAMPLE
# xy_hat defined above
x_values <- seq(3, 16, by=.5) # Design points
h <- 2 # bandwidth
# Kernel-weighted mean at each design point
y_pred_lc <- rep(NA, length(x_values))
for (i in seq_along(x_values)) {
ki <- dunif( abs(xy_hat[, 'x_hat']-x_values[i])/h, -1, 1) # the uniform kernel k_U
wi <- ki/sum(ki) # weights sum to one
y_pred_lc[i] <- sum(wi*xy_hat[, 'y_hat'])
}
# Plot
plot(y_hat ~ x_hat, data=xy_hat, pch=16, col=grey(0, .1),
main=NA, xlab='School', ylab='Wage')
lines(x_values, y_pred_lc, col=rgb(0, 0, 1, .8), lwd=2, type='o')
legend('topleft', title='Locally Constant',
legend='h=2',
lty=1, lwd=2, pch=1, col=rgb(0, 0, 1, .8), cex=.8)
Notes
When \(n\) is small, \(\hat{b}_{0}(x, h)\) is typically estimated at each unique observed value of \(\hat{x}_{i}\). For large datasets, you can select a subset or evenly spaced values of \(x\) for which to make predictions.
If \(\hat{x}_{i}\) represents time, then the locally constant regression is also called a moving average. We can create a weighted moving average by using a different distribution. I.e. using the Normal dnorm instead of the Uniform dunif.
The basic idea also generalizes to other kernels. As such, a kernel regression using uniform weights is often called a “naive kernel regression”. Typically, kernel regressions use kernels that weight nearby observations more heavily.
Dividing the kernel weights by their sum, as in the locally constant code above, makes the prediction an average of the \(\hat{y}_{i}\) in the window, even near the edges where part of the window is empty. Skipping that step gives a weighted sum rather than an average, as the example below shows.
Here is an example going into the details of the weights.
Code
## Generate Sample Data
x0_hat <- 1:5
y0_hat <- rnorm(length(x0_hat))
y0_hat
## [1] -0.8204684 0.4874291 0.7383247 0.5757814 -0.3053884
## Manually Compute Estimate at x=5
h <- 2
k5_weights <- dunif( abs(x0_hat-5)/h, -1, 1) #(x0_hat >= 3)*(x0_hat =< 7)/2
k5_weights
## [1] 0.0 0.0 0.5 0.5 0.5
y_pred_5 <- sum(k5_weights*y0_hat)
y_pred_5
## [1] 0.5043588
## Divide by the total so the weights sum to one
k5_weights_ec <- k5_weights/sum(k5_weights)
y_pred_5_ec <- sum(k5_weights_ec*y0_hat)
y_pred_5_ec
## [1] 0.3362392
Weighted Regression
Interestingly, we can also obtain the piecewise regression estimates from weighted least squares regression. For some specific design point, \(x\), we can find \(\hat{b}_{0}(x,h)\) and \(\hat{b}_{1}(x,h)\) by minimizing \[\begin{aligned}
\sum_{i}^{n}\left[ \hat{e}_{i} \right]^2 \hat{d}_{i}(x,h) &= \sum_{i}^{n}\left[ \hat{y}_{i}- b_{0}(x,h) - b_{1}(x,h) \hat{x}_{i} \right]^2 \hat{d}_{i}(x,h) \\
&= \sum_{i}^{n(x,h)}\left[\hat{y}_{i}- b_{0}\left(x,h\right) - b_{1}(x,h) \hat{x}_{i} \right]^2
\end{aligned}\]
We get nearly identical results if we instead use the uniform weights \(k_{U}\left( \hat{x}_{i}, x, h \right)\) from locally constant regression. As such we can see that \[\begin{aligned}
\sum_{i}^{n}\left[ \hat{e}_{i} \right]^2 k_{U}\left( \hat{x}_{i}, x, h \right)
&\approx \sum_{i}^{n}\left[ \hat{e}_{i} \right]^2 \hat{d}_{i}(x,h) / 2 \\
&= \sum_{i}^{n(x,h)}\left[\hat{y}_{i}- b_{0}\left(x,h\right) - b_{1}(x,h) \hat{x}_{i} \right]^2 / 2,
\end{aligned}\] The constant term \(1/2\) is irrelevant to finding the optimal solution (you can check the math yourself).
Here is the weighted regression on the running example, at the midpoint of the first bin of the two-bin piecewise regression. The code prints the weighted prediction and then the piecewise prediction.
Code
# xy_hat and preg_c defined above
## Weighted Regression
x <- 6 # midpoint of the first schooling bin
h <- 3.5 # window [2.5, 9.5] covers schooling levels 3 to 9
k <- abs(xy_hat[, 'x_hat']-x) /h
k_weights <- dunif(k, -1, 1) # the uniform kernel k_U
preg_k <- lm(y_hat ~ x_hat, data=xy_hat, weights=k_weights)
predict(preg_k, newdata=data.frame(x_hat=x))
## 1
## 4.378248
# Compare to Piecewise
predict(preg_c, newdata=data.frame(x_hat=x, xcc='(2.5,9.5]'))
## 1
## 4.378248
# Both predictions are the same, about 4.378.
Check the same result at a different design point. Use the second bin, (9.5,16.5].
Code
x <- 13 # midpoint of the second schooling bin
h <- 3.5 # window [9.5, 16.5] covers schooling levels 10 to 16
## Weighted Regression
k_weights <- dunif( abs(xy_hat[, 'x_hat']-x)/h, -1, 1)
preg_k <- lm(y_hat ~ x_hat, data=xy_hat, weights=k_weights)
predict(preg_k, newdata=data.frame(x_hat=x))
## 1
## 6.592354
# Compare to Piecewise
predict(preg_c, newdata=data.frame(x_hat=x, xcc='(9.5,16.5]'))
## 1
## 6.592354
Local Linear Regression
We can also add a slope term to the sliding window, so that each local fit can follow the trend within its window.
A local linear regression conducts a linear regression for each design point \(x\) using a subsample of data around it (Fan and Gijbels 1996). For data near \(x\), the model is \[\begin{aligned}
\hat{y}_{i} &= b_{0}(x,h) + b_{1}(x,h) \hat{x}_{i} + \hat{e}_{i},
\end{aligned}\] and the coefficients minimize a kernel-weighted sum of squared residuals: \[\begin{aligned}
\min_{b_{0}(x, h),b_{1}(x, h)}~ \sum_{i}^{n}\left[\hat{y}_{i}- b_{0}(x, h) - b_{1}(x,h) \hat{x}_{i} \right]^2 k\left( \hat{x}_{i}, x, h\right),
\end{aligned}\] where \(k(\hat{x}_{i}, x, h)\) is a kernel weight. The prediction at \(x\) is \(\hat{y}^{\mathrm{pred}}(x) = \hat{b}_{0}(x,h) + \hat{b}_{1}(x,h) x\).
Local linear regression is useful when the relationship slopes near the edges of the data, where the locally constant fit is biased, and it also reports a local slope \(\hat{b}_{1}(x,h)\) at every design point. It is the sliding-window version of piecewise regression: with the uniform kernel \(k_{U}\) and a design point at a bin midpoint, it reproduces the piecewise fit, as the weighted regression above showed. Deriving the optimal values \(\hat{b}_{0}(x, h)\) and \(\hat{b}_{1}(x,h)\) for \(k_{U}\) is left as a homework exercise.
Code
## MAIN EXAMPLE
# 'Naive' Smoother
pred_fun <- function(x0, h){
# Assign equal weight to observations within h distance to x0
# 0 weight for all other observations
ki <- abs(xy_hat[, 'x_hat']-x0)/h
ki <- dunif(ki, -1, 1) ## The uniform kernel k_U; could change, e.g. dnorm(ki)
wi <- ki/sum(ki, na.rm=TRUE) # always sum to 1 (for edge-effects)
# run regression with weighted data
llls_i <- lm(y_hat ~ x_hat, data=xy_hat, weights=wi)
y_pred <- predict(llls_i, newdata=data.frame(x_hat=x0))
}
x_values <- seq(4, 16, by=1) # Design points
# Fine Bins
y_pred_1 <- rep(NA, length(x_values))
for(i in seq_along(x_values)){
y_pred_1[i] <- pred_fun(x0=x_values[i], h=2)
}
# Coarse Bins
y_pred_2 <- rep(NA, length(x_values))
for(i in seq_along(x_values)){
y_pred_2[i] <- pred_fun(x0=x_values[i], h=6)
}
# Plot
plot(y_hat ~ x_hat, pch=16, data=xy_hat, col=grey(0, .1),
main=NA, ylab='Wage', xlab='School')
lines(x_values, y_pred_1, col=rgb(1, 0, 0, .8), lwd=2, type='o')
lines(x_values, y_pred_2, col=rgb(1, 0, 0, .8), lwd=2, lty=2, type='o')
legend('topleft', title='Local Linear',
legend=c('h=2', 'h=6'),
lty=c(1, 2), lwd=2, pch=1, col=rgb(1, 0, 0, .8), cex=.8,
seg.len=4) # long enough that the circle does not hide the dashes
Compare local linear and locally constant regressions using https://shinyserv.es/shiny/kreg/, with degree \(0\) and \(1\). Also try different kernels and different datasets.
Examine the local relationship between ‘Murder’ and ‘Urbanization’ in the USArrests dataset using local linear regression.
Code
# xy2_hat defined above
# Same smoother as pred_fun, but for any data with columns x_hat and y_hat
pred_fun_xy <- function(x0, h, xy){
ki <- abs(xy[, 'x_hat']-x0)/h
ki <- dunif(ki, -1, 1) ## The uniform kernel k_U
wi <- ki/sum(ki, na.rm=TRUE)
llls_i <- lm(y_hat ~ x_hat, data=xy, weights=wi)
y_pred <- predict(llls_i, newdata=data.frame(x_hat=x0))
}
x_values <- sort(unique(xy2_hat[, 'x_hat']))
plot(y_hat ~ x_hat, pch=16, data=xy2_hat, col=grey(0, .5),
ylab='Murder Rate', xlab='Urban Population (%)')
y_pred_1 <- numeric(length(x_values))
for (i in seq_along(x_values)) {
y_pred_1[i] <- pred_fun_xy(x_values[i], h=2, xy=xy2_hat)
}
y_pred_2 <- numeric(length(x_values))
for (i in seq_along(x_values)) {
y_pred_2[i] <- pred_fun_xy(x_values[i], h=20, xy=xy2_hat)
}
lines(x_values, y_pred_1, col=rgb(1, 0, 0, .8), lwd=2, type='o')
lines(x_values, y_pred_2, col=rgb(1, 0, 0, .8), lwd=2, lty=2, type='o')
legend('topleft', title='Local Linear',
legend=c('h=2 ', 'h=20'),
lty=c(1, 2), lwd=2, pch=1, col=rgb(1, 0, 0, .8), cex=.8,
seg.len=4) # long enough that the circle does not hide the dashes
Exercises
Comment the script you wrote for this chapter, then restart R and check that it runs from a clean session, then check the script with AI as explained in Working with AI. Write three sentences from memory on the main statistical idea of this chapter, and ask the assistant what is wrong, vague, or missing. Finish with your own questions about whatever you found hardest.
Explain the difference between a regressogram and a local linear regression. As you increase the number of bins in a regressogram, what happens to the number of observations in each bin, and how does that change the fitted values?
Using the Wages1 dataset from the Ecdat package, fit a regressogram with 4 bins and a piecewise regression with 4 bins for the relationship between wage and school. Compute the predicted value at school = 10 for each model.
Using the USArrests dataset, write R code to estimate a local linear regression of Murder on UrbanPop using uniform kernel weights with bandwidth \(h = 10\). Compute predictions at design points x_values <- seq(40, 80, by = 10) and plot them over the scatterplot.
The chapter derives the optimal locally constant coefficient \(\hat{b}_{0}(x,h)\) but leaves the local linear case as an exercise. Starting from the weighted objective \(\sum_{i} [\hat{y}_{i} - b_{0} - b_{1}\hat{x}_{i}]^2 \, k_{U}(\hat{x}_{i}, x, h)\), set the partial derivatives with respect to \(b_{0}\) and \(b_{1}\) to zero. Show that the solution is the simple-regression formula applied to the kernel-weighted data: \(\hat{b}_{1}(x,h) = \hat{c}_{XY}(x,h)/\hat{v}_{X}(x,h)\) and \(\hat{b}_{0}(x,h) = \hat{m}_{Y}(x,h) - \hat{b}_{1}(x,h)\hat{m}_{X}(x,h)\).
Recall
We let the relationship between wages and schooling in Wages1 bend by fitting local models. Regressograms fit a mean in each bin, and piecewise regressions fit a line in each bin. Locally constant and local linear regressions do the same within a sliding window, which avoids the jumps at the bin edges. The bandwidth \(h\) is the half-width that sets how local each estimate is: larger \(h\) pools more observations, and smaller \(h\) tracks the data more closely. The three-observation worked example \(\{(3, 2), (4, 4), (12, 3)\}\) split into two bins at \(h=7/2\) gave regressogram coefficients \(\hat{b}_{0}(x=6, h=7/2)=(2+4)/2=3\) and \(\hat{b}_{0}(x=13, h=7/2)=3\).
Fan, Jianqing, and Irène Gijbels. 1996.
Local Polynomial Modelling and Its Applications. Chapman; Hall.
https://doi.org/10.1201/9780203748725.
