Generalized additive model

It had been known since the 1950s (via. the Kolmogorov–Arnold representation theorem) that any multivariate function could be represented as sums and compositions of univariate functions.

${\displaystyle f({\vec {x}})=\sum _{q=0}^{2n}\Phi _{q}\left(\sum _{p=1}^{n}\phi _{q,p}(x_{p})\right)}$

Unfortunately, though the Kolmogorov–Arnold representation theorem asserts the existence of a function of this form, it gives no mechanism whereby one could be constructed. Certain constructive proofs exist, but they tend to require highly complicated (i.e. fractal) functions, and thus are not suitable for modeling approaches. Therefore, the generalized additive model[1] drops the outer sum, and demands instead that the function belong to a simpler class.

${\displaystyle f({\vec {x}})=\Phi \left(\sum _{p=1}^{n}\phi _{p}(x_{p})\right)}$

where ${\displaystyle \Phi }$ is a smooth monotonic function. Writing ${\displaystyle g}$ for the inverse of ${\displaystyle \Phi }$, this is traditionally written as

${\displaystyle g(f({\vec {x}}))=\sum _{i}f_{i}(x_{i})}$.

When this function is approximating the expectation of some observed quantity, it could be written as

${\displaystyle g(\operatorname {E} (Y))=\beta _{0}+f_{1}(x_{1})+f_{2}(x_{2})+\cdots +f_{m}(x_{m}).\,\!}$

Which is the standard formulation of a generalized additive model. It was then shown[1] that the backfitting algorithm will always converge for these functions.

The GAM model class is quite broad, given that smooth function is a rather broad category. For example, a covariate ${\displaystyle x_{j}}$ may be multivariate and the corresponding ${\displaystyle f_{j}}$ a smooth function of several variables, or ${\displaystyle f_{j}}$ might be the function mapping the level of a factor to the value of a random effect. Another example is a varying coefficient (geographic regression) term such as ${\displaystyle z_{j}f_{j}(x_{j})}$ where ${\displaystyle z_{j}}$ and ${\displaystyle x_{j}}$ are both covariates. Or if ${\displaystyle x_{j}(t)}$ is itself an observation of a function, we might include a term such as ${\displaystyle \int f_{j}(t)x_{j}(t)dt}$ (sometimes known as a signal regression term). ${\displaystyle f_{j}}$ could also be a simple parametric function as might be used in any generalized linear model. The model class has been generalized in several directions, notably beyond exponential family response distributions, beyond modelling of only the mean and beyond univariate data.[2][3][4]

The original GAM fitting method estimated the smooth components of the model using non-parametric smoothers (for example smoothing splines or local linear regression smoothers) via the backfitting algorithm.[1] Backfitting works by iterative smoothing of partial residuals and provides a very general modular estimation method capable of using a wide variety of smoothing methods to estimate the ${\displaystyle f_{j}(x_{j})}$ terms. A disadvantage of backfitting is that it is difficult to integrate with the estimation of the degree of smoothness of the model terms, so that in practice the user must set these, or select between a modest set of pre-defined smoothing levels.

If the ${\displaystyle f_{j}(x_{j})}$ are represented using smoothing splines[5] then the degree of smoothness can be estimated as part of model fitting using generalized cross validation, or by restricted maximum likelihood (REML, sometimes known as 'GML') which exploits the duality between spline smoothers and Gaussian random effects.[6] This full spline approach carries an ${\displaystyle O(n^{3})}$ computational cost, where ${\displaystyle n}$ is the number of observations for the response variable, rendering it somewhat impractical for moderately large datasets. More recent methods have addressed this computational cost either by up front reduction of the size of the basis used for smoothing (rank reduction)[7][8][9][10][11] or by finding sparse representations of the smooths using Markov random fields, which are amenable to the use of sparse matrix methods for computation.[12] These more computationally efficient methods use GCV (or AIC or similar) or REML or take a fully Bayesian approach for inference about the degree of smoothness of the model components. Estimating the degree of smoothness via REML can be viewed as an empirical Bayes method.

