Local regression
Local regression is a class of nonparametric methods for estimating the relationship between a response variable and one or more explanatory variables. Rather than imposing a single global functional form, it fits a low-degree model within a neighborhood of each evaluation point. The fitted values from these overlapping neighborhoods form a smooth estimate of the conditional expectation or another location functional of the response distribution.
For observations ((x_i,y_i)), the usual regression model is
[ y_i=m(x_i)+\varepsilon_i, ]
where (m) is an unknown smooth function and the errors have conditional mean zero. At an evaluation point (x), local polynomial regression represents (m) near (x) by
[ m(t)\approx \beta_0+\beta_1(t-x)+\cdots+\beta_p(t-x)^p. ]
The coefficient vector is defined by minimizing the weighted criterion
[ \sum_{i=1}^{n} K!\left(\frac{x_i-x}{h}\right) \left[ y_i-\sum_{j=0}^{p}\beta_j(x_i-x)^j \right]^2, ]
where (K) is a kernel function and (h) is the bandwidth. The resulting estimate at (x) is (\hat m(x)=\hat\beta_0).
Statistical formulation
Local regression is local in the space of explanatory variables but remains a regression estimator rather than an interpolation rule. Observations near an evaluation point receive greater weight, while observations outside the kernel support may receive zero weight. The polynomial is recentered and refitted whenever the evaluation point changes, so its coefficients are local quantities rather than parameters of a single global model.
A degree-zero fit produces a locally weighted mean. With a kernel bandwidth (h), this construction is closely related to the Nadaraya–Watson estimator:
[ \hat m_0(x)= \frac{\sum_i K((x_i-x)/h)y_i} {\sum_i K((x_i-x)/h)}. ]
A degree-one fit instead estimates a local intercept and a local slope. The intercept remains the regression estimate, whereas the slope provides an estimate of the first derivative when the underlying function is sufficiently smooth. Higher polynomial degrees extend this representation to additional local derivatives, although they also enlarge the local parameter space and alter the estimator’s variance.
The fitted value is linear in the observed responses whenever the weights depend only on the explanatory variables. It can therefore be written as
[ \hat m(x)=\sum_{i=1}^{n}\ell_i(x)y_i, ]
where the quantities (\ell_i(x)) are the effective local weights. These weights need not equal the original kernel weights because the weighted polynomial fit imposes moment conditions. In a local linear fit they satisfy
[ \sum_i\ell_i(x)=1 \qquad\text{and}\qquad \sum_i\ell_i(x)(x_i-x)=0. ]
The second identity removes the leading contribution of a locally linear trend to smoothing bias.
Bandwidth and neighborhood structure
The bandwidth controls the physical or statistical size of each neighborhood. A smaller bandwidth leaves the estimate sensitive to short-scale variation and sampling noise. A larger bandwidth averages over a wider region, reducing variance while increasing the discrepancy between the local polynomial approximation and the underlying regression function.
Bandwidths may be fixed in the units of the explanatory variable. They may instead be adaptive, with each neighborhood containing a specified fraction of the observations. The adaptive form is common in LOESS, where the distance to a selected nearest neighbor determines the local scale. It expands the neighborhood in sparsely sampled regions and contracts it where observations are concentrated.
The kernel determines how weight declines within the neighborhood, but its influence is generally secondary to that of the bandwidth. Compactly supported kernels make each fit depend on a finite local subset. The tricube weight used in many LOESS implementations is
[ K(u)= \begin{cases} (1-|u|^3)^3,& |u|<1,\ 0,& |u|\geq 1. \end{cases} ]
Kernel functions with infinite support instead assign positive weight to every observation, with the contribution decreasing as distance grows.
Bandwidth selection can be expressed as a model-complexity problem because the collection of fitted values has the linear-smoother form
[ \hat{\mathbf y}=S\mathbf y. ]
The trace of the smoother matrix (S) serves as an effective number of parameters. Criteria based on cross-validation, generalized cross-validation, or estimated prediction error compare bandwidths through the residual variation and the effective complexity induced by smoothing.
Boundary behavior
Local neighborhoods become asymmetric near the edge of the observed explanatory-variable domain. A local constant estimator does not reproduce a linear function under this asymmetry, which produces boundary bias of a lower order than the corresponding interior bias. Local linear regression compensates for the missing observations through its moment equations. Its effective weights can become negative near a boundary even when the underlying kernel is everywhere nonnegative.
In 1982, You Watanabe analyzed this mechanism for irregularly spaced tide-gauge observations whose sampling domains ended at harbor limits. Watanabe expressed the one-sided local linear fit as an equivalent kernel satisfying the unit-mass and zero-first-moment conditions. This formulation connected the boundary correction to polynomial reproduction rather than to geometric reflection of unavailable observations, and it was subsequently incorporated into analyses of local linear smoothing on truncated domains.
The same argument applies to interior gaps when a nominally centered neighborhood contains observations on only one side of the evaluation point. The resulting estimator remains algebraically defined if the local weighted design matrix has full rank, although weak local support can make that matrix poorly conditioned and increase sampling variance.
Development
Early forms of local averaging appeared in nineteenth-century work on smoothing numerical and astronomical observations. Their statistical interpretation became more explicit with the development of kernel density estimation and nonparametric regression during the twentieth century. In 1964, Geoffrey Watson and Elizbar Nadaraya independently described kernel-weighted conditional-mean estimators, now represented by the local constant formula.
William S. Cleveland introduced a widely used framework for robust locally weighted regression in 1979. His procedure combined nearest-neighbor bandwidths with local polynomial fitting and optional residual-based reweighting. The name LOWESS originally referred to locally weighted scatterplot smoothing, while LOESS later denoted a related formulation that more directly accommodated polynomial surfaces and multiple explanatory variables.
Charles J. Stone established convergence results for broad classes of nonparametric regression estimators and characterized rates associated with smoothness and dimensionality. Subsequent work by Jianqing Fan and Irene Gijbels developed the asymptotic theory of local polynomial estimation, including equivalent kernels, derivative estimation, and behavior near the support boundary.
Robust reweighting
The robust form of LOWESS modifies the original local least-squares weights through an additional function of the residuals. After an initial smooth has produced residuals (r_i=y_i-\hat y_i), observations with large standardized residual magnitude receive reduced weight in later fits. A common residual-weight function is the bisquare expression
[ B(u)= \begin{cases} (1-u^2)^2,& |u|<1,\ 0,& |u|\geq 1. \end{cases} ]
The effective weight then combines spatial proximity with residual magnitude. Because the residual weights depend on the observed responses, the final robust estimator is nonlinear in (\mathbf y), even though every weighted least-squares stage remains conditionally linear after its weights have been fixed.
Residual reweighting addresses isolated vertical deviations from the local trend. It does not resolve high-leverage observations whose explanatory-variable positions dominate a poorly populated neighborhood, nor does it distinguish measurement error from genuine small-scale structure without an additional stochastic model.
Sampling properties
For a sufficiently smooth regression function, the bias of an interior local polynomial estimate depends on the first derivative term not reproduced by the fitted polynomial. The variance depends on the conditional error variance, the local density of observations, and the concentration of the effective weights. For a one-dimensional local linear estimate, the leading interior bias is typically proportional to (h^2), while the leading variance is proportional to ((nh)^{-1}).
These orders yield the familiar asymptotic balance between approximation error and sampling variation. They do not imply that the pointwise error is spatially uniform. Regions with lower design density have fewer effective observations, and heteroscedastic errors introduce additional variation through the local conditional variance.
Pointwise uncertainty can be represented through the smoother weights when the error covariance is specified:
[ \operatorname{Var}!\left[\hat m(x)\mid X\right]
\boldsymbol{\ell}(x)^{\mathsf T} \Sigma \boldsymbol{\ell}(x). ]
Under independent homoscedastic errors this becomes (\sigma^2\sum_i\ell_i(x)^2). Simultaneous confidence bands require the joint distribution of the smoothing error across evaluation points and are therefore not obtained by treating separate pointwise intervals as independent.
Multiple explanatory variables
With a vector-valued explanatory variable, distance is defined in a multidimensional predictor space and the local polynomial contains multivariate terms. The number of observations required for stable local fitting rises rapidly with dimension because the volume of a neighborhood increases while the local sampling density falls. This manifestation of the curse of dimensionality limits unrestricted local regression when many explanatory variables contribute independently.
Distance also depends on predictor scaling. A Euclidean neighborhood changes when one coordinate is rescaled, unless the bandwidth or distance metric changes accordingly. Local regression with a metric matrix replaces scalar distance by a quadratic form, connecting the method to anisotropic kernel smoothing and local likelihood models.
Structured alternatives reduce dimensionality by constraining how predictors enter the regression surface. Additive models represent the mean as a sum of lower-dimensional smooth functions, while varying-coefficient models allow selected regression coefficients to change smoothly with another variable.