![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to xchap5.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 xchap5.tex CVS log | Up to [Research Unix] / researchv10dc / vol2 / index |
1.1 root 1: \ifnum\pageno=1 \input macros.tex \immediate\openout\inx=chap5.index\makecontents\pageno=105 \fi
2: \proofmodefalse
3: \draft{11/24/89 by Trevor}
4: \def \sH {\script H}
5: \def\SQ{\tilde{\bS}}
6: \def\scM {\script M}
7: \input try.tex
8: \chapter 5 {Some theory for additive models}
9: \index{additive model}%
10: \index{additive model}%
11: \Section{Introduction}
12: In the previous chapter we introduce the nonparametric additive model
13: and the backfitting procedure, a heuristic method for estimation.
14: In this chapter we provide some theoretical underpinning for these ideas.
15: The technical content and level of this chapter is somewhat
16: higher than the others, so the reader
17: interested only in applications
18: may well decide not to tackle it.
19:
20: The chapter has two main parts.
21: \index{backfitting algorithm}%
22: \index{additive model}%
23: In the first half the backfitting algorithm is justified as a method
24: for estimating the additive model.
25: Several justifications are provided.
26: \index{backfitting algorithm}%
27: We first introduce an $L_2$ version of backfitting, that is,
28: a backfitting algorithm for square integrable random variables, and
29: show how the intuitive procedure introduced in the last chapter
30: can be viewed as a data analogue of this.
31: A second, different, justification for backfitting,
32: comes from a penalized least-squares framework.
33: Unlike the $L_2$ argument, it makes no appeal to random
34: variables:
35: instead, it applies to finite-sample backfitting procedures that use linear smoothers
36: (recall the discussion of linear smoothers in sections~2.8 and 3.4.2).
37:
38: \index{estimating equations}%
39: An important by-product of these justifications
40: is the set of {\em estimating equations} that are solved by backfitting.
41: \index{linear system}%
42: This linear system can be solved noniteratively and in some
43: special cases such a direct solution is more appropriate
44: than backfitting.
45:
46: \index{reproducing-kernel hilbert-spaces}%
47: \index{estimating equations}%
48: \index{backfitting algorithm}%
49: We also describe very briefly a more technical abstract derivation using the theory of reproducing-kernel Hilbert-spaces.
50:
51: In the second half of the chapter
52: we study the existence and uniqueness of the solutions to the additive-model
53: estimating equations and
54: the convergence of the backfitting algorithm.
55: The results pertain to linear smoothers only.
56: \index{ concurvity}%
57: One of the interesting notions that arises is {\em concurvity}, the
58: analogue of collinearity.
59: \index{standard-error bands}%
60: We also discuss some theory for standard-error bands and degrees of freedom of the estimated
61: smooths,
62: and the relationship of backfitting to the Gram-Schmidt and Gauss-Seidel
63: \index{generalized additive model}%
64: techniques.
65:
66: Related theory for generalized additive models
67: is not covered here but
68: is discussed
69: in the
70: next chapter.
71: \index{additive model}%
72: We don't devote much space to asymptotic issues such as consistency or rates of convergence for
73: additive models; these are briefly mentioned in the
74: bibliographic notes.
75:
76: \index{estimating equations}%
77: \index{additive model}%
78: \index{estimating equations}%
79: \index{additive model}%
80: \index{backfitting algorithm}%
81: \Sectionskip\Section{Estimating equations for additive models}
82: \Mark{ESTIMATING EQUATIONS}
83: The additive model
84: $$\ev(Y\given \rvX)=\sum_{j=1}^p f_j(X_j)\eqn{\addm}$$
85: can be estimated by the backfitting algorithm, which we
86: give again below:
87:
88: \setbox1=\vbox{\hsize \algwidth {\setnine\parindent 20pt
89: \item{(i)} {\em Initialize:}$ \quad f_j=f_j^{\,(0)},
90: j=1,\ldots, p$
91: \item{(ii)} {\em Cycle:} $\quad j=1,\ldots, p,1,\ldots
92: p,\ldots$
93: $$f_j=\smooth_j(\vec y-\sum_{k\neq j}\vec f_k\given \vec x_j) $$
94: \item{(iii)} Continue (ii) until the individual functions don't
95: \index{backfitting algorithm}%
96: change.
97: }%end algorithm
98: \smallskip
99: } %end box1
100: \midinsert \algorithm{{\ninerm\noindent Algorithm \chapnodot 1}
101: The backfitting algorithm}{\box1} \endinsert
102:
103: In the above, the $\vec f_j$ are the $n$-vectors
104: $\{f_j(x_{1j}),\ldots, f_j(x_{nj})\}^T$, with $x_{ij}$ in the order of $y_i$.
105: \index{constant term}%
106: We have omitted the constant term $\alpha$ in \addm; we see later that
107: this does not change the resulting estimates.
108:
109: \index{scatterplot smoother}%
110: In order to justify this procedure, we need some way of introducing
111: the smoothness that is provided by the scatterplot smoothers
112: in the algorithm.
113: To put it another way,
114: if we naively
115: tried to minimize
116: $$\sum_{i=1}^n \Bigl\{y_i-\tsum_{j=1}^p f_j(x_{ij})\Bigr\}^2\eqn{\naive}$$
117: then the solution would be any set of functions $(\,f_j:j=1,\ldots,p\,)$ that interpolated the data (assuming for the moment that the $X$-values are distinct).
118: For example, $f_1(x_{i1})=y_i$ $\forall i$ and $f_j\equiv 0$
119: for $j >1$.
120: We discuss several ways of introducing smoothness.
121: One approach is explicitly to add a term to $\naive$ that penalizes for lack of
122: smoothness, the
123: {\em penalized least-squares} approach.
124: Another approach, described next, is to step back and
125: consider random variables instead of data.
126: \index{additive model}%
127: \index{conditional expectation operators}%
128: A Hilbert-space version of the additive model and backfitting can be
129: formulated,
130: with conditional expectation operators playing the role of smoothers.
131: Besides its use here,
132: this formulation is
133: mathematically
134: interesting
135: in its own right.
136: \index{conditional expectation}%
137: The data version of backfitting is then derived as an empirical version
138: of this Hilbert-space procedure, the smoothness
139: entering when one considers how best to estimate the conditional
140: expectations.
141: \index{reproducing-kernel hilbert-spaces}%
142: A third approach, based on reproducing-kernel Hilbert-spaces, is a more abstract version of the penalized least-squares approach. Later on in the chapter
143: we give yet another derivation, based on a Bayesian stochastic model.
144: \index{hilbert space}%
145:
146: \sectionskip\section{$L_2$ function spaces}
147: Let $\sH_j$ for $j=1,\ldots,p$ denote the Hilbert spaces of measurable
148: functions $\phi_j (X_j)$ with $\ev\phi_j(X_j)=0$, $\ev\phi^2_j(X_j) <
149: \infty$, and inner product
150: $\inner{\phi_j(X_j),\phi_j'(X_j)}=\ev \phi_j(X_j)\phi_j'(X_j) $. In
151: addition, denote by $\sH$ the space of arbitrary centered, square
152: integrable functions of $X_1,\ldots,X_p$. We consider the $\sH_j$
153: as subspaces of $\sH$ in a canonical way. Furthermore, denote by
154: $\sH^{add} \subset \sH$ the linear subspace of additive
155: functions: $\sH^{add}=\sH_1+\cdots +\sH_p$, which is closed under some technical assumptions.
156: These are all subspaces of $\sH_{YX}$,
157: the space of centered square integrable functions of $Y$ and $X_1,\ldots,X_p$.
158:
159: The optimization problem in this population setting is
160: to minimize
161: $$\ev\{Y-g(\rvX)\}^2\eqn{\backc}$$
162: over $g(\rvX)=\sum_j f_j(X_j)\in\sH^{add}$.
163: Of course, without the additivity
164: restriction, the solution is simply $\ev(Y\given \rvX)$; we seek the
165: closest additive approximation to this function. Since by assumption $\sH^{add}$ is a
166: closed subspace of $\sH$ this minimum exists and is unique; the
167: individual functions $f_j(X_j)$, however, may not be uniquely
168: \index{conditional expectation}%
169: \index{orthogonal projection}%
170: determined. Denote by $P_j$ the conditional expectation operator $\ev(\cdot\given
171: X_j)$; as such $P_j$ is an orthogonal projection onto $H_j$.
172:
173: The minimizer $g(\rvX)$ of \backc\ can be characterized by residuals $Y-g(\rvX)$ which are orthogonal to the space of fits:
174: $Y-g(\rvX) \perp \sH^{add}$.
175: Since $\sH^{add}$ is generated by $\sH_j\ (\subset \sH^{add})$,
176: we have equivalently:
177: $Y-g(\rvX)\perp \sH_j,\quad\forall j$
178: or:\
179: $\;P_j\{Y-g(\rvX)\}=0\quad\forall j$.
180: Component-wise this can be written as
181: $$\eqalign{f_j(X_j) &=P_j\Bigl\{Y-\sum_{k\neq j} f_k(X_k)\Bigr\}\cr &=
182: \ev\Bigl\{Y-\sum_{k\neq j}f_k(X_k)\given X_j\Bigr\}.\cr}\eqn{\backe}$$
183: \index{estimating equations}%
184: Equivalently, the following system of {\em estimating equations}
185: is necessary and sufficient for $\vec f=
186: (f_1,\ldots,f_p)$ to minimize \backc:
187: $$\pmatrix{I&P_1&P_1&\cdots&P_1\cr
188: P_2&I&P_2&\cdots&P_2\cr\vdots&\vdots&\vdots&\ddots&\vdots\cr
189: P_p&P_p&P_p&\cdots&I\cr}\pmatrix{f_1(X_1)\cr f_2(X_2)\cr\vdots\cr
190: f_p(X_p)}= \pmatrix{P_1 Y\cr P_2Y\cr \vdots\cr P_pY}\eqn{\backd}$$
191: or
192: $$\bP\vec f=\vec Q Y,$$
193: where $\vec P $ and $ \vec Q$ represent a matrix and vector of operators, respectively,
194: and operator matrix multiplication is defined in the obvious way.
195: \index{linear system}%
196:
197: The reader familiar with the Gauss-Seidel method for solving
198: linear systems of equations will recognize backfitting as a formal
199: Gauss-Seidel algorithm for solving the system $\backd$.
200: \index{conditional expectation operators}%
201: (We say {\em formal} because the elements in the left matrix
202: of $\backd$ are not real numbers or real-valued matrices but
203: \index{linear system}%
204: conditional expectation operators.)
205: Suppose we wish to solve a linear system of equations $\bA\vec z=\vec b$, with
206: $\bA=\{a_{ij}\}$ an $m\times m$ matrix, $\vec b$ an $m$ vector and $\vec z$
207: the vector of $m$ unknown coefficients.
208: The Gauss-Seidel iterative method
209: solves for each $z_i$ in turn from the
210: relation in the $i$th row
211: $\sum_{j=1}^m a_{ij} z_j=b_i$.
212: This process is repeated for $i=1,\ldots, m,1,\ldots,m,\ldots$,
213: using the latest values of each $z_j$ at each step, until convergence.
214:
215: The connection between backfitting and the Gauss-Seidel method becomes more precise
216: when we consider the corresponding data version of $\backd$ using linear smoothers.
217: \index{ smoother matrix}%
218: Recall from section~2.8 that a linear smoother can be written as
219: a {\em smoother matrix} times the response vector $\vec y$,
220: that is $\hatvec f=\bS\vec y$.
221: \index{smoothing spline}%
222: \index{kernel smoother}%
223: Examples of linear smoothers include the running-mean, locally-weighted running-line,
224: smoothing splines and kernel smoothers.
225: \index{backfitting algorithm}%
226: \index{conditional expectation}%
227: \index{linear scatterplot smoother}%
228: \index{smoother matrix}%
229: Consider then a backfitting algorithm that estimates the conditional
230: expectation operator $P_j$ by a linear scatterplot smoother
231: with smoother matrix $\bS_j$.
232: \index{estimating equations}%
233: Then the data version of the estimating equations $\backd$ is
234: the $np\times np$ system
235: $$\pmatrix{\bI&\bS_1&\bS_1&\cdots&\bS_1\cr
236: \bS_2&\bI&\bS_2&\cdots&\bS_2\cr\vdots&\vdots&\vdots&\ddots&\vdots\cr
237: \bS_p&\bS_p&\bS_p&\cdots&\bI\cr}\pmatrix{\vec f_1\cr \vec f_2\cr\vdots\cr \vec
238: f_p}= \pmatrix{\bS_1 \vec y\cr \bS_2\vec y\cr \vdots\cr \bS_p\vec
239: y}.\eqn{\backdd}$$
240:
241: In short form we write
242: $$\hat {\vec P}\vec f=\hat {\bQ}\vec y.$$
243: Backfitting is a Gauss-Seidel procedure for solving the above
244: system.
245: The only nonstandard aspect is that we solve for $n$ elements at each step instead of one (block Gauss-Seidel), although for linear smoothers this distinction can be dropped.
246:
247: Now
248: suppose we start with $\backdd$.
249: Why use an iterative procedure like backfitting to find its solution?
250: Why not use a standard, noniterative method like a QR decomposition?
251: The difficulty is that in general, $\backdd$ is an $np\times np$ system
252: and since methods like QR require $O(m^3)$ operations to solve an $m\times m$
253: system, our problem would cost $O\{(np)^3\}$ operations.
254: \index{smoothing spline}%
255: On the other hand, backfitting exploits the special structure in $\backdd$,
256: and if the smoothers can be applied in $O(n)$ computations (as is the
257: case for running-lines and smoothing splines), then backfitting
258: requires only $O(np)$ computations.
259: (This assumes that a fixed number of iterations is sufficient for convergence).
260: If, however, the effective dimension of the system $\backdd$ is really less than $np$,
261: there may be better methods than backfitting for solving the problem.
262: \index{orthogonal projection}%
263: In particular, if each $\bS_j$ is an orthogonal projection and the
264: union of the projection spaces has rank $m$, then $\backdd$ is equivalent
265: to an $m\times m$ least-squares problem (Exercise~5.3) and least-squares methods
266: are likely to be preferable if $m << n$.
267: This is the case if the $\bS_j$s produce linear or polynomial fits,
268: or if we use regression splines with a small number of knots.
269:
270: \index{estimating equations}%
271: Later in this chapter the properties of the estimating equations \backdd\ are
272: studied.
273: We find that there is an intimate connection between the
274: existence of solutions of this system and the convergence of the
275: backfitting procedure for finding these solutions.
276:
277: \sectionskip\section{Penalized least-squares}
278: In this section we provide a different justification for backfitting
279: from that given in
280: the Hilbert-space framework in the previous section.
281: \index{cubic smoothing spline}%
282: \index{penalized least-squares criterion}%
283: In section~2.10 we derive the cubic smoothing spline
284: as the minimizer over all twice continuously differentiable functions of the
285: penalized least-squares criterion
286: $$\sum_{i=1}^n\{y_i-f(x_i)\}^2+\lambda\int \{f''(x)\}^2 \, dx.\eqn{\qb}$$
287: We establish this by using the fact that the solution to $\qb$ is a cubic
288: spline,
289: and hence we simplified $\qb$ by writing it as a function of $f_i=f(x_i)$, the $n$ evaluations of the minimizing function~$f$.
290: This gave the equivalent form
291: $$(\vec y-\vec f)^T(\vec y-\vec f) +\lambda\vec f^T\bK\vec f\eqn{\nqd}$$
292: where $\bK$ is a certain quadratic penalty matrix.
293: The quantity
294: $\nqd$ is easily shown to have a minimum given by
295: $$\hat\vec f=(\bI+\lambda \bK)^{-1}\vec y.
296: \eqn{\css}$$
297: We also argued in the opposite direction for
298: other symmetric linear smoothers.
299: \index{smoother matrix}%
300: That is, given a symmetric linear smoother based on the smoother matrix
301: $\bS$, the smooth $\hat\vec f=\bS\vec y$ minimizes
302: $$(\vec y-\vec f)^T(\vec y-\vec f) +\vec f^T(\bS^- -\bI)\vec f\eqn{\newm}$$
303: over all $\vec f\in\script{R}(\bS)$ (the range of $\bS$), where $\bS^-$ is any generalized inverse of $\bS$.
304:
305: \index{additive model}%
306: In order to extend this idea to the estimation of the additive model,
307: we generalize the criterion $\qb$ in an obvious way.
308: We seek to minimize
309: $$\sum_{i=1}^n\Bigl\{y_i-\tsum_{j=1}^pf_j(x_{ij})\Bigr\}^2 +\sum_{j=1}^p
310: \lambda_j\int \{f_j''(t)\}^2 \, dt\eqn{\splinpen}$$
311: over all twice
312: continuously differentiable functions $f_j$.
313: Before deriving the solution to $\splinpen$, let's take note of some
314: of its features.
315: Notice that each function in
316: $\splinpen$ is penalized by a separate constant $\lambda_j$.
317: This in turn determines the smoothness of that function in the solution.
318: Note also that if the $\lambda_j$s are all zero (no smoothness penalty)
319: the solution to \splinpen\ is any interpolating set of functions whose evaluations satisfy $\sum_{j=1}^p f_j(x_{ij})=y_i$ for $i=1,\ldots,n$.
320: On the other hand, if each $\lambda_j$ goes to infinity, the
321: penalty term goes to infinity unless $f_j''(t)=0$ for all $j$, that is,
322: unless
323: each $f_j$ is linear.
324: Hence the problem reduces to standard linear least-squares.
325:
326: Using a straightforward
327: extension of the arguments used in the single-predictor
328: case,
329: the solution to \splinpen\ is shown to be a cubic spline in each of the predictors.
330: As before we parametrize by evaluations at the $n$ observations.
331: Thus we may rewrite $\splinpen$ as
332: $$\Bigl(\vec y-\sum_{j=1}^p \vec f_j\Bigr)^T\Bigl(\vec y-\sum_{j=1}^p\vec f_j\Bigr)+\sum_{j=1}^p
333: \lambda_j \vec f_j^T \bK_j \vec f_j \eqn{\nsplinpen}$$
334: where the $\bK_j$s are penalty matrices for each predictor, defined
335: analogously to the $\bK$ for a single predictor given in section~2.10.
336: Now if we differentiate $\nsplinpen$ with respect to the function $\vec f_k$
337: we obtain $-2(\vec y-\sum_j \vec f_j)+2\lambda_k \bK_k\vec f_k=\bf 0$
338: or
339: \index{smoother matrix}%
340: \index{cubic smoothing spline}%
341: \index{estimating equations}%
342: $$\hat\vec f_k=\Bigl(\bI+\lambda_k \bK_k\Bigr)^{-1}\Bigl(\vec y-\tsum_{j\neq k} \hat\vec f_j\Bigr).\eqn{\solsplin}$$
343: As noted earlier, $(\bI+\lambda_k \bK_k)^{-1}$ is the smoother matrix for
344: a cubic smoothing spline, and hence $\solsplin$, for $k=1,\ldots, p$,
345: are just
346: the estimating equations $\backdd$.
347:
348: \index{penalized least-squares criterion}%
349: \index{estimating equations}%
350: Arguing in the opposite direction, the minimizers of the
351: penalized least-squares criterion
352: $$\Bigl(\vec y-\sum_{j=1}^p \vec f_j\Bigr)^T\Bigl(\vec y-\sum_{j=1}^p\vec f_j\Bigr)+\sum_{j=1}^p
353: \vec f_j^T(\bS_j^{-}-\bI)\vec f_j,\eqn{\gsplinpen}$$
354: over all $\vec f_j\in\script{R}(\bS_j)$,
355: are the solutions to the estimating equations $\backdd$ (Exercise~5.1).
356:
357:
358:
359: As we do in the single predictor case (Exercise~3.6), we can interpret
360: each of the penalty terms in $\nsplinpen$ as a down-weighting of each
361: of the components of $\vec f_j$, the down-weighting determined by
362: the corresponding eigenvalue of that component
363: and
364: $\lambda_j$.
365:
366: \index{reproducing-kernel hilbert-spaces}%
367: \index{additive model}%
368: \def\sH{{\script H}}
369: \sectionskip
370: \section{Reproducing-kernel Hilbert-spaces}
371: This section describes a more abstract framework for defining and estimating general nonparametric regression models which includes additive models as a special case.
372: We present these results to give the reader a taste of this rich
373: area;
374: the level of mathematics is somewhat higher than the rest of the chapter.
375: The description is close to that of ^{Chen, Gu and Wahba (1989)}.
376:
377: \index{hilbert space}%
378: A Hilbert space $\sH$ of real-valued functions of $t\in \Omega$ is a {\em reproducing-kernel} Hilbert-space if evaluation is a continuous linear functional.
379: \index{ representers of evaluation}%
380: By the Riesz representation theorem, there exist {\em representers of evaluation}
381: $e_t\in\sH$ such that $f(t)=\langle f,e_t\rangle_\sH$ for $f\in \sH$, where $\langle \cdot,\cdot\rangle_{\script H}$ denotes the inner-product on $\sH$.
382: The consequences of these properties will become clearer as we proceed.
383:
384: The reproducing kernel itself,
385: $Q(\cdot,\cdot):\Omega\times\Omega\mapsto\R 1$, is defined by
386: $Q(s,t)=\langle e_s,e_t\rangle_\sH$, and consequently
387: $e_s=Q(s,\cdot)$, considered as a function of the second argument, with
388: the first held fixed at $s\in \Omega$.
389: We will see that the kernel $Q$, evaluated at the realizations of $t$, provides a finite dimensional basis for representing the solution to a class of optimization problems.
390:
391: Now suppose $\Omega$ is a space of vector predictors $\fat{X}=(X_1,\ldots,X_p)$ and that $\sH$ has the decomposition
392: $$
393: \sH=\sH_0+\sum_{k=1}^q\sH_k,
394: $$
395: where $\sH_0$ is spanned by $\phi_1,\ldots,\phi_M$, and $\sH_k$ has the reproducing kernel $Q_k(\cdot,\cdot)$.
396: The space $\sH_0$ is the projection component of $\sH$, that is, the
397: space of functions that are not to be penalized in the optimization.
398: In the previous section $\sH_0$ is the space of functions linear in $t$.
399:
400: We are now set up to pose the optimization problem.
401: For a given set of predictors $\vec x^1,\ldots,\vec x^n$ (with each $\vec x^i\in\Omega$), find $f=\sum_{k=0}^q f_k$ with $f_k\in \sH_k$, \ $k=0,\ldots, q$, to minimize
402: $$
403: \sum_{i=1}^n\Bigl\{y_i-\tsum_{k=0}^qf_k(\vec x^i)\Bigr\}^2+ \sum_{k=1}^q\lambda_k\norm{f_k}^2_{\sH_k}.\eqn{\rkcrit}
404: $$
405:
406: The first part of the criterion is discrete in nature, and is the reason why
407: reproducing-kernel spaces are natural for these kinds of problems.
408: We do not want small changes in the $\vec x^i$ to result in vastly different solutions; this is why it is desirable for evaluation to be continuous.
409:
410:
411: The theory of reproducing kernels guarantees that a minimizer exists, and has the form
412: $$\eqalign{
413: \hat{f}_0(\fat{X})&=\sum_{j=1}^M\beta_{j0}\phi_j(\fat{X})\cr
414: \hat{f}_k(\fat{X})&=\sum_{i=1}^n\beta_{ik}Q_k(\fat{X},\vec x^i).\cr
415: }
416: \eqn{\rksoln}
417: $$
418: Furthermore, if the projection onto $\sH_0$ is unique, then so is the solution
419: to the larger problem.
420: So even though the problem is posed in an infinite-dimensional space, the minimizing $\hat f$ is finite-dimensional, and the $Q_k$ supply bases for representing the solution.
421: The parameters are found by minimizing the finite dimensional quadratic criterion
422: $$
423: \norm{\vec y-\vec T\fat{\beta_0}-\sum_{k=1}^q\vec Q_k\fat{\beta}_k}^2 +
424: \sum_{k=1}^q\lambda_k\fat{\beta}_k^T\vec Q_k\fat{\beta}_k\eqn{\rkdata}
425: $$
426: where $\vec T$ is the $n\times M$ matrix of evaluations of $\phi_j$, with
427: $ij$th entry $T_{ij}=\phi_j(\vec x^i)$, and $\vec Q_k$ is the $n\times n$
428: { evaluated} kernel with $ij$th entry $Q_k(\vec x^i,\vec x^j)$.
429: This problem is of dimension at most $qn+M$.
430:
431: \index{additive model}%
432: At this point a number of specializations are possible:
433: \smallskip
434: {\parindent 20pt
435: \item{(i)} If $q=p$ and each of the $\sH_k$ are the canonical subspaces of $\sH$,
436: \index{cubic smoothing spline}%
437: then the additive model consists of a sum of univariate functions. Furthermore, by choosing an inner product appropriate for cubic smoothing splines, the problem reduces exactly to \splinpen, although the solution is typically represented by a different basis.
438: \index{estimating equations}%
439: \item{(ii)} The current specification has $nq+M$ parameters. If attention is restricted to $f_0$ and $f_+=\sum_{k=1}^qf_k$, with $Q_+=\sum_{k=1}^q Q_k/\lambda_k$, then the dimension of the solution is reduced to $M+n$.
440: \item{(iii)} The general problem as specified by \rksoln\ can potentially be solved more cheaply by backfitting; the system of estimating equations that characterize the minimum of \rksoln\ can be written in a form similar to \backdd:
441: $$\pmatrix{\bI&\bS_0&\bS_0&\cdots&\bS_0\cr
442: \bS_1&\bI&\bS_1&\cdots&\bS_1\cr\vdots&\vdots&\vdots&\ddots&\vdots\cr
443: \bS_q&\bS_q&\bS_q&\cdots&\bI\cr}\pmatrix{\vec f_0\cr \vec f_1\cr\vdots\cr \vec
444: f_q}= \pmatrix{\bS_0 \vec y\cr \bS_1\vec y\cr \vdots\cr \bS_q\vec
445: y}\eqn{\backrk}$$
446: where $\bS_0=\bT(\bT^T\bT)^{-1}\bT$ and $\bS_k=\bQ_k(\bQ_k+\lambda_i\bI)^{-1}$.
447: If the computational complexity of the individual operators $\bS_k$ is significantly lower than that of the full problem, savings can be made using backfitting-type algorithms.
448:
449: }\smallskip
450: \index{estimating equations}%
451: \index{backfitting algorithm}%
452: \index{additive model}%
453: This is a very brief summary of some powerful machinery; we cite a number of relevant references in the bibliographic section for more details of this approach and for pointers to the large application area.
454:
455: \Sectionskip\Section{Solutions to the estimating equations}
456: %\Mark{THE BACKFITTING ALGORITHM}
457: \section{Introduction}
458: The remainder of the chapter
459: focuses attention on additive models with linear smoothers $\smooth_1,\ldots,\smooth_p$, and the algorithms for estimating them.
460:
461: Most of the smoothers that we have discussed produce function estimates, and so we could discuss issues such as convergence in terms of these functions as well.
462: Instead we restrict attention to the evaluation of these functions at the $n$ realizations of the predictors.
463: We do this mainly for simplicity and clarity, but point out that most of the results cited here, in particular those pertaining to convergence, can be modified to include this more general case.
464: As a consequence, we usually refer to a smoother by its matrix representation $\bS$ rather than in the operator form $\smooth$.
465:
466: \index{estimating equations}%
467: The centerpiece of the discussion is the set of estimating equations \backdd.
468: Before one delves into methods for solving such a system, questions
469: of consistency and degeneracy have to be answered.
470: In other words, we
471: must confirm that the system has at least one solution, and find out
472: whether this solution is unique.
473: We first look at a few special cases which are are easy to
474: work out and illustrate the
475: main issues.
476: Later in the chapter we answer these questions in some generality.
477:
478:
479: In our discussion we sometimes assume implicitly that the same
480: smoother is used for each of the variables, but this is only for ease
481: of presentation. The results are general in nature and apply to any
482: backfitting procedure in which some linear smoother is used for each
483: of the variables.
484: In fact, there is no need even to assume that each smoother is based
485: on a single predictor:
486: for example a two-dimensional smoother or a least-squares fit on some
487: set of predictors could be included.
488: Even more generally
489: one can think of $\bS_1,\ldots,\bS_p$ as a
490: set of linear transformations, without reference to predictor variables at all.
491:
492: Throughout the chapter we assume
493: that the $\bS_j$s all reproduce constant functions.
494: \index{backfitting algorithm}%
495: Note that this
496: causes a simple kind of non\-unique\-ness in the backfitting algorithm.
497: \index{starting functions}%
498: Suppose the starting functions are all zero.
499: Then at every stage of the procedure $\hat \vec f_1$ has
500: the same mean as $\vec y$, but the other $\hat\vec f_j$s have mean $\bf 0$.
501: If the procedure started at $j=2$, however, the mean of $\vec y$ would go into
502: $\hat\vec f_2$ instead.
503: \index{starting functions}%
504: \index{starting functions}%
505: A closer look reveals that
506: nonzero starting functions cause a dependence
507: of the final iterates on the values of the starting functions.
508: \index{additive model}%
509: It is also clear that unless special constraints are built in, the constant in the additive model is not identifiable.
510: \index{ concurvity}%
511: \index{collinearity}%
512: This is a special instance of what we call {\em concurvity}, the analogue of {\em collinearity} in linear models.
513:
514: It turns out that such degeneracies do not affect convergence in any important way.
515: In this case a simple fix is possible:
516: assume that $\vec y$ has been centered to have mean $\bf 0$,
517: and replace $\bS_j$ by
518: the matrix that smooths then subtracts off the average of the smooth.
519: \index{smoother matrix}%
520: \index{centered smoother}%
521: It is easy to see that the resultant smoother matrix is
522: $(\bI-{\bf 1}{\bf 1}^T/n)\bS_j$,
523: what we call a {\em centered} smoother.
524: This ensures that at every stage of the procedure the $\hatvec f_j$s have
525: mean $\bf 0$.
526:
527:
528: \index{additive model}%
529: \sectionskip\section{Projection smoothers}
530: The additive model is introduced as a generalization of the linear regression model.
531: \index{orthogonal projection}%
532: What if we use a linear least-squares fit,
533: or any other orthogonal projection,
534: for each predictor?
535: \index{estimating equations}%
536: As mentioned in section~5.2.1, the set of estimating equations \backdd\ is
537: equivalent to the usual normal equations for linear regression
538: (Exercise~5.3).
539: \index{orthogonal projection}%
540: Thus if
541: $\bS_j$ is an orthogonal projection
542: and $\script L_{col}(\bS_j)$ denotes the subspace spanned by the columns
543: of $\bS_j$,
544: we expect
545: the backfitting solution $\hatvec y=\sum_1^p \hatvec f_j$ to converge to
546: the projection of $\vec y$ onto $V=\script L_{col}(\bS_1)\oplus
547: \script L_{col}(\bS_2)\oplus \ldots\oplus \script L_{col}(\bS_p)$.
548: Indeed if this wasn't the case we might well question the entire
549: backfitting paradigm.
550: Fortunately, it is fairly easy to show that backfitting does the
551: expected in this special case.
552: We sketch the proof here, leaving the details to the
553: reader (Exercise~5.2).
554: \index{backfitting converges}%
555: The idea is to show that the residual vector from backfitting
556: converges to the least-squares residual vector,
557: that is, the projection of $\vec y$ onto the orthogonal
558: complement of $V$.
559: After one cycle of backfitting, the residual vector
560: from backfitting, say
561: $\vec r$, is $\bC\vec y$ where
562: $$\bC=(\bI-\bS_p)(\bI-\bS_{p-1})\cdots (\bI-\bS_1).\eqn{\resv}$$
563: Thus after $m$ cycles the residual is $\vec r^{\,(m)}=\bC^m\vec y$.
564: We can split $\vec y$ into its components in the
565: projection space $V$ and its orthogonal complement,
566: that is $\vec y=\hat\vec y+\vec y^{\perp}$.
567: Now the operator $\bC$ {\em takes residuals} along each predictor in turn, and
568: hence leaves $\vec y^{\perp}$ unchanged.
569: Thus
570: $$\vec r^{\,(m)}=\bC^m\hat\vec y+\vec y^{\perp}.\eqn{\splitup}$$
571: The proof is then completed by showing that $\norm{\bC^m\hat\vec y}
572: \rightarrow 0$ and hence $\vec r^{\,(m)}\rightarrow \vec y^{\perp}$ (Exercise~5.2).
573:
574: Hence we see that the backfitting procedure provides an alternative
575: method for computing least-squares fits.
576: Examination of $\splitup$ reveals that it works by successively projecting
577: the current residual into the space orthogonal to each $\script L_{col}(
578: \bS_j)$.
579: \par
580: Figure~\zigfig\ gives a picture of this
581: in the two-predictor case.
582: Each ${\vec x}_j^\perp$ denotes the vector orthogonal to $\vec x_j$ in the span of $\script L_{col}(\vec x_1,\vec x_2)$.
583: It shows the backfitting residual $\vec r^{\,(m)}$ converging to the least-squares residual $\vec y^{\perp}$ in a zig-zag fashion.
584: At convergence the backfitting residual vector is orthogonal to
585: each $\script L_{col}(\bS_j)$ and equals the least-squares residual
586: $\vec y^{\perp}$.
587: This is a novel but not very practical way of finding the least-squares fit, for it can be very slow if the predictors are correlated.
588: In particular, one can show that
589: for two predictors
590: the difference between the $i$th iterate and the solution converges to zero
591: geometrically at rate
592: $\cos(\theta)$, where $\theta$ is the angle between
593: the two predictor vectors
594: (Exercise~5.4).
595: This is intuitively plausible from Fig.~\zigfig.
596: \index{log-linear models}%
597: There is a close connection between backfitting and the iterative proportional scaling algorithm for fitting log-linear models to contingency tables (^{Bishop \etal,~1975}).
598: \index{estimating equations}%
599: They both cycle through a system of estimating equations and update one component at a time; convergence is geometric in both cases.
600:
601: It is interesting to look at two extreme situations that can occur.
602: First, suppose that the predictors are perfectly
603: collinear
604: ($\theta=0$ in the two predictor case).
605: \index{backfitting converges}%
606: Then one can easily check that backfitting converges after a single
607: cycle, with $\hat\vec f_1=\hat\vec y$ and $\hat\vec f_j =\bf 0$ for $j>1$.
608: More interestingly, suppose that the predictors are mutually uncorrelated
609: ($\theta=90^\circ$ in the two-predictor case).
610: Then again we have convergence after a single cycle.
611: \index{gram-schmidt method}%
612: This brings up the connection of backfitting with the Gram-Schmidt method
613: for solving the least-squares normal equations.
614: \index{gram-schmidt method}%
615: Like backfitting, the Gram-Schmidt method works by regressing the
616: current residual onto each predictor in turn.
617: However there is one important difference.
618: After regressing on the $j$th predictor, the $(j+1)$th through $p$th
619: predictors are also regressed
620: on the $j$th predictor and the residual vector from each regression
621: is used in place of the predictor in the remaining steps.
622: This process orthogonalizes the predictors and because of this,
623: the Gram-Schmidt procedure converges after a single
624: cycle.
625: \index{gram-schmidt method}%
626: \index{backfitting algorithm}%
627: In the orthogonal-projection setting, clearly the Gram-Schmidt method is
628: superior to backfitting and this suggests that a sweeping-out operation be used to improve the backfitting
629: algorithm with general smoothers.
630: This is not useful, however, because the
631: resultant model would no longer be additive in the predictors.
632:
633: We note also that the above proof of the convergence of the residual
634: vector does not establish that the estimated functions converge to
635: the correct ones.
636: Indeed, if there exists strict collinearity among the predictors
637: it wouldn't be clear to which solution the functions $\hat \vec{f}_j$ produced by backfitting
638: would converge, if indeed they converge at all.
639: \index{starting functions}%
640: It turns out that backfitting always does converge to a solution representing
641: the projection of $\vec y$ onto $V$, but
642: the fixed point depends on the starting functions.
643: \index{additive model}%
644: \index{semi-parametric model}%
645:
646: \sectionskip\section{Semi-parametric models}
647: Consider an additive model in which all but one term is assumed
648: \index{backfitting algorithm}%
649: to be linear --- the so called {\em semi-parametric} model.
650: The backfitting algorithm for estimating such
651: a model can be thought
652: of as having two smoothers: one a projection $\bS_1=\bX(\bX^T\bX)^{-1}\bX^T$ producing a least-squares fit
653: $\bX\hatfat \beta$ on one or more covariates (represented by the full-rank design
654: matrix $\bX$), and the other a smoother $\bS_2$ producing an estimate $\hat\vec
655: f_2$. The backfitting steps are
656: $\vec f_1=\bS_1(\vec y-\vec f_2)=\bX(\bX^T\bX)^{-1}\bX^T(\vec y-\vec f_2)\equiv
657: \bX {\fat\beta}$, and $\vec f_2=\bS_2(\vec y-\bX{\fat\beta})$.
658: It turns out that we can solve for $\hatfat\beta$ and
659: $\hat\vec f_2$ explicitly (Exercise~2.8):
660: $$\eqalign{{\hatfat\beta}
661: &=\{\bX^T(\bI-\bS_2)\bX\}^{-1}\bX^T(\bI-\bS_2)\vec y\cr
662: { \hat\vec f}_2 &= \bS_2(\vec
663: y-\bX\hatfat\beta)\cr}\eqn{\simmm}$$
664: so that iteration is unnecessary.
665: Although $\bS_2$ is an $n\times n$ matrix, all we have to do is smooth each of
666: the $p$ columns of $\bX$, an operation that can usually be performed in $O(np)$ operations.
667: \index{smoother matrix}%
668: \index{weight matrix}%
669: This provides a computationally simple method for nonparametric
670: analysis of covariance; it is interesting that the smoother matrix
671: $\bS_2$ enters as the weight matrix for the regression on $\bX$.
672: \index{estimating equations}%
673: \index{unique solution}%
674: This manipulation also shows that in this special case, the estimating equations
675: $\backdd$ are consistent and have a unique solution as long a
676: $\bX^T(\bI-\bS_2)\bX$ is invertible.
677: In section~6.7
678: we discuss this model in more detail.
679:
680: \index{additive model}%
681: \sectionskip\section{Backfitting with two smoothers}
682: The third special case that we consider is an additive model
683: that involves two linear smoothers.
684: It turns out that one can analyse the properties of the estimating
685: equation solutions and the convergence of backfitting with some
686: fairly elementary calculations, and this exercise sheds light on the main
687: issues that arise in the more difficult $p$-predictor case.
688:
689:
690:
691: It is possible to determine general conditions under which the system \backdd\ is consistent by
692: checking whether the rank of $\hatvec P$ is the same as that of the
693: augmented matrix $[\hatvec P\colon\hatvec Q\vec y]$.
694: However, it is easier to proceed by constructing the solutions to \backdd\
695: through the backfitting procedure.
696:
697: Recall that the components of each $\vec x_j$ are in the same order as
698: the components of $\vec y$.
699: As a technical point, this means that
700: $\bS_j$ really means
701: $\bE_j^{-1}\bS_j\bE_j$ where $\bE_j$ is the permutation matrix that sorts
702: in the order of $\vec x_j$ (Exercise~5.19).
703:
704: Let
705: $\norm{\bC}=\sup_{\vec a\neq \bf 0}\norm{\bC\vec a}/\norm{\vec a}$, the
706: 2-norm of the matrix $\bC$
707: (this choice is made for convenience; any matrix norm would do).
708: The system \backdd\ can be written $$\eqalign{\vec f_1&=\bS_1(\vec y-\vec
709: f_2)\cr \vec f_2&=\bS_2(\vec y-\vec f_1).\cr}\eqn{\back}$$ Let $\vec
710: f_1^{\,(m)}$ and $\vec f_2^{\,(m)}$ denote the estimates at the $m$th stage of the
711: \index{backfitting algorithm}%
712: \index{starting functions}%
713: backfitting algorithm, with $m=0$ denoting the starting functions.
714: Backfitting consists of alternating the steps
715: $$\eqalign{\vec
716: f_1^{\,(m)}&=\bS_1(\vec y-\vec f_2^{\,(m-1)})\cr \vec f_2^{\,(m)}&=\bS_2(\vec y-\vec
717: f_1^{\,(m)}).\cr}\eqn{\backit}$$ Using induction one shows that for $m\ge 1$
718: $$\eqalign{\vec f_1^{\,(m)}&=\vec y-\sum_{j=0}^{m-1}(\bS_1\bS_2)^j(\bI-\bS_1)\vec
719: y-(\bS_1\bS_2)^{m-1}\bS_1\vec f _2^{\, (0)},\cr \vec
720: f_2^{\,(m)}&=\bS_2\sum_{j=0}^{m-1}(\bS_1\bS_2)^j(\bI-\bS_1)\vec y
721: +\bS_2(\bS_1\bS_2)^{m-1}\bS_1 \vec f_2^{\, (0)}. \cr}\eqn{\qbac}$$
722:
723: Then a sufficient condition for $\vec
724: f_1^{\,(m)}$ and $\vec f_2^{\,(m)}$ to converge is $\norm{\bS_1\bS_2} < 1$. If this is
725: the case, we can solve $\qbac$ to obtain
726: $$\eqalign{ \vec
727: f_1^{\,(\infty)}&=\{\bI-(\bI-\bS_1\bS_2)^{-1}(\bI-\bS_1)\}\vec y\cr \vec
728: f_2^{\,(\infty)}&=\bS_2(\bI-\bS_1\bS_2)^{-1}(\bI-\bS_1) \vec y\cr
729: &=\{\bI-(\bI-\bS_2\bS_1)^{-1}(\bI-\bS_2)\}\vec y.\cr }\eqn{\qconv}$$
730: The fit
731: $\hat{\vec y}$ is given by
732: $$\eqalign{\hat{\vec y}&=\vec
733: f_1^{\,(\infty)}+\vec f_2^{\,(\infty)}\cr
734: &=\{\bI-(\bI-\bS_2)(\bI-\bS_1\bS_2)^{-1}(\bI-\bS_1)\}\vec y\cr}\eqn{\qconvsum}$$ which
735: is symmetric in $\bS_1$ and $\bS_2$, as some simple calculations show.
736:
737: \index{estimating equations}%
738: \index{unique solution}%
739: This proves that if
740: $\norm{\bS_1\bS_2}<1$, the estimating equations are consistent, and have
741: a unique solution.
742: In addition, the final iterates from the
743: backfitting procedure are independent of the starting
744: values and starting order.
745:
746: Is $\norm{\bS_1\bS_2}<1$ typically?
747: If the smoothers are not centered, we
748: would have $\bS_1\bS_2\bf 1=\bf 1$ so that $\norm{\bS_1\bS_2}=1$,
749: but the centering makes $\bS_1\bS_2\bf 1=0$.
750: However, smoothers like the cubic spline
751: smoother have a second unit eigenvalue corresponding to the linear
752: \index{backfitting algorithm}%
753: \index{cubic smoothing spline}%
754: functions. Consider a backfitting algorithm with two covariates
755: $\vec x_1$ and $\vec x_2$ using cubic smoothing splines.
756: If the data show strict collinearity through the origin (i.e. $\vec
757: x_2=c\vec x_1$), we still have $\norm{\bS_1\bS_2}=1$.
758: \index{smoothing spline}%
759: With higher order splines, for example, quintic smoothing splines which also have unit quadratic
760: \index{ concurvity}%
761: eigenvectors, similar situations involving linear and quadratic functions are possible.
762: The condition $\norm{\bS_1\bS_2}=1$ is an example of {\em concurvity}, a phenomenon that we study later in
763: \index{estimating equations}%
764: the general $p$-covariate case.
765:
766: As it turns out, if $\bS_1$ and $\bS_2$ are symmetric
767: with eigenvalues in $(-1,1]$, one can prove that
768: \index{backfitting algorithm}%
769: \index{fitted value}%
770: \index{starting functions}%
771: the estimating equations are consistent.
772: Furthermore, the backfitting algorithm converges despite the presence of concurvity,
773: and
774: the fitted values
775: $\vec
776: f_1^{\,(\infty)}+\vec f_2^{\,(\infty)}$
777: are independent of the starting functions.
778: Concurvity will, however, lead to a dependence of the limits $\vec
779: f_1^{\,(\infty)}$ and $\vec f_2^{\,(\infty)}$ on the starting guess $\vec
780: f_2^{\, (0)}$.
781: That is $\vec f_+^{\,(\infty)}=
782: \vec f_1^{\,(\infty)}+\vec f_2^{\,(\infty)}$ is unique but $\vec
783: f_1^{\,(\infty)}$ and $\vec f_2^{\,(\infty)}$ are not.
784: \index{cubic smoothing spline}%
785: \index{smoother matrix}%
786: Cubic smoothing splines satisfy these conditions; however, the conditions are sufficient and not necessary.
787: Empirical evidence suggests that the results may also hold for smoothers such as locally-weighted lines, for which the smoother matrix is asymmetric and has modulus greater than one.
788:
789: We can usefully view these results as an extension of those for linear regression with a singular regression-matrix $\bX$: the fit $\hatfat{\mu}=\bX\hatfat\beta$ is unique, but $\hatfat\beta$ is not.
790:
791: \index{backfitting solutions}%
792: \index{starting functions}%
793: \index{backfitting solutions}%
794: It is not surprising that the condition $\norm{\bS_1\bS_2}=1$ leads to
795: a dependence of the backfitting solutions on the starting functions,
796: for it is immediate in this case that the backfitting solutions themselves
797: are not unique.
798: Since each $\bS_j$ is assumed to have eigenvalues in $(-1,1]$,
799: $\norm{\bS_1\bS_2}$ can only equal 1
800: if there is some vector $\vec a$ that is
801: reproduced by both $\bS_1$ and $\bS_2$.
802: But this will lead to nonuniqueness of the solution to
803: backfitting, for if $\hat\vec f_1$ and $\hat\vec f_2$ are solutions, then
804: so are $\hat\vec f_1+\vec a$ and $\hat\vec f_2-\vec a$,
805: evident upon examination of $\back$.
806:
807: \index{estimating equations}%
808: \sectionskip\section{Existence and uniqueness: $p$-smoothers}
809: The two-smoother problem has revealed important aspects of the estimating equations
810: and backfitting which can be used to motivate the more general results that
811: we discuss here.
812: In fact, the results in the general case are qualitatively the
813: same as in the simpler setting.
814: \index{symmetric smoother matrices}%
815:
816: In the two-smoother problem, we are forced to restrict attention to
817: symmetric smoother matrices with eigenvalues in $(-1,1]$.
818: The exclusion of $-1$ as an eigenvalue is necessary to avoid oscillatory behaviour in
819: the algorithm.
820: In the $p$-smoother case, it turns out that we need to assume that
821: each smoother is symmetric and has
822: eigenvalues in $[0,1]$.
823: We can give no intuitive reason for this stronger condition being
824: necessary here, but note that this is not a
825: practical limitation of the results, because any
826: reasonable symmetric smoother should satisfy this property.
827:
828: \index{estimating equations}%
829: The first issue to be settled is the existence of at least one solution
830: to the estimating equations
831: \index{estimating equations}%
832: $\backdd$. It turns out that if $\bS_1, \ldots,\bS_p$ are symmetric
833: with eigenvalues in $[0,1]$,
834: the estimating equations $\backdd$ have at least one solution for every $\vec y$.
835: Given that at least one solution exists, is it unique?
836: Nonuniqueness occurs in
837: the
838: two smoother case when $\norm{\bS_1 \bS_2}=1$.
839: \index{estimating equations}%
840: This leads us to ask: what interrelation among the $\bS_j$s in the $p$-smoother case will lead to nonuniqueness of the solution
841: to backfitting?
842: The answer lies
843: in the estimating equations $\hat\vec P \vec f=\hat\vec Q\vec y$.
844: Suppose that there is a $\vec g$ such that $\hat\vec P\vec g=\bf 0
845: $.
846: Then the system $\hat\vec P\vec f=\hat\vec Q\vec y$ has an infinite number of solutions because if
847: $\{\,\vec f_j^{\,(\infinity)}:j=1,\ldots,p\,\}$ is a solution, then
848: so is
849: $\{\,\vec f_j^{\,(\infinity)} +c\vec g_j:j=1,\ldots,p\,\}$ for any $c$.
850:
851: \index{ concurvity space}%
852: We think of this phenomenon as the analogue of collinearity, and define
853: the {\em concurvity space} of the system \backdd\
854: to be the space of functions $\vec g$
855: satisfying
856: $\hat{\vec P}\vec g = \bf 0$. Concurvity in function space is similarly
857: defined, with regard to the system \backd.
858:
859: A quick word on notation. An unsubscripted vector $\vec g$ denotes the $np$ vector of evaluations $\vec g^T=(\vec g_1^T, \ldots,\vec g_p^T)$, while an additive fit is denoted by $\vec f_+=\sum_{j=1}^p\vec f_j$.
860:
861:
862: \index{concurvity space}%
863: In the two variable case it is easy to show that
864: $\norm{\bS_1\bS_2}<1 $ and $\norm{\bS_2\bS_1}<1$ if and only the concurvity
865: space is empty.
866: Thus we see that concurvity plays an important role in the behaviour
867: of backfitting.
868: It is natural, then, to try to pin down exactly how concurvity can
869: occur.
870: We can do this, after a bit of preparation.
871: Let $\bS_j$, \ $j=1,\ldots,p$ be symmetric smoother
872: matrices with eigenvalues in $[0,1]$. Let $\scM_1(\bS_j)$
873: be the space spanned by the
874: eigenvectors of $\bS_j$ with eigenvalue +1 (that is, they pass through the smoother
875: unchanged), for $j=1,\ldots,p$.
876: Then $\hat \vec P
877: \vec g=\bf 0$ if and only if $\vec g_j\in \scM_1(\bS_j)\;\forall j$ and $\vec
878: g_+=\bf 0$.
879:
880: In other words, we have concurvity if and only if the spaces
881: $\scM_j(\bS_j)$ are linearly dependent; that is, there exist
882: $\vec g_j\in \scM_1(\bS_j)$ not all zero satisfying $\vec
883: g_+=\bf 0$. Given such a linear degeneracy, any solution $\vec
884: f_1,\ldots,\vec f_p$ of $\hat{\vec P}\vec f=\hat{\vec Q}\vec y$ leads
885: to nonuniqueness in the form of additional solutions $\vec f_1+c\vec
886: g_1,\ldots,\vec f_p+c\vec g_p$.
887:
888: The result above says that concurvity involves only functions in the eigenspaces
889: corresponding to eigenvalue +1.
890: \index{cubic smoothing spline}%
891: \index{exact concurvity}%
892: In the case of cubic smoothing splines, those eigenspaces correspond to
893: linear functions of each predictor,
894: and thus exact concurvity only exists if the predictors are exactly collinear.
895: However, approximate concurvity is of practical concern, when the predictors are clustered around some lower dimensional manifold.
896: Note that if quintic splines or
897: quadratic regression are used, the eigenspaces $\scM_1(\bS_j)$ consist of
898: the quadratic functions in the $j$th variable; hence concurvity may
899: involve truly nonlinear degeneracies between the variables.
900:
901: \sectionskip\section{Convergence of backfitting: $p$-smoothers}
902: With this definition of concurvity in hand, we can state the main result for the
903: convergence of backfitting.
904: \index{backfitting algorithm}%
905: Consider a backfitting algorithm with
906: symmetric smoothers
907: $\bS_j$, \ $j=1,\ldots,p$,
908: having eigenvalues in $[0,1]$.
909: \index{backfitting converges}%
910: \index{unique solution}%
911: \index{starting functions}%
912: Then if the $\bS_j$ do not exhibit concurvity, it can be shown that backfitting converges to
913: the unique solution of $\backdd$, independent of the starting functions.
914: \index{backfitting converges}%
915: \index{starting functions}%
916: If there is concurvity, backfitting converges to one of the
917: solutions of $\backdd$, the starting functions determining the final
918: solutions.
919:
920: \index{smoothing spline}%
921: Note that these results apply to,
922: amongst others, smoothing splines,
923: regression splines,
924: and simple linear and polynomial
925: regression.
926: \index{backfitting algorithm}%
927: \index{cubic smoothing spline}%
928: They also can be applied to a backfitting algorithm
929: that uses a mixture of these smoothers, for example a cubic smoothing spline
930: for one variable, a simple linear fit for another variable,
931: etc.
932: \index{surface smoother}%
933: The smoothers need not even be univariate; the results apply to
934: two or higher-dimensional surface smoothers as well.
935:
936: \index{estimating equations}%
937: \index{unique solution}%
938: \index{backfitting algorithm}%
939: \index{unique solution}%
940: \index{estimating equations}%
941: \index{backfitting algorithm}%
942: \sectionskip\section{Summary of the main results of the section}
943: For two smoothers $\bS_1$ and $\bS_2$:
944: \smallskip
945: {\parindent 20pt
946: \item{(i)} if $\norm{\bS_1\bS_2}<1$, then
947: the estimating equations \back\ have a unique solution \qconv\ and
948: the backfitting algorithm converges to this unique solution;
949: \item{(ii)} if $\bS_1$ and $\bS_2$ are symmetric
950: with eigenvalues in $(-1,1]$,
951: then
952: the estimating equations \back\ have at least one solution,
953: and the backfitting algorithm converges to one of the solutions.
954: This solution is dependent on the starting function
955: ${\vec f}_2^{\, (0)}$.
956:
957: \index{estimating equations}%
958: }\smallskip
959: In general for $p$ symmetric smoothers $\bS_1,\ldots,\bS_p$ with eigenvalues in $[0,1]$:
960: \smallskip
961: {\parindent 20pt
962: \item{(i)} The estimating equations $\backdd$ have at least one solution for every $\vec y$.
963: \item{(ii)} Let $\scM_1(\bS_j)$
964: be the space spanned by the
965: eigenvectors of $\bS_j$ with eigenvalue +1 (that is, they pass through the smoother
966: unchanged), for $j=1,\ldots, p$.
967: Then $\hat \vec P
968: \vec g=\bf 0$ if and only if $\vec g_j\in \scM_1(\bS_j)\;\forall j$ and $\vec
969: g_+=\bf 0$. Either of these conditions characterize the {\em concurvity} space of the
970: \index{ concurvity}%
971: \index{estimating equations}%
972: estimating equations.
973: \index{concurvity space}%
974: \index{backfitting converges}%
975: \index{unique solution}%
976: \index{starting functions}%
977: \item{(iii)} If the concurvity space is empty, backfitting converges to
978: the unique solution of $\backdd$, independent of the starting functions.
979: \index{concurvity space}%
980: \index{backfitting converges}%
981: \index{starting functions}%
982: \item{(iv)}
983: If the concurvity space is not empty, backfitting converges to one of the
984: solutions of $\backdd$, and the starting functions determine the final
985: solutions.
986:
987: \index{weighted additive model}%
988: \index{penalized least-squares criterion}%
989: \index{smoothing parameter}%
990: }
991: \Sectionskip\Section{Special topics}
992: \section{Weighted additive models}
993: Consider a weighted penalized least-squares criterion of the form
994: $$\Bigl(\vec y-\sum_j\vec f_j\Bigr)^T\bW\bigl(\vec y-\sum_j\vec f_j\Bigr)+\sum_j
995: \lambda_j \vec f_j^T \bK_j \vec f_j \eqn{\wsplinpen}$$
996: where $\bW$ is a diagonal
997: matrix of weights,
998: and $\lambda_j$ is a smoothing parameter and $\bK_j$ is a smoothing-spline penalty matrix for the $j$th predictor.
999: These weights might represent the relative precision of each observation or
1000: might arise as part of another iterative procedure,
1001: for example the local-scoring procedure
1002: described in Chapters~4 and 6.
1003: \index{estimating equations}%
1004: \index{smoothing spline}%
1005: The estimating equations for this problem have the same form as for the unweighted case, except that
1006: the smoothers are
1007: weighted smoothing splines given by
1008: $\bS_j=(\bW+\lambda_j\bK_j)^{-1}\bW$.
1009: We could generalize all the results presented so far to deal with the weighted case by simply computing norms and inner products in the metric of $\bW$.
1010: However, it is simpler
1011: to map the problem back to the unweighted case,
1012: using the transformations
1013: $\by'=\bW^{1/2}\by$,
1014: $\vec f_j'=\bW^{1/2}\vec f_j$,
1015: $\bK_j'=\bW^{-1/2}\bK_j\bW^{-1/2}$.
1016: Note that $\bS_j$ is not symmetric, but
1017: $\bW^{1/2}\bS\bW^{-1/2}$ is symmetric with eigenvalues in
1018: $[0,1]$, and unit eigenvalues
1019: corresponding to linear functions of the $j$th variable.
1020: Thus the convergence results for the unweighted case can be directly applied.
1021:
1022:
1023: \index{modified backfitting algorithm}%
1024: \index{backfitting algorithm}%
1025: \sectionskip\section{A modified backfitting algorithm}
1026: In the previous chapter we mention the possibility of modifying the backfitting
1027: algorithm to improve its efficiency.
1028: The basic idea is as follows.
1029: \index{ shrinking}%
1030: Many smoothers have a {\em projection} part and a {\em shrinking} part.
1031: \index{cubic smoothing spline}%
1032: For example, a cubic smoothing spline has unit eigenvalues that are constant
1033: and linear functions of the predictor (its projection part), and
1034: eigenvalues less than one for other eigenvectors.
1035: The idea is to combine all of the projection operations for
1036: all of the predictors into one large projection, and use only
1037: the nonprojection parts of each smoother in an iterative backfitting-type
1038: operation.
1039:
1040: This modification has several advantages.
1041: When a smoothing-spline or running-line smoother is used for
1042: several predictors,
1043: practical experience
1044: has shown that
1045: if the predictors are correlated,
1046: many iterations may be required to get the correct
1047: average slope of the functions.
1048: By performing all of the projections in one operation,
1049: all of the function slopes are simultaneously estimated.
1050: A second advantage is in collinearity/concurvity situations.
1051: \index{ concurvity}%
1052: In a backfitting problem with symmetric smoothers having
1053: eigenvalues in $[0,1]$, we have seen in the previous section that the
1054: only nonuniqueness ({\em concurvity}) occurs in the eigenspaces
1055: with eigenvalue one.
1056: By separating out the estimation of these components of the functions,
1057: the nonunique part of the solutions is
1058: conveniently allocated to the projection step.
1059: \index{starting functions}%
1060: This makes it easy to characterize the solutions and eliminates the
1061: dependence of the final solutions on the starting functions.
1062:
1063: Let us now be more specific.
1064: Let $\bG_j$ be the matrix that projects onto $\scM_1(\bS_j)$, the space of eigenvalue
1065: one for the $j$th smoother.
1066: Using $\bG_j$, we define the modified smoother matrices
1067: $$\tilde \bS_j=(\bI-\bG_j)\bS_j.$$
1068: Note that $\tilde \bS_j$ has the effect of subtracting out the component
1069: of the smoothed value that lies in $\scM_1(\bS_j)$.
1070: \index{modified backfitting algorithm}%
1071: The general form of the modified backfitting algorithm is given below.
1072: \setbox2=\vbox{\hsize \algwidth {\setnine\parindent 20pt
1073:
1074: \item{(i)} Initialize $\tilde{\vec f}_1, \ldots,\tilde{\vec f}_p$ and set $\tilde{\vec
1075: f}_+=\tilde{\vec f}_1+\cdots +\tilde{\vec f}_p$.
1076: \index{orthogonal projection}%
1077: \item{(ii)} Regress $\vec y-\tilde{\vec f}_+$ onto the space $\scM_1(\bS_1)+\cdots + \scM_1(\bS_p)$,
1078: that is,
1079: set $\vec g = \bG(\vec y-\tilde{\vec f}_+)$, where $\bG$ is the orthogonal projection
1080: onto $\scM_1(\bS_1)+\cdots +\scM_1(\bS_p)$ in $\R n$.
1081: \item{(iii)} Apply one cycle of backfitting to $\vec y-\vec g$
1082: using smoothers $\tilde{\bS}_i$; this step yields an updated
1083: additive fit $\tilde{\vec f}_+ = \tilde{\vec f} _1 + \cdots
1084: +\tilde{\vec f} _p$.
1085: \item {(iv)} Repeat steps (ii) and (iii) until
1086: convergence. The final estimate for the overall fit is $\vec f_+ =\vec g + \tilde{\vec f}_+ $.
1087:
1088: \index{modified backfitting algorithm}%
1089: }%end algorithm
1090: \smallskip
1091: } %end box 2
1092:
1093: \midinsert
1094: \algorithm{{\ninerm\noindent Algorithm \chapnodot 2} The modified backfitting algorithm}{\box2}
1095: \endinsert
1096:
1097: Note that it is not sufficient
1098: to perform the projection step only once.
1099: It must be iterated with the
1100: other steps because
1101: when $\tilde{\vec f}_+$ is changed in step (ii), the projection component
1102: $\vec g$ no longer equals $\bG(\vec y-\tilde{\vec f}_+)$.
1103: An alternative to step (iii) is to iterate it to convergence rather than cycling through once; we find that this tends to slow down convergence in terms of the number of smooths performed.
1104:
1105: In order justify this procedure, we must not
1106: only show that it does converge, but that it converges to the same solution as the
1107: original backfitting procedure.
1108: It turns out (to follow) that in the case of symmetric
1109: smoothers with eigenvalues in $[0,1]$, the modified backfitting
1110: procedure does solve the original problem.
1111: For other linear smoothers, the modified backfitting procedure might
1112: still make sense, but it solves a slightly different problem.
1113:
1114: \index{modified backfitting algorithm}%
1115: \index{modified backfitting algorithm}%
1116: We now state some convergence results about modified backfitting algorithms:
1117: \smallskip
1118: {\parindent 20pt
1119: \item{(i)}If $\bS_j$, \ $j=1,\ldots,p,$ are symmetric and have
1120: eigenvalues in $[0,1]$,
1121: then the modified
1122: backfitting algorithm converges in the sense that $\vec g$ and $ \tilde{\vec f} _1,
1123: \ldots, \tilde{\vec f} _p$ converge.
1124: \index{modified backfitting algorithm}%
1125: \item{(ii)}Suppose the modified
1126: backfitting algorithm has converged with smoothers $\tilde{\bS}_j$,
1127: \index{estimating equations}%
1128: yielding functions $\tilde{\vec f}_j$ and $\vec g_j\in \scM_1(\bS_j)$. Then the
1129: components $\vec f_j=\vec g_j + \tilde{\vec f}_j$ are solutions to
1130: the estimating equations with smoothers
1131: $\bS^*_j=\bG_j+(\bI-\bG_j)\bS_j$.
1132:
1133: \index{estimating equations}%
1134: }\smallskip
1135:
1136: Notice that if $\bS_j$ is symmetric, we have
1137: $\bS_j^*=\bS_j$ and thus the solutions to the modified algorithm solve the
1138: estimating equations with smoothers $\bS_j$.
1139:
1140: If the $\bS_j$ are symmetric and have eigenvalues in $[0,1]$ then $\tilde{\bS}_j=\bS_j-\bG_j$,
1141: \index{smoothing spline}%
1142: and $\norm{\tilde{\bS}_j }<1$. Smoothing splines belong to this
1143: \index{cubic smoothing spline}%
1144: \index{hat matrix}%
1145: class, and hence the algorithm always converges for them.
1146: If cubic smoothing splines are used for all predictors, $\bG$ is the {\em hat} matrix
1147: corresponding to the least-squares regression on $({\bf 1},{\vec x}_1,\ldots,{\vec x}_p)$.
1148: The nonlinear functions $\tilde{\vec f_j}$ are uniquely determined.
1149: Concurvity (collinearity) can show up only in the $\bG$ step, where it is dealt with in the standard linear least-squares fashion.
1150: At convergence, one may then decompose $\vec g=\sum \vec g_j$ and reconstruct final components $\vec f_j=\vec
1151: g_j+\tilde{\vec f_j}$.
1152: \index{cubic smoothing spline}%
1153: \index{multiple linear regression}%
1154: If $\bS_j$ is a cubic smoothing spline and if $\by$
1155: is centered initially, then $\vec g_j=\hat{\beta}_j\cdot\vec x_j$,
1156: where $\hat\beta_1,\ldots,\hat\beta_p$ are the coefficients from the
1157: multiple linear regression of $\vec y-\tilde{\vec f}_+$ on $\vec
1158: x_1,\ldots,\vec x_p$.
1159:
1160: \index{estimating equations}%
1161: \index{estimating equations}%
1162: \sectionskip\section{Explicit solutions to the estimating equations}
1163: By manipulating the fixed points of the modified backfitting procedures,
1164: an expression for the solutions to the estimating equations \backdd\ can be derived
1165: (Exercise~5.6).
1166: Let $\tilde \bA_j=(\bI-\tilde \bS_j)^{-1}\tilde \bS_j$, $\tilde \bA=\sum_1^p \tilde \bA_j$,
1167: and $\bB=(\bI+\tilde \bA)^{-1}\tilde \bA$.
1168: Then the solutions are $\tilde{\vec f}_+=(\bI-\bB\bG)^{-1}\bB(\bI-\bG)\vec y$ and
1169: $\vec g=\bG(\vec y-\tilde{\vec f}_+)$.
1170: These can be combined to obtain
1171: $$\eqalign{
1172: \vec f_+=&\{\bG+(\bI-\bG)(\bI-\bB\bG)^{-1}\bB(\bI-\bG)\}\vec y\cr
1173: \tilde{\vec f}_j=&(\bI-\tilde \bS_j)^{-1}(\vec y-\vec g-\tilde{\vec f}_+)\cr
1174: \vec f_j=&\vec g_j+\tilde{\vec f}_j\cr}\eqn{\mbsol}$$
1175: and the individual $\vec g_j$s are any
1176: vectors $\vec g_j\in \scM_1(\bS_j)$ such that $\sum_1^p \vec g_j=\vec g$.
1177: Interestingly, $\mbsol$ reveals that for symmetric
1178: smoother matrices with eigenvalues in $[0,1]$, a direct solution can be obtained in $O(n^3 p)$ operations,
1179: the number required for computing
1180: $(\bI-\tilde \bS_j)^{-1}$ for $j=1,\ldots, p$.
1181: \index{estimating equations}%
1182: This is less than the $O\{(np)^3\}$ operations that are needed to solve the
1183: estimating equations $\backdd$ in general.
1184:
1185: \sectionskip\section{Standard errors}
1186: From the previous sections we note that each estimated function
1187: in the additive fit is the result of a
1188: \index{additive model}%
1189: linear mapping or smoother applied to $\vec y$. This
1190: means
1191: that the variance formula developed in Chapter~3
1192: can be applied to the additive model. At convergence, we can
1193: express $\hatvec f_j$ as
1194: $\bR_j\vec y$ for some $n\times n$ matrix $\bR_j$.
1195: If the observations have independent and identically distributed errors,
1196: then $\cov(\hatvec{f}_j)=\bR_j\bR_j^T\sigma^2$ where
1197: $\sigma^2=\var(Y_i)$. As in the least-squares case, if $\hat{\vec P}$ in equation \backdd\
1198: has singular values close to $\bf 0$, this will be reflected in $\cov(\hatvec
1199: f\,)$ as large variances and covariances.
1200:
1201: \index{semi-parametric model}%
1202: Direct computation of $\bR_j$ is
1203: formidable, except in very special cases such as the semi-parametric model.
1204: Our best general approach to date is to apply the backfitting procedure to the each of the $n$ unit
1205: $n$-vectors that are the columns of $\bI_n$, the $n\times n$ identity matrix.
1206: The result of backfitting applied to the $i$th unit vector produces fitted vectors $\hatvec{f}_j^{\,i}$, $j=1,\ldots,p$, where $\hatvec{f}_j^{\,i}$ is the $i$th column
1207: of $\bR_j$.
1208: Similarly, $\hatvec{f}_+^{\,i}$ is the $i$th column of $\bR$.
1209: \index{standard-error bands}%
1210: \index{diagonal elements}%
1211: The standard-error bands in Fig.~\fone\ of Chapter~4 are constructed
1212: \index{backfitting algorithm}%
1213: using $\pm$ twice the square root of the diagonal elements of $\hat{\sigma}^2\bR_j\bR_j^T$.
1214: Since the backfitting algorithm is $O(kn)$ for $O(n)$ smoothers, this procedure is $O(kn^2)$.
1215: Now $k=pmC$, where $p$ is the number of predictors, $m$ is the number of backfitting iterations, and $C$ is the constant for the particular smoother.
1216: \index{smoothing spline}%
1217: For smoothing splines, typical numbers might be $k=5\times 5\times 35= 875$, which is likely to be larger than $n$, so this task can be tedious in practice.
1218: We have nevertheless used it in many examples, although usually not often
1219: within any single analysis.
1220:
1221: \index{additive model}%
1222: The global confidence set techniques described in Chapter~3 can be extended to apply to additive models;
1223: the procedure is similar to the univariate case.
1224: \index{additive model}%
1225: Suppose the additive model is correct, i.e., $Y_i=f_+(\vec X_i)+\varepsilon_i$,
1226: and our estimate of $\bR\vec f_+=\vec g_+ $ is $\hatvec{ f}_+=\bR\vec y$.
1227: Then an approximate pivotal for $\vec g_+$ is
1228: $$\nu(\vec g_+)=(\hatvec f_+-\bg_+)^T(\bR\bR^T\hat\sigma^2)^{-1}
1229: (\hatvec f_+-\bg_+).\eqn{\nuu}
1230: $$
1231: Assuming that we have an estimate of the dispersion parameter $\hat\sigma^2$ and the distribution $G$ of $\nu$, then we can construct a simultaneous $1-\alpha$ confidence set of all the component functions:
1232: $$C(\bg_1,\ldots,\bg_p)=\{\bg_1,\ldots,\bg_p;\nu(\bg_+) \leq G_{1-\alpha}\}.\eqn{\bogplus}$$
1233: We will not pursue this topic further here; it is an area of current research and we need to gain more experience with it.
1234:
1235: \sectionskip\section{Degrees of freedom}
1236: Each of the definitions for degrees of freedom given in Chapter~3 has a natural analogue here.
1237: The overall degrees of freedom $\df$ is simply $\tr(\bR)$,
1238: where $\bR$ is the (smoother) matrix that produces $\hatvec f_+=\bR \vec y$.
1239: \index{posterior covariance}%
1240: \index{additive model}%
1241: In addition, the posterior covariance of $\vec f_+$,
1242: \index{prior covariance}%
1243: in the Bayesian treatment of the additive model given in section~5.4.6,
1244: is proportional to $\bR$, for appropriate choice of
1245: the prior covariances.
1246:
1247: Similarly, the degrees of freedom for error is
1248: $\dferr=n-\tr(2\bR-\bR\bR^T)$.
1249: More usefully,
1250: for model comparison, we need a notion of the change in the error degrees of freedom $\Delta\dferr$ due to an individual term.
1251: Let $\bR_{(j)}$ denote that operator that produces the additive fit with
1252: the $j$th term removed.
1253: Then we define $\dferr_j$, the degrees of freedom for error due to the $j$th term:
1254: $$\dferr_j=\trace(2\bR-\bR\bR^T)-\trace(2\bR_{(j)}-\bR_{(j)}\bR_{(j)}^T).$$
1255: This is the expected increase in the residual sum of squares (up to a scale factor) if the $j$th predictor is excluded from the model, assuming its exclusion does not increase the bias.
1256: Approximate $F$ and $\chi^2$ tests that make use of $\dferr_j$ are discussed
1257: in section~6.8.
1258:
1259: \index{fitted value}%
1260: The sum of the variances of the fitted values is a meaningful concept
1261: for an additive fit and thus
1262: $\dfvar=\trace(\bR\bR^T)$.
1263: Further, the sum of the variances of the fitted component function $\hat\vec f_j$ is
1264: $\sigma^2\trace(\bR_j\bR_j^T)$;
1265: the effect of predictor correlation on this quantity is
1266: explored in Exercise~5.17.
1267:
1268:
1269:
1270:
1271:
1272: None of these definitions are attractive from a computational point of view.
1273: In particular, it would be convenient to use $\tr(\bS_j)-1$ or
1274: even $\tr(2\bS_j-\bS_j\bS_j^T)-1$ to select the amount of smoothing
1275: prior to including the $j$th predictor in a model,
1276: and as an approximation to $\dferr_j$ for model comparison.
1277: We subtract one since there is a redundant constant in $p-1$ of the $p$ terms
1278: in the model; in general we subtract the dimension of $\bigcap_j\scM_1(\bS_j)$.
1279: \index{exact concurvity}%
1280: In the extreme case of exact concurvity, it is possible to show that $\tr(2\bS_j-\bS_j\bS_j^T)-1$
1281: \index{balanced additive}%
1282: is an upper bound for $\dferr_j$ (Exercise~5.7);
1283: for a balanced additive model (the other extreme; Exercise~5.18) it is equal to $\dferr_j$.
1284: ^{Buja, Hastie and Tibshirani (1989)} carried out some small simulation experiments
1285: and found that adding up the
1286: individual degrees of freedom gave a good approximation to the
1287: true degrees of freedom.
1288: \index{smoothing parameter}%
1289: The only exceptions occurred when the predictors had extremely high
1290: correlation or when a very small smoothing parameter was used.
1291:
1292: In the examples in this book, we use the convenient
1293: approximation $\tr(\bS_j)-1$ to select
1294: the amount of smoothing, while we use the exact quantities $\dferr_j$
1295: for model comparisons via approximate $F$ or $\chi^2$ tests.
1296: Further details are given in section~6.8.
1297:
1298: \index{additive model}%
1299: \index{additive model}%
1300: \sectionskip\section{A Bayesian version of additive models}
1301: Just as in the case of a single smoother, there is a rather simple Bayesian
1302: approach to additive models.
1303: As in section~3.6, there is a functional (stochastic process) version,
1304: and a finite dimensional (sampled) version; we focus on the latter.
1305:
1306: The model is
1307: $$
1308: \vec y = \vec f_1+\cdots +\vec f_p+\fat{\varepsilon}\eqn{\bayesadd}$$
1309: with $\vec f_j \sim N({\bf 0},\sigma^2 \bQ_j)$ independently for all $j$ and independent of ${\fat\varepsilon}\sim N({\bf 0},\sigma^2 \bI)$.
1310: Note from section~3.6 that each $\bQ_j$ corresponds to the inverse of some penalty matrix $\bK_j$; also, from section~5.2.3, $\bQ_j$ can be identified with the realization of a reproducing kernel. For the moment we assume that the priors are proper and nonsingular.
1311:
1312: Straightforward derivations (Exercise~5.9) show the following:
1313: \smallskip
1314: {\parindent 20pt
1315: \item{(i)} The prior for $\vec f_+=\vec f_1+\cdots +\vec f_p$
1316: \index{posterior mean}%
1317: is $N({\bf 0},\bQ_+)$ with $\bQ_+=\sum_j\bQ_j$.
1318: \index{posterior mean}%
1319: \item{(ii)} The posterior mean for $\vec f_+$ is $\ev(\vec f_+\given \vec y)= \bQ_+(\bI+\bQ_+)^{-1}\vec y$.
1320: Similarly the posterior means for the individual functions are $\ev(\vec f_j\given \vec y)= \bQ_j(\bI+\bQ_+)^{-1}\vec y$.
1321: \index{posterior covariance}%
1322: \item{(iii)} The posterior covariance of $\vec f_+$ is $\sigma^2\bQ_+(\bI+\bQ_+)^{-1}$, and for $\vec f_j$ they are $\sigma^2(\bQ_j-\bQ_j(\bI+\bQ_+)^{-1}\bQ_j)$.
1323:
1324: }\smallskip
1325:
1326: Notice that the solution involves the inversion of an $n\times n$ {\em unstructured} matrix, which takes $O(n^3)$ operations, unless of course backfitting is used.
1327:
1328: Similar equations can be derived for partially improper priors.
1329: These give infinite variance to certain components in $\R n$, which correspond
1330: to components in $\scM_1(\bS_j)$.
1331: The simplest way to formulate the problem is along the lines taken in section~5.2.3, where the projection space is explicitly isolated.
1332: \index{smoothing parameter}%
1333: In addition, smoothing parameters are thought of as prior variances
1334: for the $\vec f_j$s; we have absorbed these into the $\bQ_j$s in the
1335: above formulation.
1336:
1337: \index{additive model}%
1338: \Sectionskip\Section{Bibliographic notes}
1339: The theory of nonparametric additive modelling is relatively recent.
1340: \index{ace algorithm}%
1341: \index{additive model}%
1342: \index{response transformation}%
1343: ^{Breiman and Friedman (1985)} proposed the ACE algorithm, a procedure
1344: more general than the additive model (allowing response transformations, and discussed in
1345: Chapter~7) and proved many results on convergence and consistency
1346: of backfitting in both Hilbert-space and data settings.
1347: Most of this chapter is based on the paper by Buja, Hastie, and Tibshirani (1989)
1348: in which some of Breiman and Friedman's convergence results were
1349: extended and the concurvity and degrees of freedom results are
1350: presented.
1351: In some cases, results stronger than those given in this chapter
1352: can be derived; see ^{Buja \etal~(1989)} and the discussions.
1353: We have traded generality for simplicity and interpretability
1354: in this chapter.
1355:
1356: Backfitting goes back a long way.
1357: ^{Friedman and Stuetzle (1981)} defined the term in the context of projection-pursuit regression,
1358: \index{time series}%
1359: while ^{Wecker and Ansley (1983)} suggested its use in the context of economic models.
1360: Some of the earlier references include ^{Papadakis (1937) } for the separation of fertility trends in the analysis of field trials, and ^{Shiskin, Young and Musgrave (1967)} in the X-11 system for decomposing time series (Chapter~8).
1361: \index{additive model}%
1362: ^{Kohn and Ansley (1989)} independently studied the properties of additive models in the stochastic setting, including convergence.
1363:
1364: Much of the statistical theory and practice of one and higher-dimensional spline models is due to Grace Wahba and her co-workers.
1365: Some relevant references are ^{ Kimeldorf and Wahba (1971)}, ^{Wahba (1978, 1980, 1986)}, ^{O'Sullivan (1983)}, ^{Gu and Wahba (1988)}, and ^{Chen, Gu and Wahba (1989)}.
1366: Section~5.2.3 is based almost entirely on this last reference.
1367: \index{additive model}%
1368: Cox (1989) ^^{Cox, D.D., 1989} described the Bayesian formulation of additive models,
1369: summarized in section~5.4.6.
1370:
1371: The notion of concurvity was introduced by ^{Buja, Donnell and Stuetzle}
1372: (1986)
1373: and was also discussed in ^{Buja \etal~(1989)}.
1374: \index{semi-parametric model}%
1375: ^{Bickel, Klaassen, Ritov, and Wellner (1990) } studied the theory
1376: for semi-parametric models,
1377: and in the two predictor projection case, showed convergence of the
1378: backfitting functions and computed the rates of convergence.
1379:
1380: The semi-parametric approach discussed in section~5.3.3 was considered
1381: by ^{Engle, Granger, Rice and Weiss (1986)}, ^{Denby (1986)}, ^{Wahba (1986)}, ^{Green and
1382: Yandell (1985)}, ^{Green (1985)},
1383: ^{Eubank (1985)},
1384: ^{Heckman (1986, 1988)},
1385: ^{Speckman (1988)}, and ^{Shiau and Wahba (1988)}.
1386: Green (1985), and Green and Yandell (1985) looked at
1387: regression and more general models
1388: with a single nonparametric term
1389: (we discuss their work in Chapter~6).
1390: ^{Heckman (1986)} proved consistency of the regression estimate in a regression
1391: model with a single cubic spline term and showed that the estimates
1392: of the regression coefficients and nonparametric function
1393: are Bayes estimates under an appropriate diffuse prior,
1394: generalizing the work of ^{Wahba (1978)}.
1395: ^{Rice (1986)} studied the convergence rates for these partially-splined models.
1396: Speckman (1988) compared the bias and variance of estimators for this model,
1397: and proposed a new estimate (Exercise~5.13) with asymptotically lower-order bias.
1398: This estimate was also suggested by ^{Denby (1986)}.
1399: ^{Shiau and Wahba (1988)} did a thorough study of bias and variance for these
1400: models, and
1401: ^{Heckman (1988) }studied minimax estimators.
1402:
1403: \index{additive model}%
1404: Stone (1982) ^^{Stone (C.J., 1982)} studied rates of convergence for additive models, with the
1405: functions estimated by polynomials or regression splines.
1406: \index{optimal rate}%
1407: \index{additive model}%
1408: He proved the interesting result
1409: that the optimal rate of convergence for an estimate
1410: of the additive model
1411: is the same as
1412: that for a single function, discussed in section~3.10.
1413: Thus an increase in the dimension $p$ does not decrease the rate of convergence,
1414: as it does if one is estimating a general (nonadditive) $p$-dimensional
1415: function.
1416:
1417: The Gauss-Seidel algorithm is discussed in most textbooks on
1418: numerical analysis; see for example ^{Golub and Van Loan (1983)}.
1419: For a description of the QR algorithm, see ^{Thisted (1988)}.
1420:
1421: \index{additive model}%
1422: \Sectionskip\Section{Further results and exercises 5}
1423: \Mark{EXERCISES \ 5}%
1424: \beginexercises
1425: \exercise
1426: Consider an additive model with
1427: symmetric but not necessarily invertible smoother matrices.
1428: \index{penalized least-squares criterion}%
1429: \index{estimating equations}%
1430: Show that the minimizers of the penalized least-squares criterion $\gsplinpen$
1431: are the solutions to the estimating equations $\backdd$.
1432:
1433: \index{orthogonal projection}%
1434: [^{Buja, Hastie and Tibshirani, 1989]}
1435: \exercise
1436: Complete the proof of convergence of backfitting in the case of
1437: orthogonal projections (section~5.3.2) by showing that
1438: $\norm{\bC^m\hat\vec y}\rightarrow 0$.
1439: (Hint: show that $\norm{\bC\vec a}\leq \norm{\vec a}$ for all $\vec a$,
1440: with equality if and only if $\vec a$ is in the orthogonal complement
1441: of $V$).
1442: \index{orthogonal projection}%
1443: \exercise
1444: Consider a backfitting procedure with orthogonal projections,
1445: and
1446: let $ \bD$ be the overall design matrix whose columns span
1447: $V=\script L_{col}(\bS_1)\oplus
1448: \script L_{col}(\bS_2)\oplus \ldots \oplus\script L_{col}(\bS_p)$.
1449: \index{estimating equations}%
1450: Show that the estimating equations $\hat {\vec P}\vec f=\hat {\bQ}\vec y$
1451: are equivalent to the least-squares normal equations
1452: $\bD^T\bD\fat \beta=\bD^T\vec y$ where $\fat \beta$ is the vector
1453: of coefficients.
1454: \exercise
1455: In a backfitting procedure with two least-squares projections $\bS_1$ and
1456: $\bS_2$ based on predictors $\vec x_1$ and $\vec x_2$, show that
1457: if $\theta\neq 0$, then
1458: the difference between the $i$th iterate and the solution converges to zero
1459: geometrically at rate
1460: $\cos{\theta}$,
1461: where
1462: $\theta$ is the angle between $\vec x_1$ and $\vec x_2$.
1463: \index{semi-parametric model}%
1464: \index{backfitting converges}%
1465:
1466: [^{Deutsch, 1983}]
1467: \exercise Prove that in the case of the semi-parametric model, backfitting converges to the solution $\simmm$.
1468: Give conditions that guarantee the solution is unique.
1469:
1470: [^{Green, 1985}]
1471: \index{estimating equations}%
1472: \exercise Prove item (ii) in section~5.4.2, stating that the modified backfitting
1473: procedure provides a solution to the estimating equations $\backdd$.
1474: \index{estimating equations}%
1475: \index{modified backfitting algorithm}%
1476: Derive the explicit solutions $\mbsol$ to the estimating equations from the fixed
1477: points of the modified backfitting algorithm.
1478:
1479: \index{additive model}%
1480: [^{Buja, Hastie and Tibshirani, 1989}]
1481: \exercise
1482: As an extreme case of concurvity, consider an additive model with $p$
1483: identical smoothers $\bS$.
1484: Assume that $\bS$ is centered and symmetric,
1485: with eigenvalues in $[0,1]$, and that the dimension of $\scM_1(\bS)$ is $q$.
1486: Show that
1487: {\parindent 20pt
1488: \item{(i)} $\trace(2\bS-\bS^2)\leq \trace(2\bR-\bR^2)\leq p\;\trace(2\bS-\bS^2)-(p-1)q$.
1489: \item{(ii)} $\trace(2\bR-\bR^2)-\trace(2\bR_{(j)}-\bR_{(j)}^2)\leq \trace(2\bS-\bS^2)-q.$
1490:
1491: }
1492: Hence conclude that
1493: {\parindent 20pt
1494: \item{(a)} including the same predictor twice or more increases the degrees of
1495: freedom (in contrast to projections such as linear regression);
1496: \item{(b)} the sum of the degrees of freedom of the individual function estimates
1497: provides an upper bound on the degrees of freedom of the fitted model;
1498: \item{(c)} $\trace(2\bS-\bS^2)-q$ provides an upper bound on the degrees of freedom contribution $\dferr_j$.
1499:
1500: \index{smoother matrix}%
1501: \index{kernel smoother}%
1502: }\smallskip
1503:
1504: [^{Buja, Hastie and Tibshirani, 1989]}
1505:
1506: \exercise Compute the smoother matrix for a running-mean and kernel smoother
1507: for a small dataset, and hence find empirically that these smoothers are
1508: not symmetric
1509: and that their eigenvalues can have moduli larger than one.
1510: \index{posterior mean}%
1511: \index{additive model}%
1512: \exercise Derive the expressions (i)--(iii) in section~5.4.6 for the Bayes posterior means and variances for a stochastic additive model.
1513: \index{smoothing spline}%
1514: \exercise Starting with equations \backdd\ with each of the smoothers a smoothing spline, derive an equivalent form for $\hatvec{f}_+$ as in the previous exercise.
1515: \index{basis function}%
1516: \index{smoothing spline}%
1517: \index{additive model}%
1518: \exercise Suppose that the $B$-spline basis functions are to be used to represent the smoothing splines in an additive model fit.
1519: \index{additive model}%
1520: \smallskip
1521: {\parindent 20pt
1522: \item{(i)}Derive closed form expressions for the $B$-spline coefficients of the additive model solution.
1523: \index{fitted function}%
1524: \item{(ii)} Derive expressions for $\hat{f}_+(\vec x^0)$ and $\hat{f}_j(x_{0j})$, the fitted functions evaluated at an arbitrary point $\vec x^0$.
1525:
1526: \index{cubic smoothing spline}%
1527: }
1528: \exercise What is the rank of $\hat {\vec P}$ in \backdd, if the smoothers are all cubic smoothing splines, and the $\vec x_j$ span a $p$ dimensional space?
1529: Identify the null space.
1530: \exercise Derive the value of $\fat\beta$ that minimizes
1531: $$\norm{\vec y-\vec X\fat{\beta}-\bS(\vec y- \vec X\fat{\beta})}\eqn{\twicing}$$
1532: \index{additive model}%
1533: and compare it with the semi-parametric estimator \simmm\ in section~5.3.3, using the same smoother $\bS$ in both cases.
1534: \exercise Suppose the same smoother $\vec S$ is used to estimate both the terms
1535: in a two-term additive model (that is, both variables are identical).
1536: \index{semi-parametric model}%
1537: Show that the backfitting residual converges to $(\bI+\bS)^{-1}(\bI-\bS)\vec y$, and that the residual sum of squares converges upwards. Can the residual sum of squares converge upwards in less structured situations?
1538: \exercise Consider a semi-parametric model with $p$ predictors (including the constant) and an additional smoother $\bS$.
1539: Derive explicit expressions for $\bR$ and $\bR_j$, and show that $\trace(\bR)\leq p+\trace(\bS)-1$.
1540: \exercise Suppose $\bS_1$ and $\bS_2$ both have the same eigenspaces; this occurs frequently in time-series applications, as in Chapter~8, where smoothing takes place in the Fourier domain. Let $\lambda_{1k}$ and $\lambda_{2k}$, $k=1,\ldots,n$,
1541: \index{constant term}%
1542: be the eigenvalue pair for the $k$th eigenvector, with $\lambda_{jk}\in [0,1]\;\forall k,j$. Assume for convenience that the constant term has been removed, and
1543: that $\lambda_{1k}\lambda_{2k}<1$.
1544: \smallskip
1545: {\parindent 20pt
1546: \item{(i)}Show that the additive-fit operators $\bR=\bR_1+\bR_2$ have eigenvalues
1547: $$\eqalign{\lambda_{k}(\bR)&=1-{(1-\lambda_{1k})(1-\lambda_{2k})\over 1-\lambda_{1k}\lambda_{2k}}\cr
1548: \lambda_{k}(\bR_1)&={\lambda_{1k}(1-\lambda_{2k})\over 1-\lambda_{1k}\lambda_{2k}}\cr
1549: \lambda_{k}(\bR_2)&={\lambda_{2k}(1-\lambda_{1k})\over 1-\lambda_{1k}\lambda_{2k}}\cr}
1550: $$
1551: \item{(ii)}Conclude from (i) that components with eigenvalue one for either smoother
1552: get totally absorbed into the corresponding term.
1553: \item{(iii)}Show that $\trace(\bR)-\trace(\bS_1)\leq\trace(\bS_2)$
1554: \item{(iv)}Show that $\trace(\bR_j\bR_j^T)\leq\trace(\bS_j\bS_j^T)$, and interpret this result in terms of pointwise variances.
1555:
1556: }
1557: \exercise Consider two extremely simple smoothers: $\bS_1=\vec u\lambda\vec u^T$ and $\bS_2=\vec v\lambda\vec v^T$, with $\lambda\in[0,1]$, and $\norm{\vec u}=\norm{\vec v}=1$.
1558: Let $\bR_j$ be the additive operators as above.
1559: Show that
1560: $$\trace(\bR_j\bR_j^T)={\lambda^2\{1-c^2\lambda(2-\lambda)\}\over(1-c^2\lambda^2)^2},$$ where $c=\langle \vec u,\vec v\rangle$.
1561: For $c=0$, $\trace(\bR_j\bR_j^T)=\lambda^2=\trace(\bS_j\bS_j^T)$, while for $c=1$, $\trace(\bR_j\bR_j^T)= \{\lambda/(1+\lambda)\}^2$.
1562: Investigate for values of $c\in (0,1)$ and a range of values of $\lambda$.
1563:
1564: \index{balanced additive}%
1565: \exercise
1566: Consider a balanced additive model as defined in Exercise~4.8.
1567: Let $\bS_0={\bf 11}^T/n$ denote the {\em mean smoother}.
1568: Using the notation of section~5.4.5,
1569: show that
1570:
1571: {\parindent 20pt
1572: \item{(i)}$\bR_j$=$\bS_j-\bS_0$, $\bR=\bS_0+\sum_{j=1}^p\bR_j$, and $\bR_j\bR_k =\bf 0$ for $j\neq k$.
1573: \item{(ii)} The full and marginal definitions of $\dferr$ and $\dfvar$ coincide, that is
1574: \index{additive model}%
1575: \index{smoothing spline}%
1576: $\tr(2\bR-\bR\bR^T)-\tr(2\bR_{(j)}-\bR_{(j)}\bR_{(j)}^T)=\tr(2\bS_j-\bS_j\bS_j^T)-1$ and $\tr(\bR_j\bR_j^T)=\tr(\bS_j\bS_j^T)-1$.
1577:
1578: }
1579: \exercise
1580: \index{estimating equations}%
1581: Suppose each of the predictors in an additive model have ties, and smoothing splines are to be used in the fit. Describe the smoother matrices $\bS_j$, as well as their eigenstructure. Show how the estimating equations \backdd\ can be reduced from the
1582: default dimension $np$ to $\sum_{j=1}^pm_j$, where $m_j$ is the number of unique values in $\vec x_j$.
1583: \endexercises
1584: \vfill\supereject
1585: \end
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.