25  Local Relationships

---

In Local Regression, we estimated nonparametric models with a single explanatory variable. We now extend those tools to multiple explanatory variables. The ideas carry over directly: instead of binning one variable, we bin each variable, and instead of fitting a line in pieces, we fit a plane in pieces.

This chapter asks three questions.

• How does the mean of \(\hat{y}\) vary with several explanatory variables when no single plane fits?
  • Multivariate regressogram: this gives each combination of bins its own mean, generalizing the univariate regressogram.
  • Piecewise regression: this bins one explanatory variable and fits a plane within each bin, where the bivariate version fits a line.
• How many bins should we use?
  • \(k\)-fold cross-validation: this model assessment check picks the number of bins by prediction error on held-out observations. It drops a group at a time, where leave-one-out cross-validation drops one observation.
• How does \(\hat{y}\) change with \(\hat{x}_{k}\) at different points in the data?
  • Finite differences and within-bin slopes: these estimate gradients the same two ways as before, now from the regressogram and the piecewise regression.

25.1 Regressograms

Setup

We want a flexible estimate of \(\mathbb{E}[Y_{i} \mid X_{i1}, X_{i2}]\) that does not commit to any functional form for either explanatory variable.

ImportantKey Definition

A multivariate regressogram cuts each explanatory variable into exclusive bins (each encoded as a dummy variable), then computes the prediction \(\hat{y}^{\mathrm{pred}}(x_{1},\ldots,x_{K})\) within each bin combination using OLS on the full set of dummies. It is the multivariate generalization of the univariate regressogram from Local Regression.

The multivariate regressogram is useful for estimating a conditional mean without committing to linearity or any particular interaction structure, because every cell gets its own prediction. This flexibility is also its main drawback. As the number of explanatory variables grows, the cells empty out and predictions in sparse regions become noisy. The simulated data below uses two predictors so we can plot it.

## MAIN EXAMPLE
## Simulate data: y_sim is a nonlinear function of x1_sim and x2_sim plus noise.
N <- 10000
e <- rnorm(N)
x1_sim <- seq(.1, 20, length.out=N)
x2_sim <- runif(N, 0, 1)
y_sim  <- 3*exp(-2*x2_sim + 1.5*x1_sim - .1*x1_sim^2)*x1_sim + e
xy_sim <- data.frame(x1_sim, x2_sim, y_sim)

## Create color palette (reused in later examples)
col_scale <- seq(min(y_sim)*1.1, max(y_sim)*1.1, length.out=401)
ycol_pal <- hcl.colors(length(col_scale), alpha=.5)
names(ycol_pal) <- sort(col_scale)

## Add legend (reused in later examples)
add_legend <- function(col_scale,
    yl=11,
    colfun=function(x){ hcl.colors(x, alpha=.5) },
    ...) {
  opar <- par(fig=c(0, 1, 0, 1), oma=c(0, 0, 0, 0), 
              mar=c(0, 0, 0, 0), new=TRUE)
  on.exit(par(opar))
  bin_mids <- hist(col_scale, plot=FALSE, breaks=yl-1)$mids
  plot(0, 0, type='n', bty='n', xaxt='n', yaxt='n')
  legend(...,
    legend=bin_mids,
    fill=colfun(length(bin_mids)),
    border=NA,
    bty='n')
}


## Plot Data
par(oma=c(0, 0, 0, 2))
plot(x1_sim ~ x2_sim, xy_sim,
    col=ycol_pal[cut(y_sim, col_scale)],
    pch=16, cex=.5,
    main=NA)
title('Raw Data', font.main=1)
add_legend(x='topright', col_scale=col_scale,
    yl=6, inset=c(0, .05), title='y')

OLS Baseline

We start with a global linear baseline, fitted without an interaction.

## MAIN EXAMPLE
## OLS (without interaction; the additive baseline)
reg <- lm(y_sim ~ x1_sim+x2_sim, data=xy_sim)