An alternative approach with particular advantages in high dimensional settings is to use boosting, although this typically requires bootstrapping for uncertainty quantification.[13][14] GAMs fit using bagging and boosting have been found to generally outperform GAMs fit using spline methods.[15]

Many modern implementations of GAMs and their extensions are built around the reduced rank smoothing approach, because it allows well founded estimation of the smoothness of the component smooths at comparatively modest computational cost, and also facilitates implementation of a number of model extensions in a way that is more difficult with other methods. At its simplest the idea is to replace the unknown smooth functions in the model with basis expansions

${\displaystyle f_{j}(x_{j})=\sum _{k=1}^{K_{j}}\beta _{jk}b_{jk}(x_{j})}$

where the ${\displaystyle b_{jk}(x_{j})}$ are known basis functions, usually chosen for good approximation theoretic properties (for example B splines or reduced rank thin plate splines), and the ${\displaystyle \beta _{jk}}$ are coefficients to be estimated as part of model fitting. The basis dimension ${\displaystyle K_{j}}$ is chosen to be sufficiently large that we expect it to overfit the data to hand (thereby avoiding bias from model over-simplification), but small enough to retain computational efficiency. If ${\displaystyle p=\sum _{j}K_{j}}$ then the computational cost of model estimation this way will be ${\displaystyle O(np^{2})}$.

Notice that the ${\displaystyle f_{j}}$ are only identifiable to within an intercept term (we could add any constant to ${\displaystyle f_{1}}$ while subtracting it from ${\displaystyle f_{2}}$ without changing the model predictions at all), so identifiability constraints have to be imposed on the smooth terms to remove this ambiguity. Sharpest inference about the ${\displaystyle f_{j}}$ is generally obtained by using the sum-to-zero constraints

${\displaystyle \sum _{i}f_{j}(x_{ji})=0}$

i.e. by insisting that the sum of each the ${\displaystyle f_{j}}$ evaluated at its observed covariate values should be zero. Such linear constraints can most easily be imposed by reparametrization at the basis setup stage,[10] so below it is assumed that this has been done.

Having replaced all the ${\displaystyle f_{j}}$ in the model with such basis expansions we have turned the GAM into a Generalized linear model (GLM), with a model matrix that simply contains the basis functions evaluated at the observed ${\displaystyle x_{j}}$ values. However, because the basis dimensions, ${\displaystyle K_{j}}$, have been chosen to be a somewhat larger than is believed to be necessary for the data, the model is over-parameterized and will overfit the data if estimated as a regular GLM. The solution to this problem is to penalize departure from smoothness in the model fitting process, controlling the weight given to the smoothing penalties using smoothing parameters. For example, consider the situation in which all the smooths are univariate functions. Writing all the parameters in one vector, ${\displaystyle \beta }$, suppose that ${\displaystyle D(\beta )}$ is the deviance (twice the difference between saturated log likelihood and the model log likelihood) for the model. Minimizing the deviance by the usual iteratively re-weighted least squares would result in overfit, so we seek ${\displaystyle \beta }$ to minimize

${\displaystyle D(\beta )+\sum _{j}\lambda _{j}\int f_{j}^{\prime \prime }(x)^{2}dx}$

where the integrated square second derivative penalties serve to penalize wiggliness (lack of smoothness) of the ${\displaystyle f_{j}}$ during fitting, and the smoothing parameters ${\displaystyle \lambda _{j}}$ control the tradeoff between model goodness of fit and model smoothness. In the example ${\displaystyle \lambda _{j}\to \infty }$ would ensure that the estimate of ${\displaystyle f_{j}(x_{j})}$ would be a straight line in ${\displaystyle x_{j}}$.

Given the basis expansion for each ${\displaystyle f_{j}}$ the wiggliness penalties can be expressed as quadratic forms in the model coefficients.[10] That is we can write

