Generalized additive model


In statistics, a generalized additive model is a generalized linear model in which the linear predictor depends linearly on unknown smooth functions of some predictor variables, and interest focuses on inference about these smooth functions.
GAMs were originally developed by Trevor Hastie and Robert Tibshirani to blend properties of generalized linear models with additive models.
The model relates a univariate response variable, Y, to some predictor variables, xi. An exponential family distribution is specified for Y along with a link function g relating the expected value of Y to the predictor variables via a structure such as
The functions fi may be functions with a specified parametric form or may be specified non-parametrically, or semi-parametrically, simply as 'smooth functions', to be estimated by non-parametric means. So a typical GAM might use a scatterplot smoothing function, such as a locally weighted mean, for f1, and then use a factor model for f2. This flexibility to allow non-parametric fits with relaxed assumptions on the actual relationship between response and predictor, provides the potential for better fits to data than purely parametric models, but arguably with some loss of interpretability.

Theoretical background

It had been known since the 1950s that any multivariate function could be represented as sums and compositions of univariate functions.
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 functions, and thus are not suitable for modeling approaches. Therefore, the generalized additive model drops the outer sum, and demands instead that the function belong to a simpler class.
where is a smooth monotonic function. Writing for the inverse of, this is traditionally written as
When this function is approximating the expectation of some observed quantity, it could be written as
Which is the standard formulation of a generalized additive model. It was then shown that the backfitting algorithm will always converge for these functions.

Generality

The GAM model class is quite broad, given that smooth function is a rather broad category. For example, a covariate may be multivariate and the corresponding a smooth function of several variables, or might be the function mapping the level of a factor to the value of a random effect. Another example is a varying coefficient term such as where and are both covariates. Or if is itself an observation of a function, we might include a term such as . 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.

GAM fitting methods

The original GAM fitting method estimated the smooth components of the model using non-parametric smoothers via the backfitting algorithm. 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 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 are represented using smoothing splines then the degree of smoothness can be estimated as part of model fitting using generalized cross validation, or by restricted maximum likelihood which exploits the duality between spline smoothers and Gaussian random effects. This full spline approach carries an computational cost, where 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 or by finding sparse representations of the smooths using Markov random fields, which are amenable to the use of sparse matrix methods for computation. These more computationally efficient methods use GCV 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. GAMs fit using bagging and boosting have been found to generally outperform GAMs fit using spline methods.

The rank reduced framework

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
where the are known basis functions, usually chosen for good approximation theoretic properties, and the are coefficients to be estimated as part of model fitting. The basis dimension is chosen to be sufficiently large that we expect it to overfit the data to hand, but small enough to retain computational efficiency. If then the computational cost of model estimation this way will be.
Notice that the are only identifiable to within an intercept term, so identifiability constraints have to be imposed on the smooth terms to remove this ambiguity. Sharpest inference about the is generally obtained by using the sum-to-zero constraints
i.e. by insisting that the sum of each the evaluated at its observed covariate values should be zero. Such linear constraints can most easily be imposed by reparametrization at the basis setup stage, so below it is assumed that this has been done.
Having replaced all the in the model with such basis expansions we have turned the GAM into a Generalized linear model, with a model matrix that simply contains the basis functions evaluated at the observed values. However, because the basis dimensions,, 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,, suppose that is the deviance for the model. Minimizing the deviance by the usual iteratively re-weighted least squares would result in overfit, so we seek to minimize
where the integrated square second derivative penalties serve to penalize wiggliness of the during fitting, and the smoothing parameters control the tradeoff between model goodness of fit and model smoothness. In the example would ensure that the estimate of would be a straight line in.
Given the basis expansion for each the wiggliness penalties can be expressed as quadratic forms in the model coefficients. That is we can write
where is a matrix of known coefficients computable from the penalty and basis, is the vector of coefficients for , and is just padded with zeros so that the second equality holds and we can write the penalty in terms of the full coefficient vector. Many other smoothing penalties can be written in the same way, and given the smoothing parameters the model fitting problem now becomes
which can be found using a penalized version of the usual iteratively reweighted least squares 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 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 is the diagonal matrix of IRLS weights at convergence, and is the GAM model matrix, then the model effective degrees of freedom is given by where
is the effective degrees of freedom matrix. In fact summing just the diagonal elements of corresponding to the coefficients of gives the effective degrees of freedom for the estimate of .

Bayesian smoothing priors