## Grid Points for Prediction
# X1 bins
l1 <- 11
bks1 <- seq(0, 20, length.out=l1)
h1 <- diff(bks1)[1]/2
mids1 <- bks1[-1]-h1
# X2 bins
l2 <- 11
bks2 <- seq(0, 1, length.out=l2)
h2 <- diff(bks2)[1]/2
mids2 <- bks2[-1]-h2
# Grid
pred_x <- expand.grid(x1_sim=mids1, x2_sim=mids2)

## OLS Predictions
y_pred_ols <- predict(reg, newdata=pred_x)
pred_df_ols  <- cbind(y_pred_ols, pred_x)

## Plot Predictions
par(oma=c(0, 0, 0, 2))
plot(x1_sim ~ x2_sim, pred_df_ols,
    col=ycol_pal[cut(y_pred_ols, col_scale)],
    pch=15, cex=2,
    main=NA)
title('OLS Predictions', font.main=1)
add_legend(x='topright', col_scale=col_scale,
    yl=6, inset=c(0, .05), title='y')

Multivariate Bins

## MAIN EXAMPLE
##################
# Multivariate Regressogram
##################

## Cut each variable into bins, then interact them: each (x1c_sim, x2c_sim) cell
## gets its own mean. With the `*` interaction the model is fully saturated.
xy_sim$x1c_sim <- cut(xy_sim$x1_sim, bks1)
xy_sim$x2c_sim <- cut(xy_sim$x2_sim, bks2)

## Regressogram (interacted bins -> one mean per (x1c_sim, x2c_sim) cell)
reg <- lm(y_sim ~ x1c_sim*x2c_sim, data=xy_sim)

## Predicted Values
## For Points in Middle of Each Bin
pred_df_rgrm <- expand.grid(
    x1c_sim=levels(xy_sim$x1c_sim),
    x2c_sim=levels(xy_sim$x2c_sim))
pred_df_rgrm$y_pred <- predict(reg, newdata=pred_df_rgrm)
pred_df_rgrm <- cbind(pred_df_rgrm, pred_x)

## Plot Predictions
par(oma=c(0, 0, 0, 2))
plot(x1_sim ~ x2_sim, pred_df_rgrm,
    col=ycol_pal[cut(pred_df_rgrm$y_pred, col_scale)],
    pch=15, cex=2,
    main=NA)
title('Regressogram Predictions', font.main=1)
add_legend(x='topright', col_scale=col_scale,
    yl=6, inset=c(0, .05), title='y')

TipTest Yourself

The regressogram divides the \((x_{1}, x_{2})\) space into a grid of bins and estimates the mean of \(y\) within each bin. Compared to OLS, the regressogram makes no assumption about functional form: it can capture nonlinear and interaction effects. The trade-off is that it requires more data to fill each bin, especially as the number of variables grows.

NoteMust Know

A regressogram prediction is just the average of \(\hat{y}_{i}\) within a bin. Consider six observations, with \(x_{1}\) and \(x_{2}\) each split into “low” and “high”.

xy0_hat <- data.frame(
    y_hat  = c(2, 4, 5, 9, 8, 12),
    x1_hat = c('low', 'low', 'low', 'high', 'high', 'high'),
    x2_hat = c('low', 'low', 'high', 'low', 'high', 'high'))

# Regressogram prediction = mean of y_hat within each (x1_hat, x2_hat) cell
aggregate(y_hat ~ x1_hat + x2_hat, data=xy0_hat, FUN=mean)
##   x1_hat x2_hat y_hat
## 1   high   high    10
## 2    low   high     5
## 3   high    low     9
## 4    low    low     3

The fitted value for any point is its cell mean: \((2+4)/2 = 3\) for the (low, low) cell, \(5\) for (low, high), \(9\) for (high, low), and \((8+12)/2 = 10\) for (high, high).

25.2 Local Regressions

Piecewise Model

A regressogram is locally constant, so every point in a bin gets the same prediction. We sometimes want slopes within each bin instead of a flat shelf.

ImportantKey Definition

A multivariate piecewise regression bins one explanatory variable and fits a separate OLS plane on all the explanatory variables within each bin.

