10.5 Nonparametric regression
A polynomial imposes one global form on the whole range of \(x\). Nonparametric regression drops the global form: it estimates the mean at each point \(x_0\) from the observations near \(x_0\) only. Let \(K\) be a kernel, a non-negative function that is largest at \(0\) and decreases away from it, and let \(h > 0\) be the bandwidth. The weight of observation \(i\) at \(x_0\) is \[ w_i(x_0) = K\Bigl (\frac {x_i - x_0}{h}\Bigr ) , \] and the kernel estimate (the Nadaraya–Watson estimate) of the mean at \(x_0\) is the weighted average \[ \tilde y(x_0) = \frac {\sum _i w_i(x_0)\,y_i}{\sum _i w_i(x_0)} . \] A common choice is the triangular kernel \(K(t) = \max (1 - |t|, 0)\), which gives weight only to the observations within \(h\) of \(x_0\).
- (i)
- The kernel estimate \(\tilde y(x_0)\) is the weighted least squares estimate of a constant: it minimises \(\sum _i w_i(x_0)(y_i - a)^2\) over \(a\).
- (ii)
- The local linear estimate, \(\tilde y(x_0) = \hat a\) where \((\hat a, \hat b)\) minimises \(\sum _i w_i(x_0)\bigl (y_i - a - b(x_i - x_0)\bigr )^2\), reproduces any straight line exactly: if \(y_i = c + dx_i\) for every \(i\), then \(\tilde y(x_0) = c + dx_0\). The kernel estimate does not, in general.
- (iii)
- Both estimates are linear in \(\mathbf {y}\): the vector of estimates at the observed \(x_i\) is \(\tilde {\mathbf {y}} = \mathbf {Sy}\) for a matrix \(\mathbf {S}\) that depends only on the \(x_i\), \(K\) and \(h\).
Proof. (i) The derivative of \(\sum w_i(y_i - a)^2\) is \(-2\sum w_i(y_i - a)\), zero at \(a = \sum w_iy_i/\sum w_i\), and the second derivative \(2\sum w_i\) is positive.
(ii) The local linear estimate is weighted least squares for the model with columns \(\mathbf {1}\) and \(\mathbf {x} - x_0\mathbf {1}\) (Theorem 6.4.1), with \(\hat {\boldsymbol \beta } = (\mathbf {X}^\top \mathbf {WX})^{-1}\mathbf {X}^\top \mathbf {Wy}\). If \(\mathbf {y} = (c + dx_0)\mathbf {1} + d(\mathbf {x} - x_0\mathbf {1})\), then \(\mathbf {y} = \mathbf {X}\boldsymbol \theta \) with \(\boldsymbol \theta = (c + dx_0, d)^\top \), so \(\hat {\boldsymbol \beta } = (\mathbf {X}^\top \mathbf {WX})^{-1}\mathbf {X}^\top \mathbf {WX}\boldsymbol \theta = \boldsymbol \theta \), and \(\hat a = c + dx_0\). For the kernel estimate, take \(x_0\) at the smallest \(x_i\) and \(d > 0\): every observation with positive weight has \(y_i \ge c + dx_0\), with strict inequality for those to the right of \(x_0\), so the weighted average exceeds \(c + dx_0\) as soon as any of them has positive weight.
(iii) In both cases \(\tilde y(x_0)\) is a fixed linear combination of the \(y_i\): \(\sum _i w_i(x_0)y_i/\sum _i w_i(x_0)\), or the first row of \((\mathbf {X}^\top \mathbf {WX})^{-1}\mathbf {X}^\top \mathbf {W}\) applied to \(\mathbf {y}\). Stacking these rows for \(x_0 = x_1, \dots , x_n\) gives \(\mathbf {S}\). □
The matrix \(\mathbf {S}\) plays the part of the hat matrix \(\mathbf {H}\). For a linear model, by Lemma 3.1.1, \[ \operatorname {tr}\mathbf {H} = \operatorname {tr}\bigl ((\mathbf {X}^\top \mathbf {X})^{-1}\mathbf {X}^\top \mathbf {X}\bigr ) = \operatorname {tr}\mathbf {I}_p = p , \] the number of parameters, so \(\operatorname {tr}\mathbf {S}\) is called the effective number of parameters of the smoother. The bandwidth controls the trade-off: a small \(h\) follows the data closely (small bias, large variance, large \(\operatorname {tr}\mathbf {S}\)), a large \(h\) averages over many points (large bias at peaks and troughs, small variance). Part (ii) explains why local linear fitting is preferred near the ends of the data, where the kernel estimate is pulled towards the interior. Fitting a local line with tricube weights, and down-weighting observations with large residuals, is the loess method used by most software. Nonparametric fits are useful for exploring the shape of a relationship and for checking a parametric model, but they give no formula, and they are unreliable near the ends of the data and wherever the data are sparse.
Example 10.5.2. For the irrigation data of Example 10.1.4, use the triangular kernel with \(h = 2\) to estimate the mean yield at \(x_0 = 4\) and at \(x_0 = 1\), by the kernel and the local linear methods, and find \(\operatorname {tr}\mathbf {S}\) for the kernel estimate.
Solution. At \(x_0 = 4\) the weights \(1 - |x_i - 4|/2\) are \(\frac 12, 1, \frac 12\) at \(x = 3, 4, 5\) and zero elsewhere, so \[ \tilde y(4) = \frac {\frac 12(23) + 26 + \frac 12(27)}{2} = 25.5 . \] The weights are symmetric about \(x_0\), so \(\sum w_i(x_i - x_0) = 0\); the local linear normal equations then decouple, and the local linear estimate is also \(25.5\). The quadratic fit gave \(26.14\) here: averaging across the peak of a concave curve pulls the smooth estimate down, the bias of a large bandwidth.
At \(x_0 = 1\) the weights are \(1\) at \(x = 1\) and \(\frac 12\) at \(x = 2\). The kernel estimate is \[ \tilde y(1) = \frac {8 + \frac 12(16)}{1.5} = 10.67 , \] well above the observation \(8\) and the quadratic fit \(8.14\), because every neighbour lies to the right, uphill. The local linear fit uses the same two points; a line through two points fits them exactly, so \(\tilde y(1) = 8\), with local slope \(8\). This is Theorem 10.5.1(ii) at work.
For \(\operatorname {tr}\mathbf {S}\) we need the weight of \(y_i\) in its own estimate. At \(x_0 = 1\) and \(x_0 = 7\) it is \(1/1.5 = \frac 23\); at each interior point the weights are \(\frac 12, 1, \frac 12\) with total \(2\), so it is \(\frac 12\). Then \[ \operatorname {tr}\mathbf {S} = 2\Bigl (\frac 23\Bigr ) + 5\Bigl (\frac 12\Bigr ) = 3.83 , \] a little more than the three parameters of the quadratic.
Questions on this section
Stuck on something here? Ask below and it stays attached to this topic.