${\displaystyle \int f_{j}^{\prime \prime }(x)^{2}dx=\beta _{j}^{T}{\bar {S}}_{j}\beta _{j}=\beta ^{T}S_{j}\beta }$,

where ${\displaystyle {\bar {S}}_{j}}$ is a matrix of known coefficients computable from the penalty and basis, ${\displaystyle \beta _{j}}$ is the vector of coefficients for ${\displaystyle f_{j}}$, and ${\displaystyle S_{j}}$ is just ${\displaystyle {\bar {S}}_{j}}$ padded with zeros so that the second equality holds and we can write the penalty in terms of the full coefficient vector ${\displaystyle \beta }$. Many other smoothing penalties can be written in the same way, and given the smoothing parameters the model fitting problem now becomes

${\displaystyle {\hat {\beta }}={\text{argmin}}_{\beta }\{D(\beta )+\sum _{j}\lambda _{j}\beta ^{T}S_{j}\beta \}}$,

which can be found using a penalized version of the usual iteratively reweighted least squares (IRLS) algorithm for GLMs: the algorithm is unchanged except that the sum of quadratic penalties is added to the working least squared objective at each iteration of the algorithm.

Penalization has several effects on inference, relative to a regular GLM. For one thing the estimates are subject to some smoothing bias, which is the price that must be paid for limiting estimator variance by penalization. However, if smoothing parameters are selected appropriately the (squared) smoothing bias introduced by penalization should be less than the reduction in variance that it produces, so that the net effect is a reduction in mean square estimation error, relative to not penalizing. A related effect of penalization is that the notion of degrees of freedom of a model has to be modified to account for the penalties' action in reducing the coefficients' freedom to vary. For example, if ${\displaystyle W}$ is the diagonal matrix of IRLS weights at convergence, and ${\displaystyle X}$ is the GAM model matrix, then the model effective degrees of freedom is given by ${\displaystyle {\text{trace}}(F)}$ where

${\displaystyle F=(X^{T}WX+\sum _{j}\lambda _{j}S_{j})^{-1}X^{T}WX}$,

is the effective degrees of freedom matrix.[10] In fact summing just the diagonal elements of ${\displaystyle F}$ corresponding to the coefficients of ${\displaystyle f_{j}}$ gives the effective degrees of freedom for the estimate of ${\displaystyle f_{j}}$.

Backfit GAMs were originally provided by the gam function in S,[22] now ported to the R language as the gam package. The SAS proc GAM also provides backfit GAMs. The recommended package in R for GAMs is mgcv, which stands for mixed GAM computational vehicle,[10] which is based on the reduced rank approach with automatic smoothing parameter selection. The SAS proc GAMPL is an alternative implementation. In Python, there is the InterpretML package, which implements a bagging and boosting approach.[23] There are many alternative packages. Examples include the R packages mboost,[13] which implements a boosting approach; gss, which provides the full spline smoothing methods;[24] VGAM which provides vector GAMs;[3] and gamlss, which provides Generalized additive model for location, scale and shape. `BayesX' and its R interface provides GAMs and extensions via MCMC and penalized likelihood methods.[25] The `INLA' software implements a fully Bayesian approach based on Markov random field representations exploiting sparse matrix methods.[12]

As an example of how models can be estimated in practice with software, consider R package mgcv. Suppose that our R workspace contains vectors y, x and z and we want to estimate the model

${\displaystyle y_{i}=\beta _{0}+f_{1}(x_{i})+f_{2}(z_{i})+\epsilon _{i}{\text{ where }}\epsilon _{i}\sim N(0,\sigma ^{2}).}$

Within R we could issue the commands

library(mgcv)  # load the package
b = gam(y ~ s(x) + s(z))

In common with most R modelling functions gam expects a model formula to be supplied, specifying the model structure to fit. The response variable is given to the left of the ~ while the specification of the linear predictor is given to the right. gam sets up bases and penalties for the smooth terms, estimates the model including its smoothing parameters and, in standard R fashion, returns a fitted model object, which can then be interrogated using various helper functions, such as summary, plot, predict, and AIC.