A piecewise regression is useful when you want something between the flat regressogram and the global OLS plane. It tracks gradients within each bin, so predictions vary smoothly within a bin and jump only at the bin edges. We bin one variable, then fit a full linear model of \(\hat{y}\) on the explanatory variables within each bin.

## MAIN EXAMPLE
##################
# Multivariate Piecewise Regression
##################

## Bin x1_sim, then fit a linear model of y_sim on (x1_sim, x2_sim) within each bin
xy_sim$x1c_sim <- cut(xy_sim$x1_sim, bks1)
preg <- lm(y_sim ~ x1c_sim/(x1_sim + x2_sim), data=xy_sim)

## Predicted Values
## For the same grid of points used above
pred_df_preg <- pred_x
pred_df_preg$x1c_sim <- cut(pred_df_preg$x1_sim, bks1)
pred_df_preg$y_pred <- predict(preg, newdata=pred_df_preg)

## Plot Predictions
par(oma=c(0, 0, 0, 2))
plot(x1_sim ~ x2_sim, pred_df_preg,
    col=ycol_pal[cut(pred_df_preg$y_pred, col_scale)],
    pch=15, cex=2,
    main=NA)
title('Piecewise Regression Predictions', font.main=1)
add_legend(x='topright', col_scale=col_scale,
    yl=6, inset=c(0, .05), title='y')

The formula y_sim ~ x1c_sim/(x1_sim + x2_sim) gives each x1_sim bin its own intercept and its own slopes on x1_sim and x2_sim. Running one regression per data subset with split() gives identical predictions, as shown for the bivariate case in Local Regression. A fully local linear regression generalizes this further, replacing hard bins with kernel weights and a separate bandwidth for each variable; see Local Regression for the bivariate version.

TipTest Yourself

Compare the three sets of predictions in this chapter (OLS, regressogram, and piecewise). Which captures the curvature in x1_sim? Which is smoothest? Refit the piecewise model with finer bins (bks1 <- seq(0, 20, length.out=21)) and describe how the prediction plot changes.

Break Points

Kinks and discontinuities in the relationship between \(\hat{y}\) and \(\hat{x}_{k}\) can be modeled with factor variables, just like any other transformation. We can then use \(F\)-tests to examine whether a break is statistically significant. Both the interaction-model \(F\)-test (anova) and the Chow test (Chow 1960) below test the same null of “no break”, and give the same answer.

data(Wages1, package='Ecdat')
xy_hat <- Wages1[, c('school', 'wage')]
colnames(xy_hat) <- c('x_hat', 'y_hat') # x_hat=years of schooling, y_hat=hourly wage
reg <- lm(y_hat ~ x_hat, data=xy_hat)

# F-test for a break at school = 12
reg2 <- lm(y_hat ~ x_hat*I(x_hat>12), data=xy_hat)
anova(reg, reg2)
## Analysis of Variance Table
## 
## Model 1: y_hat ~ x_hat
## Model 2: y_hat ~ x_hat * I(x_hat > 12)
##   Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
## 1   3292 32386                                  
## 2   3290 32146  2     239.8 12.271 4.902e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Chow test for the same break
data_splits <- split(xy_hat, xy_hat$x_hat <= 12)
resids <- sapply(data_splits, function(xy_sub){
    reg_sub <- lm(y_hat ~ x_hat, data=xy_sub)
    sum(resid(reg_sub)^2)
})
Ns <- sapply(data_splits, function(xy_sub){ nrow(xy_sub) })

# Chow F-statistic: does splitting the sample reduce RSS enough?
Rt <- (sum(resid(reg)^2) - sum(resids))/sum(resids)
Rb <- (sum(Ns) - 2*reg$rank)/reg$rank
Ft <- Rt*Rb
pf(Ft, reg$rank, sum(Ns)-2*reg$rank, lower.tail=FALSE)
## [1] 4.901941e-06
# The p-value is about 5e-06, which is evidence of a break at 12 years.
# Twelve years is the end of secondary school.

# See also
# strucchange::sctest(y_hat ~ x_hat, data=xy_hat, type='Chow', point=.5)
# segmented::segmented(reg)