Smoothing bias complicates interval estimation for these models, and the simplest approach turns out to involve a Bayesian approach. Understanding this Bayesian view of smoothing also helps to understand the REML and full Bayes approaches to smoothing parameter estimation. At some level smoothing penalties are imposed because we believe smooth functions to be more probable than wiggly ones, and if that is true then we might as well formalize this notion by placing a prior on model wiggliness. A very simple prior might be
, but we can immediately recognize this as a multivariate normal prior with mean and precision matrix. Since the penalty allows some functions through unpenalized, is rank deficient, and the prior is actually improper, with a covariance matrix given by the Moore-Penrose pseudoinverse of .
Now if this prior is combined with the GLM likelihood, we find that the posterior mode for is exactly the found above by penalized IRLS. Furthermore, we have the large sample result that
which can be used to produce confidence/credible intervals for the smooth components,.
The Gaussian smoothness priors are also the basis for fully Bayesian inference with GAMs, as well as methods estimating GAMs as mixed models that are essentially empirical Bayes methods.

Smoothing parameter estimation

So far we have treated estimation and inference given the smoothing parameters,, but these also need to be estimated. One approach is to take a fully Bayesian approach, defining priors on the smoothing parameters, and using stochastic simulation or high order approximation methods to obtain information about the posterior of the model coefficients. An alternative is to select the smoothing parameters to optimize a prediction error criterion such as Generalized cross validation or the
Akaike information criterion. Finally we may choose to maximize the Marginal Likelihood obtained by integrating the model coefficients, out of the joint density of,
Since is just the likelihood of, we can view this as choosing to maximize the average likelihood of random draws from the prior. The preceding integral is usually analytically intractable but can be approximated to quite high accuracy using Laplace's method.
Smoothing parameter inference is the most computationally taxing part of model estimation/inference. For example, to optimize a GCV or marginal likelihood typically requires numerical optimization via a Newton or Quasi-Newton method, with each trial value for the smoothing parameter vector requiring a penalized IRLS iteration to evaluate the corresponding alongside the other ingredients of the GCV score or Laplace approximate marginal likelihood. Furthermore, to obtain the derivatives of the GCV or LAML, required for optimization, involves implicit differentiation to obtain the derivatives of w.r.t. the log smoothing parameters, and this requires some care is efficiency and numerical stability are to be maintained.

Software

Backfit GAMs were originally provided by the gam function in S, 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, 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. There are many alternative packages. Examples include the R packages mboost, which implements a boosting approach; gss, which provides the full spline smoothing methods; VGAM which provides vector GAMs; 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. The `INLA' software implements a fully Bayesian approach based on Markov random field representations exploiting sparse matrix methods.
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
Within R we could issue the commands
library # load the package
b = gam + s)
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. A second example illustrates how we can control these things. Suppose that we want to estimate the model
using REML smoothing parameter selection, and we expect to be a relatively complicated function which we would like to model with a penalized cubic regression spline. For we also have to decide whether and are naturally on the same scale so that an isotropic smoother such as thin plate spline is appropriate, or whether they are really on different scales so that we need separate smoothing penalties and smoothing parameters for and 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 + te
which uses a basis size of 100 for the smooth of. 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.

Model checking

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 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 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.
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 is based on examining pattern in the residuals with respect to. This can be done using partial residuals overlaid on the plot of, or using permutation of the residuals to construct tests for residual pattern.

Model selection

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 of the fact that most penalties leave some functions un-penalized. 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. 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. In high dimensional settings then it may make more sense to attempt this task using the Lasso or Elastic net regularization. Boosting also performs term selection automatically as part of fitting.
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, 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 values for alternative models.
P-value computation for smooths is not straightforward, because of the effects of penalization, but approximations are available. AIC can be computed in two ways for GAMs. The marginal AIC is based on the Mariginal Likelihood with the model coefficients integrated out. In this case the AIC penalty is based on the number of smoothing 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. Basing AIC on the marginal likelihood in which only the penalized effects are integrated out is possible, 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 is used in the AIC, and the parameter count is taken as the effective degrees of freedom of the model.
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, however correcting the effective degrees of freedom for this problem restores reasonable performance.

Caveats

can be a problem with GAMs, 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, 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 to occasionally undersmooth substantially, particularly at moderate sample sizes, with REML being somewhat less problematic in this regard.
Where appropriate, simpler models such as GLMs may be preferable to GAMs unless GAMs improve predictive ability substantially for the application in question.