Local regression

Model definition

Local regression uses a data set consisting of observations one or more ‘independent’ or ‘predictor’ variables, and a ‘dependent’ or ‘response’ variable. The dataset will consist of a number ${\displaystyle n}$ observations. The observations of the predictor variable can be denoted ${\displaystyle x_1,\dots ,x_n}$, and corresponding observations of the response variable by ${\displaystyle Y_1,\dots ,Y_n}$.

For ease of presentation, the development below assumes a single predictor variable; the extension to multiple predictors (when the ${\displaystyle x_i}$ are vectors) is conceptually straightforward. A functional relationship between the predictor and response variables is assumed: $${\displaystyle Y_i=\mu (x_i)+{ϵ}_i}$$ where ${\displaystyle \mu (x)}$ is the unknown ‘smooth’ regression function to be estimated, and represents the conditional expectation of the response, given a value of the predictor variables. In theoretical work, the ‘smoothness’ of this function can be formally characterized by placing bounds on higher order derivatives. The ${\displaystyle {ϵ}_i}$ represents random error; for estimation purposes these are assumed to have mean zero. Stronger assumptions (e.g., independence and equal variance) may be made when assessing properties of the estimates.

Local regression then estimates the function ${\displaystyle \mu (x)}$, for one value of ${\displaystyle x}$ at a time. Since the function is assumed to be smooth, the most informative data points are those whose ${\displaystyle x_i}$ values are close to ${\displaystyle x}$. This is formalized with a bandwidth ${\displaystyle h}$ and a kernel or weight function ${\displaystyle W(\cdot )}$, with observations assigned weights $${\displaystyle w_i(x)=W\left(\frac{x_i-x}{h}\right).}$$ A typical choice of ${\displaystyle W}$, used by Cleveland in LOWESS, is ${\displaystyle W(u)=(1-|u{|}^3{)}^3}$ for ${\displaystyle |u|<1}$, although any similar function (peaked at ${\displaystyle u=0}$ and small or 0 for large values of ${\displaystyle u}$) can be used. Questions of bandwidth selection and specification (how large should ${\displaystyle h}$ be, and should it vary depending upon the fitting point ${\displaystyle x}$?) are deferred for now.

A local model (usually a low-order polynomial with degree ${\displaystyle p\le 3}$), expressed as $${\displaystyle \mu (x_i)\approx {\beta }_0+{\beta }_1(x_i-x)+\dots +{\beta }_p(x_i-x{)}^p}$$ is then fitted by weighted least squares: choose regression coefficients ${\displaystyle ({\hat{\beta }}_0,\dots ,{\hat{\beta }}_p)}$ to minimize $${\displaystyle \sum _{i=1}^n{w}_i(x){\left(Y_i-{\beta }_0-{\beta }_1(x_i-x)-\dots -{\beta }_p(x_i-x{)}^p\right)}^2.}$$ The local regression estimate of ${\displaystyle \mu (x)}$ is then simply the intercept estimate: $${\displaystyle \hat{\mu }(x)={\hat{\beta }}_0}$$ while the remaining coefficients can be interpreted (up to a factor of ${\displaystyle p!}$) as derivative estimates.

It is to be emphasized that the above procedure produces the estimate ${\displaystyle \hat{\mu }(x)}$ for one value of ${\displaystyle x}$. When considering a new value of ${\displaystyle x}$, a new set of weights ${\displaystyle w_i(x)}$ must be computed, and the regression coefficient estimated afresh.

Matrix representation of the local regression estimate