A small \(p\)-value is evidence that the line changes at the candidate break. The ratio Rt in the Chow test is the drop in the residual sum of squares, divided by the unrestricted residual sum of squares. This is the same loss-in-fit ratio as in Model Fit.

NoteMust Know

The Chow test checks whether a structural break exists at a specified point. What happens if you do not know where the break is? How does searching over many potential break points affect the validity of your test? (Hint: think about multiple testing.)

25.3 Model Selection

Leave-One-Out and K-Fold CV

A flexible model can fit the training data perfectly while predicting new observations badly, so we need a way to measure performance on data the model has not seen.

ImportantKey Definition

Cross-validation (CV) picks a model by minimizing out-of-sample prediction error. Leave-one-out CV drops one observation at a time, refits, and averages the squared prediction errors: \[ \min_{\mathbf{H}} \quad \frac{1}{n} \sum_{i=1}^{n} \left[ \hat{y}_{i} - \hat{y}^{\mathrm{pred}(-i)}(\hat{\mathbf{x}}_{i},\mathbf{H}) \right]^2, \] where \(\hat{\mathbf{x}}_{i}=[\hat{x}_{i1}, \ldots, \hat{x}_{iK}]\) collects the explanatory variables of observation \(i\), \(\hat{y}^{\mathrm{pred}(-i)}(\hat{\mathbf{x}}_{i},\mathbf{H})\) is the model’s prediction at \(\hat{\mathbf{x}}_{i}\) from a fit that excluded observation \(i\), and \(\mathbf{H}\) is the bandwidth matrix. \(k\)-fold CV drops one of \(k\) random groups at a time instead of one observation.

Cross-validation is useful because it directly penalizes overfitting: a too-flexible model predicts held-out observations poorly even though it fits in-sample data tightly. With a weighted least squares regression on three explanatory variables, the bandwidth matrix has the form \[ \mathbf{H}=\begin{pmatrix} h_{1} & 0 & 0 \\ 0 & h_{2} & 0 \\ 0 & 0 & h_{3} \\ \end{pmatrix}, \] where each \(h_{k}\) is the bandwidth for variable \(\hat{x}_{k}\). In the discrete-bin setting below, the role of the bandwidth is played by the number of bins \(L\): more bins means finer resolution (a smaller effective bandwidth) and fewer observations per bin.

There are many types of cross-validation (Arlot and Celisse 2010; Bates et al. 2023). Generalized cross-validation adjusts for the degrees of freedom, whereas the npreg function in R uses least-squares cross-validation (Racine 2019, 74) by default. You can refer to extensions on a case by case basis. In applied work, Adamson (2025) uses leave-one-out cross-validation (Allen 1974; Stone 1974; Hall and Racine 2015) to choose the bandwidth of geographically weighted regressions on ancient Greek city-states. These are local regressions whose kernel weights depend on the distance between locations (Brunsdon et al. 1996; Fotheringham et al. 2003).

TipTest Yourself

Cross-validation works by repeatedly holding out one observation, fitting the model on the rest, and measuring how well the model predicts the held-out point. A model that overfits (e.g., too many bins) will have low in-sample error but high cross-validation error. The bandwidth \(\mathbf{H}\) that minimizes cross-validation error balances flexibility against overfitting.

To make this concrete, we use cross-validation to choose the number of x1_sim bins for the piecewise regression above. We split the data into 5 folds, fit the model on 4 of them, and measure the squared prediction error on the held-out fold. With almost no noise, finer bins tend to keep lowering the error, so we first add more noise to y_sim (standard deviation \(2500\)).

## MAIN EXAMPLE
## 5-fold cross-validation for the number of x1_sim bins
n_folds <- 5
fold_id <- sample(rep(1:n_folds, length.out=nrow(xy_sim)))
bin_counts <- c(2, 6, 10, 14, 18, 22)

## Noisier copy of the data, so that too many bins overfit
xy_cv_sim <- xy_sim
xy_cv_sim$y_sim <- xy_sim$y_sim + rnorm(nrow(xy_sim), sd=2500)
# The noise in y_sim alone is tiny compared to its curvature.
# Without more, the error would keep falling across all these bin counts.

