Backfitting algorithm explained

In statistics, the backfitting algorithm is a simple iterative procedure used to fit a generalized additive model. It was introduced in 1985 by Leo Breiman and Jerome Friedman along with generalized additive models. In most cases, the backfitting algorithm is equivalent to the Gauss - Seidel method, an algorithm used for solving a certain linear system of equations.

Algorithm

Additive models are a class of non-parametric regression models of the form:

Y _{i} = \alpha + \sum ^{p} _{j =1} f _{j(X} _{j} ) + \epsilon _{i}

where each

X _{1,} X _{2,} \ldots, X _{p} is a variable in our p -dimensional predictor X , and Y is our outcome variable. \epsilon represents our inherent error, which is assumed to have mean zero. The f _{j} represent unspecified smooth functions of a single X _{j} . Given the flexibility in the f _{j} , we typically do not have a unique solution: \alpha is left unidentifiable as one can add any constants to any of the f _{j} and subtract this value from \alpha . It is common to rectify this by constraining\sum ^{N} _{i =1} f _{j(X} _{ij} ) = 0 for all j

leaving

\alpha = 1 /N \sum ^{N} _{i =1} y _{i}

necessarily.

The backfitting algorithm is then: Initialize

\hat{\alpha} = 1 /N \sum ^{N} _{i =1} y _{i,} \hat{f _{j}} \equiv 0 ,\forall j Do until \hat{f _{j}} converge: For each predictor j: (a) \hat{f _{j}} \leftarrow Smooth [\lbrace y _{i} - \hat{\alpha} - \sum _{k} \hat{f _{k}(x} _{ik} ) \rbrace ^{N} _{1} ] (backfitting step) (b) \hat{f _{j}} \leftarrow \hat{f _{j}} - 1 /N \sum ^{N} _{i =1} \hat{f _{j}(x} _{ij} ) (mean centering of estimated function)

where

Smooth is our smoothing operator. This is typically chosen to be a cubic spline smoother but can be any other appropriate fitting operation, such as:

• local polynomial regression
• kernel smoothing methods
• more complex operators, such as surface smoothers for second and higher-order interactions

In theory, step (b) in the algorithm is not needed as the function estimates are constrained to sum to zero. However, due to numerical issues this might become a problem in practice.^[1]

Motivation

If we consider the problem of minimizing the expected squared error:

min E[Y - (\alpha + \sum ^{p} _{j =1} f _{j(X} ^{2} ^{j))]}

There exists a unique solution by the theory of projections given by:

f _{i(X} _{i)} = E[Y - (\alpha + \sum ^{p} _{j ≠ i} f _{j(X} _{j))} | X _{i]}

for i = 1, 2, ..., p.

This gives the matrix interpretation:

\begin{pmatrix} I & P _{1} & … & P _{1} \\ P _{2} & I & … & P _{2} \\ \vdots & & \ddots & \vdots \\ P _{p} & … & P _{p} & I \end{pmatrix} \begin{pmatrix} f _{1(X} _{1)\\ f} _{2(X} _{2)\\ \vdots} \\ f _{p(X} _{p) \end{pmatrix} = \begin{pmatrix} P} _{1} Y\\ P _{2} Y\\ \vdots \\ P _{p} Y \end{pmatrix}

where

P _{i( ⋅ )} = E( ⋅ |X _{i)} . In this context we can imagine a smoother matrix, S _{i} , which approximates our P _{i} and gives an estimate, S _{i} Y , of E(Y|X) \begin{pmatrix} I & S _{1} & … & S _{1} \\ S _{2} & I & … & S _{2} \\ \vdots & & \ddots & \vdots \\ S _{p} & … & S _{p} & I \end{pmatrix} \begin{pmatrix} f _{1\\ f} _{2\\ \vdots} \\ f _{p \end{pmatrix} = \begin{pmatrix} S} _{1} Y\\ S _{2} Y\\ \vdots \\ S _{p} Y \end{pmatrix}

or in abbreviated form

\hat{S}f = QY

An exact solution of this is infeasible to calculate for large np, so the iterative technique of backfitting is used. We take initial guesses

f ^{( 0 )} _{j} and update each f ^{(\ell)} _{j} in turn to be the smoothed fit for the residuals of all the others:\hat{f ^{(\ell)} _{j}} \leftarrow Smooth [\lbrace y _{i} - \hat{\alpha} - \sum _{k} \hat{f _{k}(x} _{ik} ) \rbrace ^{N} _{1} ]

Looking at the abbreviated form it is easy to see the backfitting algorithm as equivalent to the Gauss - Seidel method for linear smoothing operators S.

Explicit derivation for two dimensions

Following,^[2] we can formulate the backfitting algorithm explicitly for the two dimensional case. We have:

f _{1} = S _{1(Y-f} _{2),} f _{2} = S _{2(Y-f} _{1)}

If we denote

\hat{f} ^{(i)} _{1} as the estimate of f _{1} in the ith updating step, the backfitting steps are\hat{f} ^{(i)} _{1} = S _{1[Y} - \hat{f} ^{(i -1 )} _{2} ], \hat{f} ^{(i)} _{2} = S _{2[Y} - \hat{f} ^{(i)} _{1} ]