This simple example has used several default settings which it is important to be aware of. For example a Gaussian distribution and identity link has been assumed, and the smoothing parameter selection criterion was GCV. Also the smooth terms were represented using `penalized thin plate regression splines', and the basis dimension for each was set to 10 (implying a maximum of 9 degrees of freedom after identifiability constraints have been imposed). A second example illustrates how we can control these things. Suppose that we want to estimate the model

${\displaystyle y_{i}\sim {\text{Poi}}(\mu _{i}){\text{ where }}\log \mu _{i}=\beta _{0}+\beta _{1}x_{i}+f_{1}(t_{i})+f_{2}(v_{i},w_{i}).}$

using REML smoothing parameter selection, and we expect ${\displaystyle f_{1}}$ to be a relatively complicated function which we would like to model with a penalized cubic regression spline. For ${\displaystyle f_{2}}$ we also have to decide whether ${\displaystyle v}$ and ${\displaystyle w}$ are naturally on the same scale so that an isotropic smoother such as thin plate spline is appropriate (specified via `s(v,w)'), or whether they are really on different scales so that we need separate smoothing penalties and smoothing parameters for ${\displaystyle v}$ and ${\displaystyle w}$ as provided by a tensor product smoother. Suppose we opted for the latter in this case, then the following R code would estimate the model

b1 = gam(y ~ x + s(t,bs="cr",k=100) + te(v,w),family=poisson,method="REML")

which uses a basis size of 100 for the smooth of ${\displaystyle t}$. The specification of distribution and link function uses the `family' objects that are standard when fitting GLMs in R or S. Note that Gaussian random effects can also be added to the linear predictor.

These examples are only intended to give a very basic flavour of the way that GAM software is used, for more detail refer to the software documentation for the various packages and the references below.[10][24][3][22][13][25]

As with any statistical model it is important to check the model assumptions of a GAM. Residual plots should be examined in the same way as for any GLM. That is deviance residuals (or other standardized residuals) should be examined for patterns that might suggest a substantial violation of the independence or mean-variance assumptions of the model. This will usually involve plotting the standardized residuals against fitted values and covariates to look for mean-variance problems or missing pattern, and may also involve examining Correlograms (ACFs) and/or Variograms of the residuals to check for violation of independence. If the model mean-variance relationship is correct then scaled residuals should have roughly constant variance. Note that since GLMs and GAMs can be estimated using Quasi-likelihood, it follows that details of the distribution of the residuals beyond the mean-variance relationship are of relatively minor importance.

One issue that is more common with GAMs than with other GLMs is a danger of falsely concluding that data are zero inflated. The difficulty arises when data contain many zeroes that can be modelled by a Poisson or binomial with a very low expected value: the flexibility of the GAM structure will often allow representation of a very low mean over some region of covariate space, but the distribution of standardized residuals will fail to look anything like the approximate normality that introductory GLM classes teach us to expect, even if the model is perfectly correct.[26]

The one extra check that GAMs introduce is the need to check that the degrees of freedom chosen are appropriate. This is particularly acute when using methods that do not automatically estimate the smoothness of model components. When using methods with automatic smoothing parameter selection then it is still necessary to check that the choice of basis dimension was not restrictively small, although if the effective degrees of freedom of a term estimate is comfortably below its basis dimension then this is unlikely. In any case, checking ${\displaystyle f_{j}(x_{j})}$ is based on examining pattern in the residuals with respect to ${\displaystyle x_{j}}$. This can be done using partial residuals overlaid on the plot of ${\displaystyle {\hat {f}}_{j}(x_{j})}$, or using permutation of the residuals to construct tests for residual pattern (as in the `gam.check' function in R package `mgcv').