cv_mse <- sapply(bin_counts, function(L){
    bks <- seq(0, 20, length.out=L+1)
    fold_mse <- sapply(1:n_folds, function(k){
        xy_train_sim <- xy_cv_sim[fold_id != k, ]
        xy_test_sim  <- xy_cv_sim[fold_id == k, ]
        xy_train_sim$x1c_sim <- cut(xy_train_sim$x1_sim, bks)
        xy_test_sim$x1c_sim  <- cut(xy_test_sim$x1_sim, bks)
        reg_k <- lm(y_sim ~ x1c_sim/(x1_sim + x2_sim), data=xy_train_sim)
        y_pred_k <- predict(reg_k, newdata=xy_test_sim)
        mean((xy_test_sim$y_sim - y_pred_k)^2, na.rm=TRUE)
    })
    mean(fold_mse)
})

## Plot CV error against model complexity
plot(bin_counts, cv_mse, type='o', pch=16,
    xlab='Number of x1 bins', ylab='5-fold CV mean squared error',
    main=NA)
L_star <- bin_counts[which.min(cv_mse)]
abline(v=L_star, col=rgb(1, 0, 0, .8), lwd=2)
title(paste0('CV-optimal: ', L_star, ' bins'), font.main=1)

# The error is smallest at 10 bins, between too few and too many.

Too few bins cannot capture the curvature in x1_sim (high bias); too many bins split the data so finely that each local regression is noisy (high variance). The cross-validated error is usually smallest in between.

TipTest Yourself

The cross-validation above selects bins for the piecewise model. Adapt the code to instead select bins for the regressogram lm(y_sim ~ x1c_sim*x2c_sim), cutting both x1_sim and x2_sim into L bins. What happens to the cross-validated error once L is large enough that some grid cells are empty? (Hint: empty cells give NA predictions.)

25.4 Hypothesis Testing

Finite Differences

A gradient describes how \(\hat{y}^{\mathrm{pred}}\) changes with each \(x_{k}\) at a point. We first summarize gradients, then test them.

ImportantKey Definition

A gradient of a fitted multivariate model is the partial derivative \(\partial \hat{y}^{\mathrm{pred}}/\partial x_{k}\) evaluated at a point. For a regressogram with no slope coefficients, it must be approximated by a finite difference between adjacent bins, \[\begin{aligned} \hat{b}_{k}(\mathbf{x}) &= \frac{ \hat{y}^{\mathrm{pred}}(x_{1},...,x_{k}+ \frac{d}{2}...,x_{K})-\hat{y}^{\mathrm{pred}}(x_{1},...,x_{k}-\frac{d}{2}...,x_{K})}{d}, \end{aligned}\] where the step \(d\) is the bin width and \(x_{k}\) is the edge between two bins, while a piecewise or local linear regression has within-bin slopes that are themselves gradient estimates.

The gradient is useful because, unlike the single OLS coefficient, it varies with \(\mathbf{x}\). This lets us see whether the marginal effect of \(\hat{x}_{k}\) is the same at small, medium, and large values of the other regressors. We can compute it two ways:

1. For regressograms, approximate the gradient with a small finite difference between adjacent bins (formula above).

2. For split-sample regressions or local linear regressions, read the within-bin slope coefficients \(\hat{b}_{k}(\mathbf{x})\) directly as gradient estimates.

After computing gradients, you can summarize them in various plots: a histogram of \(\hat{b}_{k}(\mathbf{x})\), a scatterplot of \(\hat{b}_{k}(\mathbf{x})\) against \(x_{k}\), or the sorted gradients shown with their confidence intervals (Chaudhuri and Marron 1999; Henderson et al. 2012).

