Return to foo CVS log

Up to [Research Unix] / researchv10dc / vol2 / index

Annotation of researchv10dc/vol2/index/foo, revision 1.1

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

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.