When smoothing parameters are estimated as part of model fitting then much of what would traditionally count as model selection has been absorbed into the fitting process: the smoothing parameters estimation has already selected between a rich family of models of different functional complexity. However smoothing parameter estimation does not typically remove a smooth term from the model altogether, because most penalties leave some functions un-penalized (e.g. straight lines are unpenalized by the spline derivative penalty given above). So the question of whether a term should be in the model at all remains. One simple approach to this issue is to add an extra penalty to each smooth term in the GAM, which penalizes the components of the smooth that would otherwise be unpenalized (and only those). Each extra penalty has its own smoothing parameter and estimation then proceeds as before, but now with the possibility that terms will be completely penalized to zero.[27] In high dimensional settings then it may make more sense to attempt this task using the Lasso (statistics) or Elastic net regularization. Boosting also performs term selection automatically as part of fitting.[13]

An alternative is to use traditional Stepwise regression methods for model selection. This is also the default method when smoothing parameters are not estimated as part of fitting, in which case each smooth term is usually allowed to take one of a small set of pre-defined smoothness levels within the model, and these are selected between in a stepwise fashion. Stepwise methods operate by iteratively comparing models with or without particular model terms (or possibly with different levels of term complexity), and require measures of model fit or term significance in order to decide which model to select at each stage. For example, we might use p-values for testing each term for equality to zero to decide on candidate terms for removal from a model, and we might compare Akaike information criterion (AIC) values for alternative models.

P-value computation for smooths is not straightforward, because of the effects of penalization, but approximations are available.[1][10] AIC can be computed in two ways for GAMs. The marginal AIC is based on the Mariginal Likelihood (see above) with the model coefficients integrated out. In this case the AIC penalty is based on the number of smoothing parameters (and any variance parameters) in the model. However, because of the well known fact that REML is not comparable between models with different fixed effects structures, we can not usually use such an AIC to compare models with different smooth terms (since their un-penalized components act like fixed effects). Basing AIC on the marginal likelihood in which only the penalized effects are integrated out is possible (the number of un-penalized coefficients now gets added to the parameter count for the AIC penalty), but this version of the marginal likelihood suffers from the tendency to oversmooth that provided the original motivation for developing REML. Given these problems GAMs are often compared using the conditional AIC, in which the model likelihood (not marginal likelihood) is used in the AIC, and the parameter count is taken as the effective degrees of freedom of the model.[1][21]

Naive versions of the conditional AIC have been shown to be much too likely to select larger models in some circumstances, a difficulty attributable to neglect of smoothing parameter uncertainty when computing the effective degrees of freedom,[28] however correcting the effective degrees of freedom for this problem restores reasonable performance.[2]

Overfitting can be a problem with GAMs,[21] especially if there is un-modelled residual auto-correlation or un-modelled overdispersion. Cross-validation can be used to detect and/or reduce overfitting problems with GAMs (or other statistical methods),[29] and software often allows the level of penalization to be increased to force smoother fits. Estimating very large numbers of smoothing parameters is also likely to be statistically challenging, and there are known tendencies for prediction error criteria (GCV, AIC etc.) to occasionally undersmooth substantially, particularly at moderate sample sizes, with REML being somewhat less problematic in this regard.[30]

Where appropriate, simpler models such as GLMs may be preferable to GAMs unless GAMs improve predictive ability substantially (in validation sets) for the application in question.

• Additive model
• Backfitting algorithm
• Generalized additive model for location, scale, and shape (GAMLSS)
• Residual effective degrees of freedom
• Semiparametric regression

• gam, an R package for GAMs by backfitting.
• gam, Python module in statsmodels.gam module.
• InterpretML, a Python package for fitting GAMs via bagging and boosting.
• mgcv, an R package for GAMs using penalized regression splines.
• mboost, an R package for boosting including additive models.
• gss, an R package for smoothing spline ANOVA.
• INLA software for Bayesian Inference with GAMs and more.
• BayesX software for MCMC and penalized likelihood approaches to GAMs.
• Doing magic and analyzing seasonal time series with GAM in R
• GAM: The Predictive Modeling Silver Bullet