You may also be interested in a particular gradient or a single summary statistic. For example, a bivariate regressogram can estimate the marginal effect of \(\hat{x}_{1}\) at the means, \(\hat{b}_{1}(\mathbf{x}=[\hat{m}_{X_{1}}, \hat{m}_{X_{2}}])\), where \(\hat{m}_{X_{k}}\) is the mean of \(\hat{x}_{k}\). You may also be interested in the mean of the marginal effects (sometimes said simply as “average effect”), which averages the marginal effect over all datapoints in the dataset: \(1/n \sum_{i=1}^{n} \hat{b}_{1}(\mathbf{x}=\hat{\mathbf{x}}_{i})\), or the median marginal effect. Such statistics are single numbers that can be presented similar to an OLS regression table, where each row corresponds to a variable and each cell has two elements: “mean gradient (sd gradient)”. A single average can hide gradients that change sign (Henderson et al. 2012; Adamson 2021).

We illustrate both methods on the models fit earlier in this chapter. Method 1 differences the regressogram predictions; method 2 reads the slopes directly off the piecewise regression.

## MAIN EXAMPLE
##################
# Method 1: finite-difference gradients from the regressogram
##################

## Arrange the regressogram predictions into an x1_sim-by-x2_sim grid
grid_y_pred <- matrix(pred_df_rgrm$y_pred,
    nrow=length(mids1), ncol=length(mids2))

## Centered finite difference at each interior x1_sim bin edge, with d equal to the bin width
## The points x1 - d/2 and x1 + d/2 are the midpoints of the two adjacent bins
d1 <- diff(bks1)[1]
grad_x1 <- apply(grid_y_pred, 2, diff)/d1
x1_edges <- mids1[-1] - d1/2

## Summarize the gradient distribution
par(mfrow=c(1, 2))
hist(grad_x1, breaks=20, border=NA, freq=FALSE,
    main=NA, xlab=expression(hat(b)[1](x)))
title('Distribution of x1-gradients', font.main=1)
plot(rep(x1_edges, ncol(grad_x1)), grad_x1,
    pch=16, col=grey(0, .5),
    xlab=expression(x[1]), ylab=expression(hat(b)[1](x)),
    main=NA)
title('Gradient vs x1', font.main=1)

Slope Coefficients

## MAIN EXAMPLE
##################
# Method 2: slope coefficients from the piecewise regression
##################

## The within-bin slope on x1_sim is itself a gradient estimate
coef_preg <- coef(preg)
slope_x1 <- coef_preg[grep(':x1_sim$', names(coef_preg))]

## Mean and standard deviation of the x1_sim-gradient
c(mean=mean(slope_x1, na.rm=TRUE), sd=sd(slope_x1, na.rm=TRUE))
##        mean          sd 
##   0.8222098 422.1854069

The mean gradient is a single-number summary, much like an OLS slope. The standard deviation across bins measures how much the gradient changes over the data, which a single OLS coefficient cannot show. To test whether a mean gradient differs from zero, bootstrap it exactly as in Gradients. To test whether all gradients equal zero, use Racine’s statistic instead.

NoteMust Know

Suppose the regressogram predicts \(\hat{y}^{\mathrm{pred}}(x_{1}-1) = 8\) in one \(x_{1}\)-bin and \(\hat{y}^{\mathrm{pred}}(x_{1}+1) = 14\) in the next, so the bin midpoints are \(2\) units apart. The estimated \(x_{1}\)-gradient at the edge \(x_{1}\) between the two bins is \((14 - 8)/2 = 3\): moving one unit in \(x_{1}\) raises the prediction by about \(3\).

25.5 Exercises

1. 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.

2. A regressogram bins each explanatory variable and estimates a separate mean within each bin combination. How does this differ from a standard OLS regression with the same variables? What advantage does it have when the true relationship is nonlinear, and what is its main drawback as the number of explanatory variables grows?

3. Using the simulated data from this chapter (x1_sim <- seq(.1, 20, length.out = 1000), x2_sim <- runif(1000, 0, 1), y_sim <- 3*exp(-2*x2_sim + 1.5*x1_sim - .1*x1_sim^2)*x1_sim + rnorm(1000)), build a regressogram by cutting x1_sim into 5 bins and x2_sim into 5 bins. Compute the mean of y_sim within each of the 25 bin combinations. Which bin combination has the highest predicted value?

