![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to xchap2.tex CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10dc / vol2 / index|
Return to xchap2.tex CVS log | Up to [Research Unix] / researchv10dc / vol2 / index |
1.1 root 1: \ifnum\pageno=1 \input macros.tex \pageno=200\immediate\openout\inx=chap2.index\makecontents\draft{\date (Rob)} \fi
2: %\proofmodefalse
3: \def\gsmooth{{s}}
4: \def\vgsmooth{{\vec s}}
5: \def\symNN{N^S}
6: \def\NN{N}
7: \def\pepscat{1.1}
8: \def\pepadd{1.2}
9: \def\X{\bf X}
10: \input try.tex
11: \chapter 2 {Smoothing}
12: \index{smoothing}%
13: \index{what is a smoother?}%
14: \medskip
15: \Section{What is a smoother?}
16: A smoother is a tool for summarizing the dependence of a response
17: measurement $Y$ on one or more predictor measurements $X_1,\ldots,
18: X_p$.
19: \index{ smoother}%
20: It produces an estimate of this dependence that is less variable than $Y$
21: itself; hence the name {\sl smoother}.
22: \index{ nonparametric nature}%
23: \index{rigid form}%
24: An important property of a smoother is its {\sl nonparametric} nature; that is,
25: it doesn't assume a rigid
26: form for the dependence.
27: \index{parametric regression}%
28: For this reason, a smoother is
29: often referred to as a
30: tool for nonparametric regression.
31: \index{running mean}%
32: \index{moving average}%
33: The running mean (or moving average) is a simple example of a smoother,
34: while a regression line is not thought of as a smoother because of its rigid
35: parametric form.
36: \figname{\spansizes}
37: \pageinsert
38: \vbox to \vsize{
39: \centerline{\psfig{file=figures/spansizes.p,width=4in}}
40: \figurecaption{\spansizes}{Bias-variance tradeoff. A sample of size 80
41: was generated from $Y=sin(X)+\varepsilon$ with $X\sim N(0,1)$ and
42: \index{smoothing}%
43: \index{smoothing}%
44: $\varepsilon\sim N(0,1/9)$. The top two panels represent one realization
45: of the model with a smooth included; panel (a) uses a small neighbourhood
46: and thus
47: does little smoothing, while (b) uses a large neighbourhood and does a
48: lot of smoothing. In terms of approximate degrees of freedom
49: (a concept introduced in Chapter~3), the smoother in (a) uses
50: $\df =15$, whereas for (b) $\df=4$.
51: The $\sin$ function is included in both panels.
52: The lower two panels show the bias-variance behaviour of this smoother.
53: \index{fitted smooth}%
54: The light shaded regions show
55: twice the pointwise standard errors of the fitted smooth values, and
56: thus reflect the variance of the smooths under repeated sampling of
57: 80 pairs from the model. The darker shaded region shows the bias,
58: which is the difference between the generating curve and the average
59: value of the smooth curves under repeated sampling. In panel (c),
60: the bias is virtually nonexistent, while the standard-error bands are
61: quite wide. In panel (d), there is a considerable cost in bias
62: \index{ smooth}%
63: for a moderate decrease in variance.}*
64: \vfill
65: }
66: \endinsert
67: \par\figname{\allsmooths}
68: We call the estimate produced by a smoother a {\sl smooth}.
69: \index{smoothing}%
70: \index{single predictor}%
71: \index{scatterplot smoothing}%
72: The most important example of smoothing involves a single predictor
73: setting, usually referred to as {\sl scatterplot
74: smoothing}.
75:
76: Smoothers have two main uses.
77: The first use is for description.
78: \index{scatterplot smoother}%
79: For example, a scatterplot smoother can be used to enhance the
80: visual appearance of the scatterplot of $Y$ vs $X$,
81: to help our eyes pick out the trend in the plot.
82: Figure~\pepscat\ in Chapter~1 shows a plot
83: of \name{log(C-peptide)} versus \name{age} for the diabetes data.
84: \index{scatterplot smoother}%
85: It seems that \name{log(C-peptide)} has a strong dependence on \name{age} and a scatterplot smoother
86: will be helpful in describing this relationship.
87: \index{scatterplot smooth}%
88: Figure~\allsmooths\ shows a number of scatterplot smooths of these data; each described in turn in this chapter.
89: \index{building block}%
90: The second use of a smoother is as an estimator of the dependence of $\ev(Y)$ on
91: the predictors, and thus as a building block for the estimation of additive
92: models, discussed in the remainder of the book.
93:
94:
95: In this chapter we give a brief overview of some useful smoothers.
96: \index{scatterplot smoother}%
97: Most of our discussion concerns scatterplot smoothers; at the end of the
98: chapter we give a brief description of multiple-predictor smoothers.
99: \index{ categorical}%
100:
101: The simplest smoother occurs in the case of a {\sl categorical}
102: predictor,
103: for example, \name{sex} (male, female) or \name{colour} (red, blue, green etc.)
104: To smooth $Y$ we can simply average the values of $Y$ in each category.
105: \index{scatterplot smooth}%
106: This satisfies our requirements for a scatterplot smooth: it captures the
107: trend of $Y$ on $X$ and is smoother than the $Y$-values themselves.
108: \index{ smoothing}%
109: \index{smoothing}%
110: While the reader might not normally think of this as {\sl smoothing}, this simple
111: averaging process is the conceptual basis for
112: smoothing in the most general setting, that of an ordered (noncategorical)
113: predictor.
114: The problem here is often the lack of replicates at each
115: predictor value.
116: \index{ local averaging}%
117: \index{target value}%
118: Most smoothers attempt to mimic category averaging through
119: {\sl local averaging}, that is,
120: averaging the
121: $Y$-values of observations having predictor values close to a target
122: value.
123: \index{ neighbourhoods around}%
124: \index{target value}%
125: The averaging is done in {\sl neighbourhoods} around the target value.
126: \index{scatterplot smoothing}%
127: There are two main decisions to be made in scatterplot
128: smoothing:\smallskip
129: {\parindent 20pt
130: \item{(i)}how to average the
131: response values in each neighbourhood, and
132: \item{(ii)} how big to take the neighbourhoods.
133:
134: }\smallskip
135: \pageinsert
136: \vbox to \vsize{
137: \centerline{\psfig{file=figures/allsmooths.p,width=\hsize}}
138: \figurecaption{\allsmooths}{A variety of smooths applied to the scatterplot of \name{log(C-peptide)}
139: \index{smoothing parameter}%
140: versus \name{age}. In each case the smoothing parameter was chosen so
141: that degrees of freedom in the fit is about five. The linear fit of
142: course has only two degrees of freedom. The parametric bin and quartic
143: \index{three interior knots}%
144: \index{smoothing parameter}%
145: polynomial fits have five parameters. The (cubic) regression spline has one interior knot at the median and hence five parameters. The natural cubic regression spline has three interior knots and two endpoint constraints, which also results in five parameters. For the other smooths, the smoothing parameter $\lambda$ was chosen such that
146: \index{ brand}%
147: $\df\approx 5$, where \df\ is a measure of the approximate degrees of freedom (a concept introduced in Chapter~3).}*
148: \vfill
149: }
150: \endinsert
151: \par
152: The question of how to average within a neighbourhood is really the
153: question of which {\sl brand} of smoother to use, because smoothers
154: differ mainly in their method of averaging.
155: In this chapter we describe a number of different smoothers and
156: compare them informally.
157: Formal recommendations on how to choose among smoothers are
158: difficult to make because
159: few systematic comparisons have been made so far in the
160: literature.
161:
162: The question of how big to make the neighbourhoods is discussed in
163: the next chapter.
164: The issue underlying this question is very important, however, so we briefly
165: discuss it here.
166: \index{ smoothing parameter}%
167: The size of the neighbourhood is typically expressed in terms of an
168: adjustable {\sl smoothing parameter}.
169: Intuitively, large neighbourhoods will produce an estimate with low
170: variance but potentially high bias, and conversely for small
171: neighbourhoods (Fig.~\spansizes).
172: \index{ fundamental tradeoff between bias and variance}%
173: \index{smoothing parameter}%
174: Thus there is a {\sl fundamental tradeoff between bias and variance},
175: governed by the smoothing parameter.
176: This issue is exactly analogous to the question of how many predictors to put
177: in a regression equation.
178: \index{smoothing parameter}%
179: In the
180: next
181: chapter we discuss this important issue and the practical
182: question of how to choose the smoothing parameter, based on the data,
183: to trade bias against variance in an optimal way.
184: \index{fitted smooth}%
185: We also discuss some other important topics, such as
186: linear and nonlinear smoothers, the
187: incorporation of observation weights into --- and inference
188: for --- the
189: fitted smooth.
190:
191: \index{ robustified}%
192: All the smoothers we shall describe can easily be {\em robustified} by replacing the averaging or least-squares operation by a more robust procedure.
193: Details can be found in Chapters~3 and 9.
194:
195: \index{bibliographic notes}%
196:
197: Bibliographic notes on smoothers are given at the end
198: of the next chapter.
199:
200: \index{scatterplot smoothing: definition}%
201: \index{scatterplot smoothing}%
202: \Sectionskip
203: \Section{Scatterplot smoothing: definition}
204: \Mark{SCATTERPLOT SMOOTHING}
205: Suppose we have response measurements $\vec y=(y_1,\ldots,y_n)^T$ at design
206: points $\vec x=(x_1,\ldots,
207: x_n)^T$. We assume that each of $\vec y$ and $\vec x$ represent
208: measurements of variables $Y$ and $X$.
209: In most cases it is useful to think of
210: $Y$, and sometimes $X$, as having been generated by some random mechanism, but this is not
211: necessary for the discussion here.
212: In particular we don't need to assume
213: that the pairs $(x_i,y_i)$ are a random sample from some joint distribution.
214: For
215: example, the $X$-values might be preset dose levels of a drug.
216:
217:
218:
219: Since $Y$ and $X$ are noncategorical we don't expect to find
220: many replicates at any given value of $X$.
221: For convenience we assume
222: that the data are sorted by $X$
223: and for the present discussion that there are no tied $X$-values, so that
224: $x_1 < \cdots < x_n$.
225: A trivial remedy when there are ties is to use weighted smoothers, which we
226: \index{scatterplot smooth}%
227: discuss in Chapter~3.
228:
229:
230: We denote a scatterplot smooth of $\vec y$ against $\vec x$ by $\gsmooth=\smooth(\vec y\given \vec x)$, where $\gsmooth=\gsmooth(x)$ is a function of $X$.
231: Usually the recipe that defines $\gsmooth(x_0)$, which is the function $\smooth(\vec y\given \vec x)$ evaluated at $x_0$, will be defined for all $x_0$, but at other times is defined only at $x_1,\ldots,x_n$, the sample values of $X$.
232: In this latter case some kind of interpolation is necessary in order to
233: obtain estimates at other $X$-values.
234:
235: \index{parametric regression}%
236: \Sectionskip
237: \Section{Parametric Regression}
238: A regression line, estimated for example by least-squares,
239: provides an estimate of the dependence of $\ev(Y)$ on $X$.
240: \index{rigid form}%
241: It does so, however,
242: by assuming a rigid form for this dependence,
243: and thus it
244: may or may not be appropriate for a given set of data.
245: If the dependence is linear or close to it, the regression line provides a
246: concise and useful summary.
247: For the data in Fig.~\allsmooths, however, the regression line is clearly inappropriate and creates a misleading impression.
248:
249: \index{ infinitely smooth}%
250: \index{scatterplot smoother}%
251: In a sense the regression line is an {\sl infinitely smooth} function, and
252: not surprisingly many of the scatterplot smoothers that we
253: discuss approach the linear regression line as one extreme.
254: The other extreme
255: is usually some function that interpolates the data.
256:
257:
258: Other nonlocal parametric fits, such as polynomial regression estimates, share
259: the same pros and cons as the regression line.
260: They are useful
261: if they are appropriate for the data at hand but potentially misleading otherwise.
262: \index{parametric fitting}%
263: \index{scatterplot smoothing}%
264: Although parametric fitting certainly isn't the solution to the scatterplot
265: smoothing problem,
266: we will still make use of it for
267: comparative purposes and for possible summaries of
268: \index{regression splines}%
269: \index{parametric fitting}%
270: smooths estimated by other methods.
271: Regression splines are a less rigid
272: \index{bin smoother}%
273: \index{bin smoother}%
274: form of parametric fitting which are closer in spirit to a smoother; we discuss them later in this chapter.
275: \Sectionskip
276: \Section{Bin smoothers}
277: A bin smoother mimics a categorical smoother by
278: partitioning the predictor values into a number of disjoint and exhaustive regions,
279: then averaging the response in each region.
280: \index{data points}%
281: Formally, we choose cutpoints
282: $c_0 < \cdots < c_K$
283: where $c_0=-\infinity$ and $c_K=\infinity$,
284: and define
285: $$R_i=\{j;c_i\leq x_j < c_{i+1}\};\qquad i=1,\ldots, K$$
286: the indices of
287: the data points in each region.
288: Then $\gsmooth=\smooth(\vec y\given \vec x)$ is given by $\gsmooth(x_0)=\ave_{j\in R_i} (y_j)$ if
289: $x_0 \in R_i$.
290: Typically one chooses five regions, for example, and picks the cutoff points
291: so that
292: there is approximately an equal number of points in each region.
293: The bin smooth shown in Fig.~\allsmooths\ was constructed in this way and illustrates
294: the limitation of the method.
295: The estimate is not very smooth because it jumps at each cut point.
296: Unless we have a prior reason to believe such discontinuities exist, this
297: estimate is not very appealing.
298: One way to make it more smooth is through the use of
299: overlapping regions, as we'll see in the next section.
300: \index{running-mean and running-lines smoothers}%
301: \index{target value}%
302: \Sectionskip
303: \Section{Running-mean and running-lines smoothers}
304: Assume that our target value $x_0$ equals one of the $x_j$s, say $x_i$. If we had replicates at $x_i$,
305: we could simply use the average of the $Y$-values at $x_i$ as our estimate
306: $\gsmooth(x_i)$.
307: Now
308: \index{ close}%
309: we are assuming that we don't have replicates, so instead we can average $Y$-values corresponding to $X$-values close to $x_i$.
310: How do we pick points that
311: are {\sl close} to $x_i$?
312: A simple way is to choose $x_i$ itself, as well as the $k$ points to the left and $k$ points to the
313: right of $x_i$ that are closest in $X$-value to $x_i$.
314: \index{ symmetric nearest neighbourhood}%
315: This is called a {\sl symmetric nearest neighbourhood} and
316: the indices of these points
317: will be denoted by
318: \index{ running mean}%
319: $\symNN(x_i)$. Then we could define
320: the {\sl running mean}
321: $$\gsmooth(x_i)=\ave_{j\in
322: \symNN(x_i)} (y_j).\eqn{\runm}$$
323: If it is not possible to take $k$ points to the left or right of
324: $x_i$, we take as many as we can.
325: \index{symmetric nearest neighbourhood}%
326: A formal definition of a symmetric nearest neighbourhood is
327: $$\symNN(x_i)= \{ \max(i-k ,1),\ldots,i-1,i,i+1,\ldots,\min(i+ k
328: ,n)\}.\eqn{\neigh}$$
329:
330: \index{symmetric nearest neighbours}%
331: \index{target point}%
332: It is not obvious how to define the symmetric nearest neighbours at target points $x_0$ other than the $x_i$ in the sample.
333: We could simply interpolate linearly
334: between the fit of the two values of $X$ in the sample adjacent to $x_0$.
335: \index{ nearest neighbourhood}%
336: Alternatively we could ignore symmetry and take the $r$ closest points to $x_0$, regardless
337: of which side they are on; this is called a {\sl nearest neighbourhood}.
338: This handles arbitrary $x_0$ in a simple and clean way.
339: The pros and cons of these two types of neighbourhoods are discussed later.
340: \figname{\rlsmooth}
341: \midinsert
342: \vbox{
343: \centerline{\psfig{file=figures/rlsmooth.p,width=3in}}
344: \index{symmetric nearest neighbourhood}%
345: \vskip .05in
346: \index{target point}%
347: \index{target point}%
348: \index{running mean}%
349: \index{running mean}%
350: \figurecaption{\rlsmooth}{The running-line fit computes linear fits locally in the symmetric nearest neighbourhoods. The fitted value is the point on the line at the target point of the neighbourhood. Typically in the interior the target point is close to the neighbourhood mean, so the running-line behaves like a running mean. On the boundaries the neighbourhoods become asymmetric and the slope corrects the first-order bias of the running mean.
351: \index{ moving average}%
352: }*
353: }
354: \endinsert
355: \par
356: This simple smoother
357: is
358: also called a {\sl moving average}, and is popular for evenly-spaced time-series data.
359: Although it is valuable for theoretical calculation because of its simplicity,
360: in practice it does not work very well. It tends to be so wiggly that it hardly deserves the name {\sl smoother.}
361: Apart from looking unpleasant, it tends to flatten
362: out trends near the endpoints and hence can be severely biased.
363: \index{running mean}%
364: For example, Fig.~\allsmooths\ shows a running mean smooth with $2k+1=11$
365: or 25\% of the 43 observations.
366:
367:
368: \index{running mean}%
369: \index{least-squares line}%
370: A simple generalization of the running mean alleviates the bias problem: we compute
371: \index{data points}%
372: a least-squares line instead of a mean in each neighbourhood.
373: The {\sl running
374: lines smoother} is defined by
375: $$\gsmooth(x_0)=\hat\alpha(x_0) +\hat
376: \beta(x_0) x_0\eqn{\runl}$$
377: where $\hat\alpha(x_0)$ and $ \hat \beta(x_0) $ are the least-squares
378: estimates for the data points in $\symNN(x_0)$.
379: Figure~\rlsmooth\ shows the local line computed in two different
380: neighbourhoods, one in the interior, the other near the boundary. As we
381: would expect, the fit in the interior is dominated largely by the mean and the slope
382: plays a small role, whereas near the boundary the slope is important for
383: picking up the trend in the asymmetric neighbourhood.
384: \index{running-lines smooth}%
385: A running-lines smooth with a neighbourhood size of $2k+1=13$ or 30\% for the diabetes data, is shown in Fig.~\allsmooths.
386: It seems to capture the trend in the data quite nicely
387: but is still somewhat jagged.
388:
389:
390: \index{running-lines smooth}%
391: The parameter $k$ controls the appearance of the running-lines smooth.
392: Large values of $k$ tend to produce smoother curves while small values
393: tend to produce more jagged curves.
394: It is more convenient to think not in terms of $k$ but instead in terms of
395: $w=(2k+1)/n$, the proportion of points in each neighbourhood, called the {\it
396: span}.
397: We denote by $[\symNN(x_i)]$ the number of points in $\symNN(x_i)$.
398: \index{running-lines smooth}%
399: \index{least-squares line}%
400: In the extreme case, if $w= 2$, so that each neighbourhood contains all of
401: the data (note that $w=1$ won't work because of endpoint effects), the running-lines smooth is the least-squares line.
402: On the other hand, if $w=1/n$, each neighbourhood consists only of one data point
403: and hence the smoother interpolates the data.
404: In the next chapter we discuss the quantitative effects of varying the span and a
405: data-based criterion for choosing it.
406:
407:
408: \index{running-lines smoother}%
409: The running-lines smoother can be computed in O(n) operations (Exercise~2.2),
410: and there is a simple, effective method for span selection.
411: \index{building block}%
412: These can be important advantages when it is used as a building block in the
413: iterative algorithms discussed later in this book.
414: \index{smoothing}%
415: On the other hand, it tends to produce curves that are quite jagged so that
416: \index{running-lines smooth}%
417: \index{weighted least-squares}%
418: a second stage of smoothing might be necessary.
419:
420:
421:
422: One way to improve the appearance of the running-lines smooth is through
423: the use of a {\sl weighted} least-squares fit in each
424: \index{running-lines smoother}%
425: neighbourhood.
426: The running-lines smoother can produce jagged output because points in
427: a given neighbourhood are given equal (nonzero) weight while points
428: outside of the neighbourhood are given zero weight.
429: Thus as the neighbourhoods move from left to right, there are discrete
430: changes in the weight given to the leftmost and rightmost points.
431: We can alleviate this problem by
432: giving the highest weight to $x_i$, and weights smoothly decreasing as
433: we move further away from $x_i$.
434: \index{locally-weighted running-lines}%
435: \index{loess}%
436: ^{Cleveland's (1979)} implementation of a locally-weighted running-lines smoother, {\sl loess},
437: is popular and (essentially) works in this way.
438: A full description of locally-weighted running-line smoothers is given later in this Chapter; Fig.~\allsmooths\ shows one applied to the diabetes data.
439: \index{unweighted running-lines smoother}%
440: It is smoother than the unweighted running-lines smoother.
441: \index{running-lines smoother}%
442: The cost of using smooth weights is a computational one, as the
443: shortcuts available for the unweighted running-lines
444: \index{running-lines smoother}%
445: smoother (Exercise~2.2) do not work for locally
446: weighted running-lines smoothers.
447: \index{unweighted running-lines smoother}%
448: \index{running-lines smoother}%
449: Henceforth we will refer to the unweighted running-lines smoother
450: simply as the running-lines smoother.
451: \index{kernel smoother}%
452: \index{kernel smoother}%
453: \index{kernel}%
454: \index{target value}%
455: \Sectionskip
456: \Section{Kernel smoothers}
457: A kernel smoother uses an explicitly defined set of local weights, defined by the {\sl kernel}, to
458: produce the estimate at each target value.
459: \index{kernel smoother}%
460: \index{target point}%
461: Usually a kernel smoother uses weights that decrease in a smooth fashion
462: as one moves away from the target point.
463: \index{kernel smoother}%
464: \figname{\kernelsm}
465: \topinsert
466: \vbox{
467: \centerline{\psfig{file=figures/kernelsm.p,width=3in}}
468: \medskip
469: \figurecaption{\kernelsm}{The Gaussian kernel smoother uses the Gaussian density function to assign weights to neighbouring points.
470: The points within the shaded neighbourhood are those that receive weights that are effectively nonzero for the computation of $\ev(Y\given x_0)$.}*
471: }
472: \endinsert
473: \par
474: The weight given to the $j$th point in producing the estimate at $x_0$
475: is defined by
476: $$S_{0j}={c_0\over \lambda}
477: d\left(\big\vert {{x_0-x_j}\over{\lambda}} \big\vert\right)\eqn{\kernel}$$
478: where $d(t)$ is an even function defined for $\abs{t}>0$.
479: The parameter $\lambda$ is the window width,
480: and the constant
481: $c_0$ is usually chosen so that the weights sum to unity, although
482: there are slight variations on this.
483: \index{standard gaussian density}%
484: A natural candidate for $d$ is the standard Gaussian density:
485: this gives the so-called { Gaussian kernel} smoother.
486: Other popular kernels, with some theoretical justification (^{Eubank, 1988}), are the Epanechnikov kernel
487: $$d(t)=\cases{ {3\over 4} (1-t^2),& for $\abs{t}\leq 1$;\cr 0& otherwise\cr}$$
488: which minimizes (asymptotic) mean squared error, and the minimum
489: variance kernel
490: $$d(t)=\cases{ {3\over 8} (3-5t^2),& for $\abs{t}\leq 1$;\cr 0& otherwise\cr}$$
491: which minimizes the asymptotic variance of the estimate.
492: \index{ metric distance}%
493: \index{metric}%
494: \index{rank distance}%
495: Note that the weight given to an observation is a
496: \index{kernel smoother}%
497: function only of its {\sl metric} distance from $x_0$, while the weights used by the nearest-neighbour smoothers are typically a function of both {\sl metric} and {\sl rank} distance.
498:
499: Kernel smoothers also exhibit biased endpoint behaviour.
500: Special kernels have been developed to overcome this bias; a simple approach is to use kernel weights in a locally-weighted straight-line fit.
501:
502: \index{kernel smooth}%
503: A Gaussian-kernel smooth for the diabetes data is shown in
504: \index{fitted smooth}%
505: Fig.~\allsmooths, where once again
506: the parameter $\lambda$ was chosen so that the approximate degrees of freedom of the
507: fitted smooth was five, the same as that of the other smooths in
508: Fig.~\allsmooths.
509:
510: \index{computational issues}%
511: \index{kernel smooth}%
512: \index{weight function}%
513: \sectionskip
514: \section{Computational issues}
515: One can visualize the action of the kernel smooth as sliding the weight function along the $x$-axis in short steps, each time computing the weighted mean of $y$.
516: The smooth is thus similar to a convolution between the kernel and an empirical step function defined on the data.
517: This is indeed the case, although the practical details obscure the resemblance.
518: For example, the kernel is usually truncated at the ends of the data, unless the data itself is periodic.
519: \index{kernel smooth}%
520: Typically the kernel smooth is computed as
521: $$\gsmooth(x_0)={\sum_{i=1}^n
522: d({x_0-x_i \over \lambda})y_i\over \sum_{i=1}^n d({x_0-x_i \over \lambda})},
523: \eqn{\NadWat}$$
524: and so both the numerator and denominator are convolutions.
525: If the $x$-values are evenly spaced (and preferably a power of two in number), great savings can be made by using the FFT (fast Fourier transform) in performing the calculations.
526: \index{spaced data}%
527: \index{fine grid}%
528: Although evenly spaced data are rarely the case in general regression contexts, reasonable approximations can often reduce the data to a fine grid, where the FFT can be put to work.
529: The details are laid out in Exercise~2.9.
530:
531: \index{running medians and enhancements}%
532: \index{smoothing}%
533: \index{running mean}%
534: \Sectionskip\Section{Running medians and enhancements}
535: A somewhat different approach to smoothing is based on
536: improvements to the simple
537: running mean smoother.
538: \index{running mean}%
539: First,
540: the running mean is replaced by a running median to make the smoother
541: resistant to outliers in the data.
542: \index{ hanning}%
543: \index{twicing}%
544: Then the appearance of the estimate is enhanced by applying compound operations
545: known as {\em Hanning}, {em splitting} and {\em twicing}, in various combinations.
546: \index{spaced data}%
547: \index{time series}%
548: \index{smoothing}%
549: We will not discuss these smoothers here because they are mainly
550: useful for evenly spaced data (sequences), especially time series, and they don't
551: typically provide an adequate amount of smoothing
552: for our purposes.
553: \index{interested reader}%
554: \index{bibliographic notes}%
555: The interested reader can pursue the references given in the Bibliographic notes.
556: \index{spaced data}%
557: \index{regression smoothers}%
558: We focus on smoothers for unevenly spaced data (scatterplots), sometimes
559: called {\em regression smoothers}.
560: In our view,
561: another disadvantage of the enhanced running median smoothers is that they are
562: highly nonlinear functions of the response, and thus it is difficult to assess
563: the amount of fitting that they do.
564: This is an important consideration for the inferential stage of a data analysis.
565: \index{ linear}%
566: \index{ degrees of freedom}%
567: Many of the smoothers that we discuss are {\em linear} (section~3.4.2), and this facilitates
568: an approximate assessment of their {\em degrees of freedom} (section~3.5).
569:
570: \index{time series}%
571: We return to the time series setting for our discussion of the Fourier analysis
572: of smoothers in section~3.7, and a seasonal decomposition procedure in
573: section~8.5.
574: \index{equivalent kernel}%
575: \figname{\hatplot}
576: \pageinsert
577: \vbox to \vsize{
578: \centerline{\psfig{file=figures/hatplot.p,width=\hsize}}
579: \medskip
580: \index{target point}%
581: \index{equivalent kernel}%
582: \figurecaption{\hatplot}{The equivalent kernels for a variety of smoothers, based on the data used in Fig.~\allsmooths. The arrows indicate the target point in each case. The points are plotted at the pairs $\{x_j,\,S_{0j}\}$ }*
583: \vfill
584: }
585: \endinsert
586: \Sectionskip
587: \Section{Equivalent kernels}
588: The smoothers described in this chapter are defined in quite different ways, with
589: \index{ equivalent kernels}%
590: more different ways to come.
591: Their {\sl equivalent kernels} are one way to compare
592: \index{ linear}%
593: \index{smoothing parameter}%
594: them on common ground.
595: All the smoothers studied in this chapter are {\sl linear} in $Y$, which means that the fit at a point $x_0$ can be written as $\gsmooth(x_0)=\sum_{j=1}^n S_{0j}y_j$, and the $S_{0j}$ depend on all the $x_i$ and on the smoothing parameter $\lambda$.
596: \index{ equivalent kernel}%
597: Thus the weight in the fit at $x_0$ which is associated with the point $x_j$ is $S_{0j}$, and this sequence of weights is known as the {\em equivalent kernel} at $x_0$.
598: A plot of $S_{0j}$ against $x_j$ shows which observations have an influence on the fit at $x_0$.
599: \index{running mean}%
600: Looking at Fig.~\hatplot, we see that the running mean smoother has weights
601: either zero or $1/m$, where $m$ is the number of observations in the neighbourhood.
602: We have used the same diabetes data as was used in Fig.~\allsmooths.
603:
604: \index{equivalent kernel}%
605: For each smoother we have computed the equivalent kernel centered at two points: one near the boundary and one at the centre.
606: The weights at the boundary are larger and warn us that end effects may be a problem (as indeed they are).
607:
608: \index{running mean}%
609: For the running-line smoother the weights change linearly within the window (Exercise~2.3).
610: For both the running mean and line smoothers the weights drop off
611: abruptly to zero outside the neighbourhood, and account for their jagged appearance.
612: \index{ loess smooth}%
613: The {\em loess} smooth on the other hand has a strictly local neighbourhood yet the weights die down smoothly to zero.
614:
615: \index{smoothing}%
616: \index{smoothing parameter}%
617: \index{equivalent kernel}%
618: \index{ equivalent degrees of freedom}%
619: We need to calibrate the smoothers so that they are doing approximately the same amount of smoothing, since the value of the smoothing parameter will clearly widen or narrow the equivalent kernels.
620: We do this using the {\sl equivalent degrees of freedom}, which we describe in the next chapter.
621: \index{smoothing}%
622: \index{equivalent kernel}%
623: All the smoothers in Fig.~\hatplot\ have been calibrated to do about five degrees of freedom worth of smoothing.
624:
625: In the top panel we see the equivalent kernel for the quartic polynomial fit.
626: As we might expect it spreads its influence everywhere.
627:
628: We will return to Fig.~\hatplot\ as we encounter the remaining smoothers.
629:
630: \index{regression splines}%
631: \Sectionskip
632: \Section {Regression splines}
633: Polynomial regression has limited appeal due to the global nature of the
634: fit, while in contrast the smoothers we have seen so far have an explicit local nature.
635: \index{regression splines}%
636: \index{piecewise}%
637: Regression splines offer a compromise by representing the fit as a {\sl piecewise}
638: \index{ knots}%
639: polynomial.
640: The regions that define the pieces are separated by a sequence
641: \index{piecewise polynomials}%
642: of {\sl knots} or breakpoints, $\xi_1,\ldots,\xi_K$.
643: \index{bin smoother}%
644: \index{piecewise cubic polynomials}%
645: In addition, it is customary to enforce the piecewise polynomials to join smoothly at these knots.
646: Although many different configurations are possible (the bin smoother was one), a popular choice
647: consists of piecewise cubic polynomials constrained to be continuous and have continuous first and second derivatives at the knots.
648: Apparently our eyes are skilled at picking up second order and lower discontinuities, but not higher.
649: \figname{\BS}
650: \midinsert
651: \vbox{
652: \centerline{\psfig{file=figures/BS.p,width=\hsize}}
653: \figurecaption{\BS}{Piecewise cubic fits to some simulated data. The order of continuity increases from 0 --- discontinuous as in figure (a), to 3 --- continuous second derivatives as in figure (d).
654: The fitted values are connected within each region, leaving gaps at the knots to make the orders of continuity more visible.
655: The vertical lines indicate the location of the two knots.
656: }*
657: }
658: \endinsert
659: \par
660: \index{interior knots}%
661: To illustrate this point, Fig.~\BS\ shows four piecewise cubic fits to some simulated data, with 2
662: interior knots and different orders of continuity.
663: Apart from being smoother, the fit in panel (d) required fewer parameters (6 altogether) than the others.
664:
665: By allowing more knots, the family of curves becomes more flexible.
666: \index{multiple regression}%
667: \index{basis functions}%
668: \index{piecewise cubic polynomials}%
669: For any given set of knots, the smooth is computed by multiple regression on an appropriate set of basis vectors.
670: These vectors are the basis functions representing the particular family of piecewise cubic polynomials evaluated at the observed values of the predictor $X$.
671: This approach is attractive because
672: of its computational and statistical simplicity; for example standard
673: parametric inferential methods can be used to test the importance
674: of any of the parameters.
675:
676: \index{natural splines}%
677: \index{piecewise polynomials}%
678: \index{cubic splines}%
679: \index{boundary knots}%
680: A variant of polynomial splines are the natural splines; although they are defined for all piecewise polynomials of odd degree, we discuss the natural {\sl cubic} splines.
681: These are cubic splines with the additional constraint that the function is linear beyond the boundary knots.
682: To enforce this condition we have to impose the two constraints in each of the boundary regions: $f'''=f''=0$, which reduces the dimension of the space from $K+4$ to $K$ if there are $K$ knots.
683: \index{interior knots}%
684: \index{boundary knots}%
685: In practice it is common to supply an additional knot at each extreme of the data, and impose the linearity beyond them.
686: \index{natural splines}%
687: \index{regression splines}%
688: Then with $K$ interior knots (and two boundary knots), the dimension of the space of fits is $K+2$.
689: Natural splines have less flexibility at the boundaries, but this tends to be a plus since the fitted values of regular regression splines have high variance near the boundary.
690: \index{interior knots}%
691: \index{computational aspects}%
692: Furthermore, for the same number of parameters, we get two more interior knots.
693:
694: \sectionskip
695: \section{Computational aspects}
696: The main difficulty when working with regression splines is
697: \index{ number}%
698: \index{interior knots}%
699: selecting the number and position of the knots.
700: A very simple approach (referred to as cardinal splines) requires a single parameter, the {\sl number} of interior knots.
701: The positions are then chosen uniformly over the range of the data.
702: \index{three interior knots}%
703: A slightly more adaptive version will place the knots at appropriate quantiles of the predictor variable; e.g., three interior knots would be placed at the three quartiles.
704: More adventurous schemes use data driven criteria to select the number and positions of the knots.
705: The main challenge is to come up with a sensible procedure, while avoiding a combinatorial nightmare.
706: We describe some specific proposals in Chapter~9.
707:
708: \index{basis functions}%
709:
710: \index{interior knots}%
711: \index{boundary knots}%
712: Another computational aspect concerns the choice of basis functions for representing the splines (a vector space) for a given set of knots.
713: \index{basis functions}%
714: Suppose these interior knots are denoted by $\xi_1<\cdots<\xi_K$, and for notational simplicity we augment the set with two boundary knots $\xi_0$ and $\xi_{K+1}$.
715: A simple choice of basis functions for piecewise-cubic splines, known as the truncated power series basis, derives from the parametric expression for the smooth
716: $$\gsmooth(x)= \beta_0+\beta_1x+\beta_2x^2+\beta_3x^3+\sum_{j=1}^K\theta_j(x-\xi_j)_+^3.\eqn{\tpseries}$$ where the $a_+$ denotes the positive part of $a$.
717: \index{third derivative}%
718: Evidently (Exercise~2.4) \tpseries\ has the required properties:\smallskip
719: {\parindent 20pt
720: \item{(i)} $\gsmooth$ is cubic polynomial in any subinterval $[\xi_j,\xi_{j+1})$,
721: \item{(ii)} $\gsmooth$ has two continuous derivatives, and
722: \item{(iii)} $\gsmooth$ has a third derivative that is a step function with jumps at $\xi_1,\ldots,\xi_K$.
723:
724: \index{basis functions}%
725: }\smallskip
726:
727: We can write \tpseries\ as a linear combination of $K+4$ basis functions $P_j(x)$: $P_1(x)=1$, $P_2(x)=x$, and so on.
728: Each of these functions must also satisfy the three conditions and be linearly independent in order to qualify as a basis.
729: \index{piecewise cubics}%
730: Even without \tpseries\ it is clear that we require $K+4$ parameters to represent piecewise cubics; four for each of $(K+1)$ cubics, less three per interior knot due to the constraints.
731: \index{matrix}%
732: To actually smooth some data pairs $\{x_i,y_i\}$, we would construct a regression matrix with $K+4$ columns, each column corresponding to a function $P_j$ evaluated at the $n$ values of $x$.
733:
734:
735: Although \tpseries\ has algebraic appeal, it is not the recommended form for computing the regression spline.
736: Even though a particular cubic piece is designed to accommodate one interval, it is evaluated at all points to the right of its knot and the numbers usually get large.
737: \index{basis functions}%
738: \figname{\BSbasis}
739: \midinsert
740: \vbox{
741: \centerline{\psfig{file=figures/BSbasis.p,width=\hsize}}
742: \index{interior knots}%
743: \figurecaption{\BSbasis}{The six $B$-spline basis functions for representing the fits in Fig.~\BS(d). The two interior knots are indicated by the vertical broken lines.
744: The points indicate the occurrence of data, and are $\{x_i,B_j(x_i)\}$ for $j=1,\ldots,6$.
745: \index{basis functions}%
746: }*
747: }
748: \endinsert
749: \par
750: The $B$-spline basis functions provide a numerically superior alternative basis to the truncated power series.
751: \index{matrix}%
752: Their main feature is that any given basis function $B_j(x)$ is nonzero over a span of at most five distinct knots.
753: In practice this means that their evaluation rarely gets out of hand, and the resulting regression matrix is banded.
754: \index{piecewise cubics}%
755: Of course the $B_j$ are themselves piecewise cubics, and we need $K+4$ of them if we want to span the space.
756: Their definition is rather simple (in terms of divided differences), but unhelpful for fully understanding their properties, so we omit it here and refer the reader to ^{de~Boor~(1978)}.
757: Figure~\BSbasis\ shows the evaluated $B$-splines used to compute the fits in Fig.~\BS(d).
758:
759: \index{regression splines}%
760: \index{when the knots are given}%
761: In summary, regression splines are attractive because of their computational neatness, {\sl when the knots are given}.
762: \index{regression splines}%
763: Standard linear model estimation can be applied, which will be found to be very convenient in later chapters when we use regression splines in additive models.
764:
765: However the difficulty of choosing the number and position of the knots
766: tends to be a drawback of this approach.
767: \index{equivalent kernel}%
768: \index{three interior knots}%
769: When a small number of knots is used, the smoother can show some disturbing nonlocal behaviour.
770: Figure~\hatplot\ shows the equivalent kernels for a natural cubic regression spline with three interior knots.
771: The right-hand kernel does not look much different from that for the quartic polynomial! With more knots this global influence would be dampened, but we often don't have that many \df\ to spare.
772:
773: \index{regression splines}%
774: \index{smoothing parameter}%
775: Another problem with regression splines is that the smoothness of the
776: estimate can not easily be varied continuously as a function of a single smoothing
777: parameter; this is an attractive, if only approximate, property of most of the
778: \index{cubic smoothing spline}%
779: other smoothers that we discuss.
780:
781: \Sectionskip
782: \Section{Cubic smoothing splines}
783: This smoother is not constructed explicitly like those described so far, but instead emerges as the solution to an optimization
784: problem.
785:
786:
787: Consider the following problem: among all functions
788: $f(x)$ with absolutely continuous first derivatives and integrable second derivatives, find the one that minimizes
789: $$\sum_{i=1}^n\{y_i-f(x_i)\}^2 +\lambda\int_{-\infty}^\infty \{f''(t)\}^2 \, dt\eqn{\qb}$$ where
790: $\lambda$ is a fixed constant.
791: \index{scatterplot smoother}%
792: This criterion satisfies the requirements for a scatterplot smoother that we
793: mentioned earlier.
794: The first term measures closeness to the data while the
795: second term penalizes curvature in the function.
796: Now it is one thing to state
797: a criterion like $\qb$ and quite another to find the optimizing function.
798: \index{ natural cubic spline}%
799: Remarkably, it can be shown that $\qb$ has an explicit, unique minimizer and
800: that minimizer is a {\sl natural cubic spline} with knots at the unique values of $x_i$
801: (Exercise~2.6).
802: \index{interior knots}%
803: At face value it seems that the family is overparametrized, since there are as many as $n-2$ interior knots.
804: \index{ effective dimension}%
805: This would result in $n+2$ parameters, although the constraints on each end bring it down to $n$. We'll see however that the coefficients are estimated in a constrained way as well, and this can bring the {\sl effective} dimension down dramatically.
806:
807:
808: \index{smoothing}%
809: A cubic
810: \index{running-lines smooth}%
811: smoothing spline fitted to the diabetes data is shown in Fig.~\allsmooths.
812: It looks
813: much like the
814: running-lines smooth
815: but is less jagged.
816:
817:
818: \index{running-lines smooth}%
819: The parameter $\lambda$ plays the same role as the span in the running-lines smooth.
820: Large values of $\lambda$ produce smoother curves while
821: \index{least-squares line}%
822: smaller values produce more wiggly curves. At the one extreme, as
823: $\lambda\rightarrow \infinity$, the penalty term dominates, forcing
824: $f''(x)=0$, and thus the solution is the least-squares line.
825: At the other extreme, as $\lambda\rightarrow 0$, the penalty term becomes unimportant and the
826: solution tends to an interpolating twice-differentiable function.
827: We discuss methods for choosing $\lambda$ in Chapter~3.
828: \index{computational aspects}%
829: \index{natural cubic spline}%
830: \index{interior knots}%
831:
832: \sectionskip
833: \section{Computational aspects}
834: Using the fact that the solution to $\qb$ is a natural cubic spline with $n-2$ interior knots, we can represent it in terms of a basis for this space of fits.
835: \index{basis functions}%
836: For computational convenience, we will use the unconstrained $B$-spline basis, and write
837: $\gsmooth(x)=\sum_1^{n+2} \gamma_i B_i(x)$, where $\gamma_i$ are
838: \index{natural splines}%
839: coefficients and the $B_i$ are the cubic $B$-spline basis functions.
840: As written here, $s$ lies in an $(n+2)$-dimensional space; however, the natural splines are a subspace.
841: We can now simply replace $f$ in \qb\ by $s$ and perform the integration.
842: \index{matrix}%
843: \index{matrix}%
844: Then defining the $n\times (n+2)$ matrix $ \bB$ and $(n+2)\times (n+2)$ matrix $\fat{\Omega}$ by
845: $$
846: B_{ij}=B_j(x_i)$$
847: and
848: $$\fat{\Omega}_{ij}=\int_{-\infty}^\infty B_i''(x)B_j''(x) \, dx,$$
849: we can rewrite the criterion $\qb$ as
850: $$(\vec y-\bB\fat\gamma)^T(\vec y-\bB\fat\gamma) +\lambda\fat \gamma^T\fat{\Omega}
851: \fat \gamma\eqn{\qd}$$
852: \index{ sorted values}%
853: Although at face value it seems that there are no boundary derivative constraints, it turns out that the penalty term automatically imposes them.
854:
855: Setting the derivative with respect to $\gamma$ equal to zero gives
856: $$(\bB^T\bB+\lambda\fat\Omega)\hatfat\gamma=\bB^T\vec y\eqn{\eqq}$$
857: \index{matrix}%
858: Since the columns of $\bB$ are the evaluated $B$-splines, in order from left to right and evaluated at the {\sl sorted} values of $X$, and the cubic $B$-splines have local support, $\bB$ is lower 4-banded.
859: Consequently the matrix $\vec M=(\bB^T\bB+\lambda\fat{\Omega})$ is
860: 4-banded and hence its Cholesky decomposition $\vec M=\bL\bL^T$ can be computed easily.
861: One
862: then solves $\bL\bL^T\hatfat\gamma=\bB^T\vec y$ by back-substitution to give $\hatfat\gamma$
863: and hence the solution $\hat\gsmooth$ in O(n) operations.
864:
865:
866: For theoretical purposes, it is convenient to rewrite the solution
867: vector $\vgsmooth$ ($\gsmooth$ evaluated at each of the $n$ values of $X$ in the sample) in another form.
868: \index{ natural-spline basis}%
869: \index{matrix}%
870: Let $\bN$ be an $n\times n$ nonsingular {\sl natural-spline} basis matrix for representing the solution (Exercise~2.5).
871: Denote by $\hatfat{\beta}$ the transformed version of $\hatfat\gamma$ corresponding to the change in basis. Then we can write
872: $$\vgsmooth=\bN\hatfat{\beta}=\bN(\bN^T\bN+\lambda\fat{\Omega})^{-1}\bN^T\vec y=
873: (\bI+\lambda \bK)^{-1}\vec y\eqn{\qf}$$
874: where $\bK={\bN^{-T}}\fat{\Omega} \bN^{-1}$.
875: \index{cubic smoothing spline}%
876: In terms of the candidate fitted vector $\vec f$ and $\bK$, the cubic smoothing spline $\vgsmooth $ minimizes
877: $$(\vec y-\vec f)^T(\vec y-\vec f) +\lambda\vec f^T\bK
878: \vec f\eqn{\nqd}$$ over all vectors $\vec f$.
879: It is appropriate to call the term $\vec f^T\bK\vec f$ a roughness penalty, since it can be shown to be a quadratic form in second differences.
880: \index{cubic smoothing spline}%
881: \index{local averaging}%
882: The cubic smoothing spline, because of the implicit way it is defined, doesn't appear
883: to use {\sl local averaging}.
884: \index{equivalent kernel}%
885: \index{target point}%
886: However, the equivalent-kernel weights in Fig.~\hatplot\ show that it does possess local behaviour quite similar to kernels or locally-weighted lines.
887: \index{equivalent kernel}%
888: \index{equivalent kernel}%
889: \index{smoothing splines}%
890: The equivalent kernel is nowhere nonzero, but is close to zero far from the target point.
891: Other examples of equivalent kernels are given in ^{Buja, Hastie and Tibshirani (1989)}, and
892: an asymptotic form for the equivalent kernel of smoothing splines is derived by ^{Silverman (1984)}.
893: \index{smoothing}%
894: \index{kernel}%
895: One might say, then, that a cubic smoothing spline is approximately a {\sl kernel}
896: \index{smoothing parameter}%
897: \index{weight function}%
898: smoother.
899:
900: The smoothing parameter $\lambda$ controls the shape of the kernel
901: or weight function.
902: \index{target point}%
903: As one would expect,
904: when $\lambda$ is decreased
905: the weights tend to concentrate more around the target point, while as
906: \index{smoothing splines}%
907: $\lambda$ is increased, they spread out more.
908: We discuss the properties of smoothing splines in more detail in Chapters~5 and section~9.3.6.
909:
910: \index{locally weighted running-line smoothers}%
911: \index{loess}%
912: \Sectionskip
913: \Section{Locally weighted running-line smoothers}
914: Here we describe in more detail the locally-weighted smoother of ^{Cleveland (1979)}, currently called {\sl loess} in the S statistical-computing language.
915: Although any polynomial can be fitted locally, we discuss local lines.
916:
917: \index{nearest neighbours}%
918: \index{nearest neighbours}%
919: A locally-weighted straight-line smooth $\gsmooth(x_0)$ using $k$ nearest neighbours is computed in a number of steps:\smallskip
920: {\parindent 20pt
921: \item{(i)} The $k$ nearest neighbours of $x_0$ are identified, denoted by
922: $\NN(x_0)$.
923: \item{(ii)}$\Delta(x_0)=\max_{N(x_0)}\abs{x_0-x_i}$ is computed, the distance of the furthest near-neighbour from $x_0$.
924: \index{ tri-cube weight}%
925: \item{(iii)}Weights $w_i$ are assigned to each point in $\NN(x_0)$, using the {\sl tri-cube} weight function:
926: $$W\biggl({\abs{x_0-x_i}\over {\Delta(x_0)}}\biggr)
927: \eqn{\clev}$$
928: where $$W(u)=(1-u^3)^3 \quad\for\quad 0\leq u < 1 \eqn{\tri}$$
929: and zero otherwise.
930: \index{weighted least-squares fit}%
931: \item{(iv)} $\gsmooth(x_0)$ is the fitted value at $x_0$ from the weighted least-squares fit of $y$ to $x$ confined to $\NN(x_0)$ using the weights computed in (iii).
932:
933: \index{target point}%
934: }\smallskip
935: \figname{\loessm}
936: \topinsert
937: \vbox{
938: \centerline{\psfig{file=figures/loessm.p,width=3in}}
939: \medskip
940: \index{target point}%
941: \figurecaption{\loessm}{The loess smoother first computes the $\lambda$~\% nearest-neighbours to the target point, indicated here by the shaded region.
942: \index{target point}%
943: A tricube kernel, centered at the target point, becomes zero exactly at the furthest neighbour. The smooth at the target point is the fitted value from the locally-weighted linear fit.}*
944: }
945: \endinsert
946: \par
947: Figure~\loessm\ depicts the situation, using the diabetes data once again.
948: Notice that the kernel is truncated at one end due to the metric asymmetry of the neighbourhood.
949:
950: Cleveland also discusses the use of a robust regression within each
951: \index{smoothing}%
952: neighbourhood, to protect against outliers.
953: This effectively works by repeatedly smoothing the data, and at each iteration down-weighting points with large residuals.
954: See Chapters~3 and 9 for more details.
955:
956: A locally-weighted running-line smooth
957: \index{ span}%
958: \index{data points}%
959: \index{smoothing parameter}%
960: is shown in Fig.~\allsmooths.
961: The number of nearest-neighbours, usually expressed
962: as a percentage or {\sl span} of the data points, is the smoothing parameter.
963: In Figs~\allsmooths\ and \loessm\ a span of $21/43$ was used, chosen to make
964: the degrees of freedom approximately equal to that of the other smooths.
965:
966: \index{nearest neighbourhood}%
967: \index{symmetric nearest neighbourhood}%
968: \index{target point}%
969: \index{nearest neighbourhood}%
970: In principle, nearest neighbourhoods
971: are preferable to
972: symmetric nearest neighbourhoods because in a neighbourhood with
973: a fixed number of points, the average distance of the points to the target point
974: is less in the nearest neighbourhood (unless the predictors are evenly spaced).
975: In general this should result in less bias.
976: \index{symmetric nearest neighbourhood}%
977: \index{running-lines smoother}%
978: We use symmetric nearest neighbourhoods in a running-lines smoother, however,
979: because of their superior performance at the left and right endpoints.
980: \index{nearest neighbourhood}%
981: \index{data points}%
982: \index{target point}%
983: A nearest neighbourhood at the endpoint contains the same number of data
984: points
985: as a neighbourhood in the middle of the data,
986: so that when a least-squares fit is applied there, too much
987: weight is given to points far away from the target point.
988: \index{symmetric nearest neighbourhood}%
989: A symmetric nearest neighbourhood, on the other hand, contains only about
990: half of the number of points as a full neighbourhood, and thus
991: automatically reduces the weight given to far way points.
992: \index{nearest neighbourhood}%
993: \index{loess}%
994: Nearest neighbourhoods work satisfactorily with {\sl loess} at the endpoints, however,
995: because the tri-cube function does the job of down-weighting far away points.
996:
997: \index{kernel smoother}%
998: The locally-weighted smoothers are popular, since they enjoy the best of both worlds.
999: They share the ability of near-neighbour smoothers to adapt their bandwidth
1000: to the local density of the predictors, while they have the smoothness features
1001: of kernel smoothers.
1002:
1003: Although in principal it requires $O(n^2)$ operations to compute a locally-weighted smooth, the current implementations compute the fit over a grid of values of $x$, and use interpolation elsewhere.
1004:
1005: \index{smoothers for multiple predictors}%
1006: \index{single predictor}%
1007: \index{scatterplot smoother}%
1008: \Sectionskip
1009:
1010: \Section{Smoothers for multiple predictors}
1011: So far we have talked about smoothers for a single predictor, that is,
1012: scatterplot smoothers.
1013: What if we have more than one predictor, say $X_1,\ldots, X_p?$.
1014: Then our problem is one of fitting a $p$-dimensional surface to $Y$.
1015: \index{multiple regression}%
1016: The multiple regression of $Y$ on $X_1,\ldots, X_p$ provides a
1017: simple, but very limited, estimate of the surface.
1018: \index{running mean}%
1019: \index{kernel smoother}%
1020: On the other hand,
1021: it is easy conceptually to generalize the running mean, locally-weighted
1022: running-lines, and kernel smoothers to this setting.
1023: \index{ nearest neighbours}%
1024: The first two smoothers require a definition of {\sl nearest neighbours}
1025: of a point in $p$-space.
1026: \index{ nearest}%
1027: \index{euclidean distance}%
1028: {\sl Nearest} is determined by a distance measure and for this
1029: the most obvious choice is Euclidean distance.
1030: \index{symmetric nearest neighbours}%
1031: Note that the notion of symmetric nearest neighbours is no longer meaningful
1032: when $p > 1$.
1033: \index{running mean}%
1034: \index{target point}%
1035: Having defined a neighbourhood,
1036: the generalization of the running mean estimates
1037: the surface at the target point by averaging the
1038: response values in the neighbourhood.
1039:
1040: \index{euclidean distance}%
1041: \index{matrix}%
1042: This
1043: highlights an important detail: what shape neighbourhood do we use?
1044: Nearest-neighbourhoods defined in terms of Euclidean distance are typically
1045: spherical, but one could imagine more general
1046: neighbourhoods defined in terms of a covariance matrix $\Sigma$ of the predictors.
1047: \index{target point}%
1048: Such a metric would consider points lying on an ellipse centered at the target
1049: point to be equidistant from the target.
1050: This generalization is not a purely
1051: academic issue --- it may be important if the covariates are measured in different
1052: units or are correlated. In linear
1053: regression, the coefficients automatically scale each covariate
1054: \index{smoothing}%
1055: correctly.
1056: Various strategies are used in practice, an obvious one being to standardize the individual variables prior to smoothing.
1057: This might not always be the best choice, however.
1058: For example, if the underlying surface changes more rapidly with one variable
1059: than the other, one may wish to have the neighbourhoods thinner in the direction of the first variable.
1060: \index{running mean}%
1061: This can be achieved by differential scaling of the variables.
1062: \figname{\bivar}\par
1063: The kernel smoother is generalized in an analogous way
1064: to the running mean.
1065: \index{predictor space}%
1066: Given two $p$-vectors
1067: $\vec x^0$ and $\vec x^i$
1068: in predictor space,
1069: the weight given to this $i$th point for the fit at $\vec x^0$
1070: is $S_{0i}=(c_0/\lambda)
1071: d({{\norm{\vec x^0-\vec x^i}}/\lambda})$
1072: where $d(t)$ is an inverse distance measure
1073: and $\norm{\cdot}$ is a norm, for example squared distance.
1074: The constant $c_0$ is usually chosen to make the weights sum to unity, so that the smoother
1075: reproduces the constant function.
1076: The multi-predictor Gaussian kernel uses
1077: $d(t)$ equal to the standard Gaussian density.
1078: \index{weight function}%
1079: The resulting weight function is depicted in Fig.~\bivar.
1080: \index{weight function}%
1081: \index{kernel smoother}%
1082: \midinsert
1083: \vbox{
1084: \centerline{\psfig{file=figures/bivar.p,width=3in}}
1085: \figurecaption{\bivar}{ The weight function for a bivariate Gaussian kernel smoother.
1086: \index{predictor space}%
1087: }*
1088: }
1089: \endinsert
1090: The choice of norm is important in this setting, and determines the shape of
1091: the neighbourhoods in predictor space.
1092: \index{matrix}%
1093: One might argue, for example, that the norm should incorporate
1094: covariance information by replacing the
1095: squared error norm with the ``Mahalanobis'' distance
1096: $$\norm{\vec x^0-\vec x^i}=(\vec x^0-\vec x^i)^T\Sigma^{-1}(\vec x^0-\vec x^i)$$
1097: where
1098: $\Sigma$ is the covariance matrix of the predictors.
1099:
1100:
1101: Locally-weighted lines generalize equally easily to two or more dimensions.
1102: The simplest generalization hardly needs any different wording than that used to describe the univariate situation.
1103: Once a multivariate distance is defined,
1104: nearest-neighbourhoods and tri-cube weights get assigned in exactly the same way.
1105: One computes a local plane rather than a line.
1106:
1107: \index{smoothing}%
1108: \index{thin-plate spline}%
1109:
1110: The cubic smoothing spline is more difficult to generalize to two or higher
1111: dimensions: the so-called {\sl thin-plate spline} is one such
1112: generalization.
1113: It derives from generalizing the second derivative penalty for smoothness to a two dimensional Laplacian penalty of the form
1114: $$\int\int\left\{\left({\partial^2 f\over \partial x_1^2}\right)^2 + \left({\partial^2 f\over \partial x_1\partial x_2}\right)^2+\left({\partial^2 f\over \partial x_2^2}\right)^2\right\}\, dx_1dx_2.$$
1115: The computations are $O(n^3)$ in number compared to the $O(n)$ for univariate splines.
1116: \index{interested reader}%
1117: \index{bibliographic notes}%
1118: Details are not given here; the interested reader can find references in the
1119: bibliographic notes at the end of Chapter~3.
1120:
1121: \index{ tensor product}%
1122: \index{regression splines}%
1123: Another generalization is known as multivariate {\sl tensor product} splines.
1124: \index{basis functions}%
1125: \index{basis functions}%
1126: These are also useful for generalizing univariate regression splines.
1127: \index{regression splines}%
1128: The basic idea is to construct two-dimensional basis functions by multiplying together one-dimensional basis functions in the respective predictors.
1129: \index{smoothing}%
1130: We discuss these in more detail for multivariate regression splines in Chapter~9.
1131: \index{basis functions}%
1132: For the smoothing version of tensor product splines, the effective basis set has dimension $n^2$ formed by the tensor product of the univariate bases.
1133: The fit will be a damped regression onto this space of basis functions.
1134: \index{multi-predictor smoothers}%
1135: See the references at the end of Chapter~3 for details.
1136:
1137: We are not devoting much space to multi-predictor smoothers because we don't
1138: feel that they are very useful for more than two or three predictors.
1139: Indeed, their many shortcomings (e.g.~difficulty of interpretation and computation)
1140: provide an impetus for studying additive models, the central
1141: topic of this book.
1142: \index{parametric regression}%
1143: These comments refer to the generic multivariate smoothers as described here.
1144: \index{smoothing}%
1145: We are not referring to some of the adaptive multivariate nonparametric regression methods, which might also be termed surface-smoothers, and which were designed to overcome some of these objectional aspects.
1146: We discuss some of the problems associated with multi-predictor smoothing in detail at the
1147: beginning of Chapter~4.
1148: \index{multi-predictor smoothers}%
1149: Despite these criticisms, multi-predictor smoothers (especially the bi-predictor) are useful
1150: \index{scatterplot smoother}%
1151: \index{running-lines smoother}%
1152: \index{smoothing}%
1153: \index{locally-weighted running-lines}%
1154: \index{kernel smoother}%
1155: in certain settings, as later examples show.
1156:
1157: \Sectionskip
1158: \Section{Discussion and bibliography}
1159: In this chapter we have described a number of scatterplot smoothers,
1160: including the running-mean and running-lines smoothers, cubic regression and smoothing
1161: splines, locally-weighted running-lines and kernel smoothers.
1162: We have also touched on surface smoothers.
1163: Our discussion has focused on smoothers with which we have some experience, and any imbalance in our introduction simply reflects this.
1164: \index{ per~se}%
1165: \index{building block}%
1166: In addition, the emphasis in the book is not on smoothers {\sl per~se}, but
1167: on
1168: smoothers as building blocks in the estimation of additive models.
1169:
1170: \index{smoothing parameter}%
1171: \index{locally-weighted running-lines}%
1172: \index{smoothing}%
1173: \index{kernel smoother}%
1174:
1175: Some recent theoretical results suggest that for appropriately chosen
1176: smoothing parameters, there are not likely to be large differences between
1177: locally-weighted running-lines, cubic smoothing-splines and kernel smoothers
1178: (^{Silverman, 1984}; ^{Muller, 1987}).
1179: In order to assess the finite-sample operating characteristics of these
1180: and other smoothers, a comprehensive Monte Carlo study, akin to the
1181: Princeton robustness study, would be very useful.
1182: A recent report by ^{Breiman}
1183: and Peters (1988) has this flavour.
1184:
1185: For a detailed bibliography, see the notes at the end of Chapter~3.
1186: Here we list some introductory references we have found useful.
1187: \index{smoothing splines}%
1188: ^^{Eubank, R.L.} Eubank (1988) covers most of the material in various levels of detail, especially smoothing splines.
1189: \index{smoothing splines}%
1190: ^^{Silverman, B.W.} Silverman's (1985) paper on smoothing splines not only gives an easy to read introduction but also touches on a variety of applications.
1191: ^^{Cleveland, W. S.} Cleveland (1979) introduces the locally-weighted running-line smoother for univariate problems, and
1192: ^^{Devlin, S.}^^{Grosse, E. H.} ^{Cleveland, Devlin and Grosse (1988)} describe the multivariate version.
1193:
1194: \Sectionskip
1195: \Section{Further results and exercises 2}
1196: \Mark{EXERCISES \ 2}%
1197: \beginexercises%
1198: \exercise For the diabetes data, compute one of each of the smooths described
1199: in this Chapter
1200: for the variable \name{base deficit}.
1201: Compare qualitatively the resulting estimates in terms of notable features
1202: and smoothness.
1203:
1204: \exercise {\sl Updating formula for running-line smooth.}
1205: Suppose we add a point $(x_{j+1},y_{j+1})$ to a neighbourhood
1206: containing $j$ points.
1207: If the means, variance of $X$, and covariance of $X$ and $Y$ for the first $j$ points are
1208: $\bar x_j$, $\bar y_j$, $S^x_j$, and $S^{xy}_j$, show that
1209: $$\bar x_{j+1}=(j\bar x_j+x_{j+1})/(j+1)$$
1210: $$\bar y_{j+1}=(j\bar y_j+y_{j+1})/(j+1)$$
1211: $$(j+1)S^x_{j+1}=jS^x_j+{{j+1}\over j}(x_{j+1}-\bar x_{j+1})^2$$
1212: $$(j+1)S^{xy}_{j+1}=jS^{xy}_j+{{j+1}\over j}(x_{j+1}-\bar x_{j+1})(y_{j+1}-
1213: \bar y_{j+1}).\eqn{\upd}$$
1214:
1215: What are the equivalent equations for deleting a point?
1216: Together these
1217: can be used to update the least-squares slope and intercept, and hence the entire fit can be computed in $O(n)$ operations.
1218:
1219: \exercise Derive an expression for $S_{ij}$, the coefficient of $y_j$ in the
1220: expression $\gsmooth(x_i)=\sum_{j=1}^n S_{ij}y_j$ for the running-line fit at $x_i$,
1221: using $\symNN(x_i)$ to denote the set of indices of the $k$ symmetric nearest-neighbours to $x_i$.
1222:
1223:
1224: \exercise Show that the truncated power series representation $$\gsmooth(x)= \beta_0+\beta_1x+\beta_2x^2+\beta_3x^3+\sum_{j=1}^k\theta_j(x-\xi_j)_+^3$$ satisfies the three conditions of a cubic spline given in the text.
1225:
1226: \index{natural splines}%
1227: \index{matrix containing}%
1228: \index{basis functions}%
1229: \index{interior knots}%
1230: \exercise {\sl Basis for natural splines.} Suppose $B$ is an $n\times (K+4)$ matrix containing the evaluations of the cubic $B$-spline basis functions with $K$ interior knots evaluated at the $n$ values
1231: \index{matrix containing}%
1232: \index{basis functions}%
1233: of $X$.
1234: Let $C$ be the $2\times (K+4)$ matrix containing the second derivatives
1235: \index{matrix}%
1236: \index{interior knots}%
1237: \index{boundary knots}%
1238: of the basis functions at the boundary points $x_1$ and $x_n$.
1239: Show how
1240: \index{ derivation of smoothing splines}%
1241: to derive $N$ from $B$, an $n\times (K+2)$ basis matrix for the natural cubic splines with the same interior knots and boundary knots at the extremes of the $X$.
1242:
1243: \exercise {\sl Derivation of smoothing splines; ^{Reinsch (1967)}.} Consider the following optimization problem: minimize
1244: $$\int_{-\infty}^\infty \{f''(x)\}^2 \, dx
1245: \quad \hbox{subject to}\quad \sum_{i=1}^n \{y_i-f(x_i)\}^2 \leq \sigma
1246: \eqn{\sopt}$$
1247: over all twice differentiable
1248: functions $f(x)$.
1249: To solve this problem, consider the Lagrangian functional
1250: $$F(f)=\int_{-\infty}^\infty\{f''(x)\}^2 \, dx + \rho \left[ \sum_{i=1}^n \{y_i-f(x_i)\}^2-\sigma+z^2\right]\eqn{\lageq}$$
1251: where $z$ is an auxiliary variable.
1252: \smallskip
1253: {\parindent 20pt
1254: \item{(i)} To find the minimum of $F(f)$, compute $d F\{f+\delta h(x)\}\delta$ where
1255: $h(x)$ is any twice differentiable function.
1256: (This is a standard technique in the calculus of variations).
1257: Integrate by parts the resulting expression (twice) and conclude that the solution
1258: is a piecewise cubic polynomial with knots at the $x_i$.
1259: \index{third derivative}%
1260: Show that the solution and its first two derivatives are continuous at the
1261: knots, but that the third derivative may be discontinuous there.
1262: Show also that the second derivative is zero outside of the range of the
1263: $x_i$.
1264:
1265: \index{least-squares line}%
1266: \item{(ii)} Compute $d F(f)\rho$ and $\partial F(f)/\partial \sigma$ and conclude that either $\hat \rho=0$
1267: implying that the least-squares line is a solution to $\sopt$, or
1268: $\hat\rho\neq 0$, in which case the solution is nonlinear.
1269:
1270: \item{(iii)} Establish an equivalence between the $\sopt$ and the penalized least-squares
1271: problem $\qb$ given in the text.
1272: }
1273:
1274: \exercise Suppose each observation has associated with it a weight $w_i$.
1275: \index{cubic smoothing spline}%
1276: \smallskip
1277: {\parindent 20pt
1278: \item{(i)} Derive the appropriate cubic smoothing spline using a weighted residual sum of squares.
1279: \index{smoothing}%
1280: \item{(ii)} Suppose there are ties in the $x$-values.
1281: Suggest a way to overcome this for smoothing splines.
1282:
1283: \index{ semi-parametric regression}%
1284: }
1285:
1286: \index{matrix}%
1287: \exercise {\sl Semi-parametric regression; ^^{Green, P.J.}^^{Jennison, C. }^^{Seheult, A.} Green, Jennison, and Seheult (1985).} Suppose we have a set of $n$ observations of $p$ predictors arranged in
1288: the $n \times p$ matrix $\bf X$, an additional covariate $\vec z$, and a response vector $\vec y$.
1289: We wish to fit the model $y_i=\vec x^i{\fat\beta} +f(z_i) +\varepsilon_i$ by penalized least-squares.
1290: \index{ estimating equations}%
1291: \index{smoothing}%
1292: Construct an appropriate penalized residual sum of squares, and show that the minimizers must satisfy the following pair of {\sl estimating} equations:
1293: $$\eqalign{ {\fat\beta}&=(\X^T\X)^{-1}\X^T(\vec y-\vec f)\cr
1294: f &= \smooth(\vec y -\X{\fat\beta}\given \vec z)\cr}$$
1295: where $\smooth$ is an appropriate smoothing-spline operator, and $\vec f$ represents the function $f$ evaluated at the $n$ values $z_i$.
1296: \index{ efficient kernel smoothing}%
1297: \exercise {\sl Efficient kernel smoothing; ^{Silverman (1982)}, ^{H\"ardle (1986)}}
1298: \index{kernel smooth}%
1299: Consider the ^^{Nadaraya, E.A.}^^{Watson, G.S.}Nadaraya-Watson form of the kernel smooth given by equation \NadWat:
1300: $$\gsmooth(x_0)={\sum_{i=1}^n\rho({x_0-x_i \over \lambda})y_i\over \sum_{i=1}^n\rho({x_0-x_i \over \lambda})}.$$
1301: \index{fine grid}%
1302: We wish to use the FFT to compute both the numerator and denominator.
1303: With that end in mind we define a fine grid that covers the range of $X$.
1304: Describe how you would approximate the data on the grid and use the FFT to perform the computations.
1305: Pay particular attention to:\smallskip
1306: {\parindent 20pt
1307: \item{(i)} what happens at the boundaries, and
1308: \item{(ii)} a method for overcoming discontinuities due to the approximation.
1309:
1310: }
1311: \exercise Extend the methods described in the previous exercise in order to compute a kernel-weighted running-line fit.
1312: \index{matrix}%
1313: \index{smoothing}%
1314: \exercise Suppose $\bS$ is an $n\times n$ smoother matrix for
1315: smoothing against the $n$ unique predictor values $\vec x=(x_1,\ldots, x_n)^T$.
1316: That is, if $s(x_i)=\sum_{j=1}^n S_{ij}y_j$, then $\bS$ has $ij$th entry $S_{ij}$.
1317: For all of the
1318: smoothers that we have described, $\bS\vec 1=\vec 1$, and for most $\bS\vec
1319: x=\vec x$. For which is the second statement true, and what are the
1320: \index{smoothing splines}%
1321: implications of these two results? For smoothing splines the symmetric
1322: statements are also true: $\bS^T\vec 1=\vec 1$ and $\bS^T\vec x=\vec
1323: x$, since $\bS$ is symmetric in this case. Is this true for any of the
1324: other smoothers? What are the implications of these conditions? How
1325: would you modify $\bS$ for the other smoothers in order that the first
1326: \index{equivalent kernel}%
1327: or both of these latter conditions are satisfied. Try some examples
1328: and see what happens to the local support properties of the equivalent
1329: kernels.
1330: \endexercises
1331: \vfill\supereject
1332: \end
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.