As with all least squares estimates, the estimated regression coefficients can be expressed in closed form (see Weighted least squares for details): $${\displaystyle \hat{\mathit{\beta }}={\left({\mathbf{X}}^{\mathsf{\text{T}}}\mathbf{W}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathsf{\text{T}}}\mathbf{W}\mathbf{y}}$$ where ${\displaystyle \hat{\mathit{\beta }}}$ is a vector of the local regression coefficients; ${\displaystyle \mathbf{X}}$ is the ${\displaystyle n\times (p+1)}$ design matrix with entries ${\displaystyle (x_i-x{)}^j}$; ${\displaystyle \mathbf{W}}$ is a diagonal matrix of the smoothing weights ${\displaystyle w_i(x)}$; and ${\displaystyle \mathbf{y}}$ is a vector of the responses ${\displaystyle Y_i}$.

This matrix representation is crucial for studying the theoretical properties of local regression estimates. With appropriate definitions of the design and weight matrices, it immediately generalizes to the multiple-predictor setting.

Selection issues: bandwidth, local model, fitting criteria

Implementation of local regression requires specification and selection of several components:

1. The bandwidth, and more generally the localized subsets of the data.
2. The degree of local polynomial, or more generally, the form of the local model.
3. The choice of weight function ${\displaystyle W(\cdot )}$.
4. The choice of fitting criterion (least squares or something else).

Each of these components has been the subject of extensive study; a summary is provided below.

Localized subsets of data; Bandwidth

The bandwidth ${\displaystyle h}$ controls the resolution of the local regression estimate. If h is too small, the estimate may show high-resolution features that represent noise in the data, rather than any real structure in the mean function. Conversely, if h is too large, the estimate will only show low-resolution features, and important structure may be lost. This is the bias-variance tradeoff; if h is too small, the estimate exhibits large variation; while at large h, the estimate exhibits large bias.

Careful choice of bandwidth is therefore crucial when applying local regression. Mathematical methods for bandwidth selection require, firstly, formal criteria to assess the performance of an estimate. One such criterion is prediction error: if a new observation is made at ${\displaystyle \tilde{x}}$, how well does the estimate ${\displaystyle \hat{\mu }(\tilde{x})}$ predict the new response ${\displaystyle \tilde{Y}}$?

Performance is often assessed using a squared-error loss function. The mean squared prediction error is $${\displaystyle \begin{array}{cc}\mathrm{E}{\left[\tilde{Y}-\hat{\mu }(\tilde{x})\right]}^2& =\mathrm{E}{\left[\tilde{Y}-\mu (x)+\mu (x)-\hat{\mu }(\tilde{x})\right]}^2\\ & =\mathrm{E}{\left[\tilde{Y}-\mu (x)\right]}^2+\mathrm{E}{\left[\mu (x)-\hat{\mu }(\tilde{x})\right]}^2.\end{array}}$$ The first term ${\displaystyle E{\left(\tilde{Y}-\mu (x)\right)}^2}$ is the random variation of the observation; this is entirely independent of the local regression estimate. The second term, ${\displaystyle \mathrm{E}{\left[\mu (x)-\hat{\mu }(\tilde{x})\right]}^2}$ is the mean squared estimation error. This relation shows that, for squared error loss, minimizing prediction error and estimation error are equivalent problems.

In global bandwidth selection, these measures can be integrated over the ${\displaystyle x}$ space ("mean integrated squared error", often used in theoretical work), or averaged over the actual ${\displaystyle x_i}$ (more useful for practical implementations). Some standard techniques from model selection can be readily adapted to local regression:

1. Cross Validation, which estimates the mean-squared prediction error.
2. Mallow's Cp and Akaike's Information Criterion, which estimate mean squared estimation error.
3. Other methods which attempt to estimate bias and variance variance components of the estimation error directly.

Any of these criteria can be minimized to produce an automatic bandwidth selector. Cleveland and Devlin prefer a graphical method (the M-plot) to visually display the bias-variance trade-off and guide bandwidth choice.

One question not addressed above is, how should the bandwidth depend upon the fitting point ${\displaystyle x}$? Often a constant bandwidth is used, while LOWESS and LOESS prefer a nearest-neighbor bandwidth, meaning h is smaller in regions with many data points. Formally, the smoothing parameter, ${\displaystyle \alpha }$, is the fraction of the total number n of data points that are used in each local fit. The subset of data used in each weighted least squares fit thus comprises the ${\displaystyle n\alpha }$ points (rounded to the next largest integer) whose explanatory variables' values are closest to the point at which the response is being estimated.

More sophisticated methods attempt to choose the bandwidth adaptively; that is, choose a bandwidth at each fitting point ${\displaystyle x}$ by applying criteria such as cross-validation locally within the smoothing window. An early example of this is Jerome H. Friedman's "supersmoother", which uses cross-validation to choose among local linear fits at different bandwidths.

Degree of local polynomials

Most sources, in both theoretical and computational work, use low-order polynomials as the local model, with polynomial degree ranging from 0 to 3.

The degree 0 (local constant) model is equivalent to a kernel smoother; usually credited to Èlizbar Nadaraya (1964) and G. S. Watson (1964). This is the simplest model to use, but can suffer from bias when fitting near boundaries of the dataset.

Local linear (degree 1) fitting can substantially reduce the boundary bias.

Local quadratic (degree 2) and local cubic (degree 3) can result in improved fits, particularly when the underlying mean function ${\displaystyle \mu (x)}$ has substantial curvature, or equivalently a large second derivative.

In theory, higher orders of polynomial can lead to faster convergence of the estimate ${\displaystyle \hat{\mu }(x)}$ to the true mean ${\displaystyle \mu (x)}$, provided that ${\displaystyle \mu (x)}$ has a sufficient number of derivatives. See C. J. Stone (1980). Generally, it takes a large sample size for this faster convergence to be realized. There are also computational and stability issues that arise, particularly for multivariate smoothing. It is generally not recommended to use local polynomials with degree greater than 3.

As with bandwidth selection, methods such as cross-validation can be used to compare the fits obtained with different degrees of polynomial.

Weight function

As mentioned above, the weight function gives the most weight to the data points nearest the point of estimation and the least weight to the data points that are furthest away. The use of the weights is based on the idea that points near each other in the explanatory variable space are more likely to be related to each other in a simple way than points that are further apart. Following this logic, points that are likely to follow the local model best influence the local model parameter estimates the most. Points that are less likely to actually conform to the local model have less influence on the local model parameter estimates.

Cleveland (1979) sets out four requirements for the weight function:

1. Non-negative: ${\displaystyle W(x)>0}$ for ${\displaystyle |x|<1}$.
2. Symmetry: ${\displaystyle W(-x)=W(x)}$.
3. Monotone: ${\displaystyle W(x)}$ is a nonincreasing function for ${\displaystyle x\ge 0}$.
4. Bounded support: ${\displaystyle W(x)=0}$ for ${\displaystyle |x|\ge 1}$.

Asymptotic efficiency of weight functions has been considered by V. A. Epanechnikov (1969) in the context of kernel density estimation; J. Fan (1993) has derived similar results for local regression. They conclude that the quadratic kernel, ${\displaystyle W(x)=1-x^2}$ for ${\displaystyle |x|\le 1}$ has greatest efficiency under a mean-squared-error loss function. See "kernel functions in common use" for more discussion of different kernels and their efficiencies.

Considerations other than MSE are also relevant to the choice of weight function. Smoothness properties of ${\displaystyle W(x)}$ directly affect smoothness of the estimate ${\displaystyle \hat{\mu }(x)}$. In particular, the quadratic kernel is not differentiable at ${\displaystyle x=\pm 1}$, and ${\displaystyle \hat{\mu }(x)}$ is not differentiable as a result. The tri-cube weight function, $${\displaystyle W(x)=(1-|x{|}^3{)}^3;|x|<1}$$ has been used in LOWESS and other local regression software; this combines higher-order differentiability with a high MSE efficiency.

One criticism of weight functions with bounded support is that they can lead to numerical problems (i.e. an unstable or singular design matrix) when fitting in regions with sparse data. For this reason, some authors choose to use the Gaussian kernel, or others with unbounded support.

Choice of fitting criterion

As described above, local regression uses a locally weighted least squares criterion to estimate the regression parameters. This inherits many of the advantages (ease of implementation and interpretation; good properties when errors are normally distributed) and disadvantages (sensitivity to extreme values and outliers; inefficiency when errors have unequal variance or are not normally distributed) usually associated with least squares regression.

These disadvantages can be addressed by replacing the local least-squares estimation by something else. Two such ideas are presented here: local likelihood estimation, which applies local estimation to the generalized linear model, and robust local regression, which localizes methods from robust regression.

Local likelihood estimation

In local likelihood estimation, developed in Tibshirani and Hastie (1987), the observations ${\displaystyle Y_i}$ are assumed to come from a parametric family of distributions, with a known probability density function (or mass function, for discrete data), $${\displaystyle Y_i\sim f(y,\theta (x_i)),}$$ where the parameter function ${\displaystyle \theta (x)}$ is the unknown quantity to be estimated. To estimate ${\displaystyle \theta (x)}$ at a particular point ${\displaystyle x}$, the local likelihood criterion is $${\displaystyle \sum _{i=1}^n{w}_i(x)\log \left[f\left(Y_i,{\beta }_0+{\beta }_1(x_i-x)+\cdots +{\beta }_p{\left(x_i-x\right)}^p\right)\right].}$$ Estimates of the regression coefficients (in, particular, ${\displaystyle {\hat{\beta }}_0}$) are obtained by maximizing the local likelihood criterion, and the local likelihood estimate is $${\displaystyle \hat{\theta }(x)={\hat{\beta }}_0.}$$

When ${\displaystyle f(y,\theta (x))}$ is the normal distribution and ${\displaystyle \theta (x)}$ is the mean function, the local likelihood method reduces to the standard local least-squares regression. For other likelihood families, there is (usually) no closed-form solution for the local likelihood estimate, and iterative procedures such as iteratively reweighted least squares must be used to compute the estimate.

Example (local logistic regression). All response observations are 0 or 1, and the mean function is the "success" probability, ${\displaystyle \mu (x_i)=Pr(Y_i=1|x_i)}$. Since ${\displaystyle \mu (x_i)}$ must be between 0 and 1, a local polynomial model should not be used for ${\displaystyle \mu (x)}$ directly. Insead, the logistic transformation $${\displaystyle \theta (x)=\log \left(\frac{\mu (x)}{1-\mu (x)}\right)}$$ can be used; equivalently, $${\displaystyle \begin{array}{cc}1-\mu (x)& =\frac{1}{1+e^{\theta (x)}};\\ \mu (x)& =\frac{e^{\theta (x)}}{1+e^{\theta (x)}}\end{array}}$$ and the mass function is $${\displaystyle f(Y_i,\theta (x_i))=\frac{e^{Y_i\theta (x_i)}}{1+e^{\theta (x_i)}}.}$$

An asymptotic theory for local likelihood estimation is developed in J. Fan, Nancy E. Heckman and M.P.Wand (1995); the book Loader (1999) discusses many more applications of local likelihood.

Robust local regression

To address the sensitivity to outliers, techniques from robust regression can be employed. In local M-estimation, the local least-squares criterion is replaced by a criterion of the form $${\displaystyle \sum _{i=1}^n{w}_i(x)\rho \left(\frac{Y_i-{\beta }_0-\cdots -{\beta }_p(x_i-x{)}^p}{s}\right)}$$ where ${\displaystyle \rho (\cdot )}$ is a robustness function and ${\displaystyle s}$ is a scale parameter. Discussion of the merits of different choices of robustness function is best left to the robust regression literature. The scale parameter ${\displaystyle s}$ must also be estimated. References for local M-estimation include Katkovnik (1985) and Alexandre Tsybakov (1986).

The robustness iterations in LOWESS and LOESS correspond to the robustness function defined by $${\displaystyle {\rho }^{\prime }(u)=u(1-u^2/6{)}^2;|u|<1}$$ and a robust global estimate of the scale parameter.

If ${\displaystyle \rho (u)=|u|}$, the local ${\displaystyle L_1}$ criterion $${\displaystyle \sum _{i=1}^n{w}_i(x)\left|Y_i-{\beta }_0-\dots -{\beta }_p(x_i-x{)}^p\right|}$$ results; this does not require a scale parameter. When ${\displaystyle p=0}$, this criterion is minimized by a locally weighted median; local ${\displaystyle L_1}$ regression can be interpreted as estimating the median, rather than mean, response. If the loss function is skewed, this becomes local quantile regression. See Keming Yu and M.C. Jones (1998).

A recent alternative is to modify the local regression weights rather than the loss function. Shulman (2025) proposes robust local polynomial regression with similarity kernels, in which the kernel weighting is generalized to incorporate both predictor and response variables; in one version, the local least-squares criterion is reweighted by an estimate of the conditional density ${\displaystyle {\hat{f}}_{Y\mid X}(Y_i\mid X_i)}$, so that observations with low estimated local conditional density are downweighted. This provides robustness to outliers and high-leverage points without the multiple robustness iterations used in methods such as LOWESS and LOESS.

Advantages

As discussed above, the biggest advantage LOESS has over many other methods is the process of fitting a model to the sample data does not begin with the specification of a function. Instead the analyst only has to provide a smoothing parameter value and the degree of the local polynomial. In addition, LOESS is very flexible, making it ideal for modeling complex processes for which no theoretical models exist. These two advantages, combined with the simplicity of the method, make LOESS one of the most attractive of the modern regression methods for applications that fit the general framework of least squares regression but which have a complex deterministic structure.

Although it is less obvious than for some of the other methods related to linear least squares regression, LOESS also accrues most of the benefits typically shared by those procedures. The most important of those is the theory for computing uncertainties for prediction and calibration. Many other tests and procedures used for validation of least squares models can also be extended to LOESS models .

Disadvantages

LOESS makes less efficient use of data than other least squares methods. It requires fairly large, densely sampled data sets in order to produce good models. This is because LOESS relies on the local data structure when performing the local fitting. Thus, LOESS provides less complex data analysis in exchange for greater experimental costs.

Another disadvantage of LOESS is the fact that it does not produce a regression function that is easily represented by a mathematical formula. This can make it difficult to transfer the results of an analysis to other people. In order to transfer the regression function to another person, they would need the data set and software for LOESS calculations. In nonlinear regression, on the other hand, it is only necessary to write down a functional form in order to provide estimates of the unknown parameters and the estimated uncertainty. Depending on the application, this could be either a major or a minor drawback to using LOESS. In particular, the simple form of LOESS can not be used for mechanistic modelling where fitted parameters specify particular physical properties of a system.

Finally, as discussed above, LOESS is a computationally intensive method (with the exception of evenly spaced data, where the regression can then be phrased as a non-causal finite impulse response filter). LOESS is also prone to the effects of outliers in the data set, like other least squares methods. There is an iterative, robust version of LOESS [Cleveland (1979)] that can be used to reduce LOESS' sensitivity to outliers, but too many extreme outliers can still overcome even the robust method; other non-iterative robust local regression methods have also been proposed.