By induction we get

\hat{f} ^{(i)} _{1} = Y - \sum ^{i -1} _{\alpha =0} (S _{1} S ^{\alpha(I - S} _{1 )Y} - (S _{1} S ^{i -1} _{2 )} S _{1\hat{f}} ^{( 0 )} ^{2}

and

\hat{f} ^{(i)} _{2} = S _{2} \sum ^{i -1} _{\alpha =0} (S _{1} S ^{\alpha(I - S} _{1 )Y} + S _{2(S} _{1} S ^{i -1} _{2 )} S _{1\hat{f}} ^{( 0 )} ^{2}

If we set

\hat{f} ^{( 0 )} _{2} = 0 then we get\hat{f} ^{(i)} _{1} = Y - S ^{-1} _{2} \hat{f} ^{(i)} _{2} = [I - \sum ^{i -1} _{\alpha =0} (S _{1} S ^{\alpha(I - S} _{1 )]Y} \hat{f} ^{(i)} _{2} = [S _{2} \sum ^{i -1} _{\alpha =0} (S _{1} S ^{\alpha(I - S} _{1 )]Y}

Where we have solved for

\hat{f} ^{(i)} _{1} by directly plugging out from f _{2} = S _{2(Y-f} _{1)} .

We have convergence if

\|S _{1} S _{2\|} < 1 . In this case, letting \hat{f} ^{(i)} _{1} , \hat{f} ^{(i)} _{2} \xrightarrow{} \hat{f} ^{( inf ty)} _{1} , \hat{f} ^{( inf ty)} _{2} :\hat{f} ^{( inf ty)} _{1} = Y - S ^{-1} _{2} \hat{f} ^{( inf ty)} _{2} = Y - (I - S _{1} S ^{-1} _{2 )} (I - S _{1)} Y \hat{f} ^{( inf ty)} _{2} = S _{2} (I - S _{1} S ^{-1} _{2 )} (I - S _{1)} Y

We can check this is a solution to the problem, i.e. that

\hat{f} ^{(i)} _{1} and \hat{f} ^{(i)} _{2} converge to f _{1} and f _{2} correspondingly, by plugging these expressions into the original equations.

Issues

The choice of when to stop the algorithm is arbitrary and it is hard to know a priori how long reaching a specific convergence threshold will take. Also, the final model depends on the order in which the predictor variables

X _{i} are fit.

As well, the solution found by the backfitting procedure is non-unique. If

b is a vector such that \hat{S}b = 0 from above, then if \hat{f} is a solution then so is \hat{f} + \alpha b is also a solution for any \alpha \in R . A modification of the backfitting algorithm involving projections onto the eigenspace of S can remedy this problem.

Modified algorithm

We can modify the backfitting algorithm to make it easier to provide a unique solution. Let

l{V} _{1(S} _{i)} be the space spanned by all the eigenvectors of S_i that correspond to eigenvalue 1. Then any b satisfying \hat{S}b = 0 has b _{i} \in l{V} _{1(S} _{i)} \forall i =1 , ... ,p and \sum ^{p} _{i =1} b _{i} = 0 . Now if we take A to be a matrix that projects orthogonally onto l{V} _{1(S} _{1)} + ... + l{V} _{1(S} _{p)} , we get the following modified backfitting algorithm:

Initialize

\hat{\alpha} = 1 /N \sum ^{N} _{1} y _{i,} \hat{f _{j}} \equiv 0 ,\forall i, j , \hat{f _{+}} = \alpha + \hat{f _{1}} + ... + \hat{f _{p}} Do until \hat{f _{j}} converge: Regress y - \hat{f _{+}} onto the space l{V} _{1(S} _{i)} + ... + l{V} _{1(S} _{p)} , setting a = A(Y - \hat{f _{+})} For each predictor j: Apply backfitting update to (Y - a) using the smoothing operator (I - A _{i)S} _{i} , yielding new estimates for \hat{f _{j}}

References

• 10.2307/2288473 . Estimating optimal transformations for multiple regression and correlations (with discussion) . 2288473 . Breiman, L. . Friedman, J. H. . amp . Journal of the American Statistical Association . 80 . 391 . 580–619 . 1985.
• Generalized Additive Models . Hastie, T. J. . Tibshirani, R. J. . amp . Monographs on Statistics and Applied Probability . 43 . 1990.
• Web site: Backfitting . Härdle, Wolfgang . June 9, 2004 . 2015-08-19 . etal . dead . https://web.archive.org/web/20150510225240/http://sfb649.wiwi.hu-berlin.de/fedc_homepage/xplore/ebooks/html/spm/spmhtmlnode37.html . 2015-05-10 .

External links

• R Package for GAM backfitting
• R Package for BRUTO backfitting

Notes and References

1. [Trevor Hastie|Hastie, Trevor]
2. Härdle, Wolfgang; et al. (June 9, 2004). "Backfitting". Archived from the original on 2015-05-10. Retrieved 2015-08-19.