4. Using lm(), fit two models to the mtcars dataset with mpg as the outcome: (a) a standard linear model mpg ~ wt + hp, and (b) a regressogram model where wt and hp are each cut into 4 bins using cut() and interacted (mpg ~ cut(wt,4) * cut(hp,4)). Compare the \(\hat{R}^2\) of each model. Plot the predicted values from both models against the observed mpg.

Further Reading

• https://cran.r-project.org/web/packages/np/vignettes/np.pdf – the np package vignette for kernel and local-linear nonparametric regression in R.
• https://www.econometrics-with-r.org/11-rwmad.html – nonparametric methods and their connection to regression, including bin-smoothers.

Recall

This chapter generalized regressograms, piecewise regression, and cross-validation to the multivariate setting, and showed two ways to compute gradients: finite differences from the regressogram, and within-bin slopes from the piecewise model. The 5-fold CV exercise on y_sim ~ x1c_sim/(x1_sim + x2_sim) illustrated the bias-variance trade-off. We tried \(\{2, 6, 10, 14, 18, 22\}\) bins for x1_sim. After we added more noise to y_sim, the cross-validated MSE fell to its minimum at the CV-optimal \(L^{*}\) and rose again as we kept adding bins.

Adamson, Jordan. 2021. “The Scope of Political Jurisdictions and Violence: Theory and Evidence from Africa.” Public Choice 186 (3): 467–90. https://doi.org/10.1007/s11127-019-00763-8.

Adamson, Jordan. 2025. “Trade and the Rise of Ancient Greek City-States.” Journal of Economic Behavior & Organization 235: 107035. https://doi.org/10.1016/j.jebo.2025.107035.

Allen, David M. 1974. “The Relationship Between Variable Selection and Data Augmentation and a Method for Prediction.” Technometrics 16 (1): 125–27. https://doi.org/10.1080/00401706.1974.10489157.

Arlot, Sylvain, and Alain Celisse. 2010. “A survey of cross-validation procedures for model selection.” Statistics Surveys 4 (none): 40–79. https://doi.org/10.1214/09-SS054.

Bates, Stephen, Trevor Hastie, and Robert Tibshirani. 2023. “Cross-Validation: What Does It Estimate and How Well Does It Do It?” Journal of the American Statistical Association 0 (0): 1–12. https://doi.org/10.1080/01621459.2023.2197686.

Brunsdon, Chris, A. Stewart Fotheringham, and Martin E. Charlton. 1996. “Geographically Weighted Regression: A Method for Exploring Spatial Nonstationarity.” Geographical Analysis 28 (4): 281–98. https://doi.org/10.1111/j.1538-4632.1996.tb00936.x.

Chaudhuri, Probal, and J. S. Marron. 1999. “SiZer for Exploration of Structures in Curves.” Journal of the American Statistical Association 94 (447): 807–23. https://doi.org/10.1080/01621459.1999.10474186.

Chow, Gregory C. 1960. “Tests of Equality Between Sets of Coefficients in Two Linear Regressions.” Econometrica 28 (3): 591–605. https://doi.org/10.2307/1910133.

Fotheringham, A. S., C. Brunsdon, and M. Charlton. 2003. Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. Wiley. https://books.google.de/books?id=9DZgV1vXOuMC.

Hall, Peter G., and Jeffrey S. Racine. 2015. “Infinite Order Cross-Validated Local Polynomial Regression.” Journal of Econometrics 185 (2): 510–25. https://doi.org/10.1016/j.jeconom.2014.06.003.

Henderson, Daniel J., Subal C. Kumbhakar, and Christopher F. Parmeter. 2012. “A Simple Method to Visualize Results in Nonlinear Regression Models.” Economics Letters 117 (3): 578–81. https://doi.org/10.1016/j.econlet.2012.07.040.

Racine, Jeffrey S. 2019. An Introduction to the Advanced Theory and Practice of Nonparametric Econometrics: A Replicable Approach Using r. Cambridge University Press. https://doi.org/10.1017/9781108649841.

Stone, M. 1974. “Cross-Validatory Choice and Assessment of Statistical Predictions.” Journal of the Royal Statistical Society: Series B (Methodological) 36 (2): 111–33. https://doi.org/10.1111/j.2517-6161.1974.tb00994.x.