Statistics models classIn statistics, a generalized additive model (GAM) is a generalized linear model in which the linear response variable 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. They can be interpreted as the discriminative generalization of the naive Bayes generative model.
The model relates a univariate response variable, "Y", to some predictor variables, "x""i". An exponential family distribution is specified for Y (for example normal, binomial or Poisson distributions) along with a link function "g" (for example the identity or log functions) relating the expected value of "Y" to the predictor variables via a structure such as
 formula_1
The functions "f""i" may be functions with a specified parametric form (for example a polynomial, or an un-penalized regression spline of a variable) 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 "f"1("x"1), and then use a factor model for "f"2("x"2). 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 (via the Kolmogorov–Arnold representation theorem) that any multivariate continuous function could be represented as sums and compositions of univariate functions,
formula_2.
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 drops the outer sum, and demands instead that the function belong to a simpler class,
formula_3.
where formula_4 is a smooth monotonic function. Writing formula_5 for the inverse of formula_4, this is traditionally written as
 formula_7.
When this function is approximating the expectation of some observed quantity, it could be written as
 formula_1
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 formula_9 may be multivariate and the corresponding formula_10 a smooth function of several variables, or formula_10 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 formula_12 where formula_13 and formula_9 are both covariates. Or if formula_15 is itself an observation of a function, we might include a term such as formula_16 (sometimes known as a signal regression term). formula_17 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 (for example smoothing splines or local linear regression 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 formula_18 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 formula_18 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 (REML, sometimes known as 'GML') which exploits the duality between spline smoothers and Gaussian random effects. This full spline approach carries an formula_20 computational cost, where formula_21 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) 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 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. 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 
formula_22
where the formula_23 are known basis functions, usually chosen for good approximation theoretic properties (for example B splines or reduced rank thin plate splines), and the formula_24 are coefficients to be estimated as part of model fitting. The basis dimension formula_25 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 formula_26 then the computational cost of model estimation this way will be formula_27.
Notice that the formula_17 are only identifiable to within an intercept term (we could add any constant to formula_29 while subtracting it from formula_30 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 formula_17 is generally obtained by using the sum-to-zero constraints
formula_32 
i.e. by insisting that the sum of each the formula_17 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 formula_34 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 formula_35 values. However, because the basis dimensions, formula_25, 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, formula_37, suppose that formula_38 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 formula_39 to minimize
formula_40 
where the integrated square second derivative penalties serve to penalize wiggliness (lack of smoothness) of the formula_17 during fitting, and the smoothing parameters formula_42 control the tradeoff between model goodness of fit and model smoothness. In the example formula_43 would ensure that the estimate of formula_18 would be a straight line in formula_35.
Given the basis expansion for each formula_17 the wiggliness penalties can be expressed as quadratic forms in the model coefficients. That is we can write
formula_47, 
where formula_48 is a matrix of known coefficients computable from the penalty and basis, formula_49 is the vector of coefficients for formula_17, and formula_51 is just formula_48 padded with zeros so that the second equality holds and we can write the penalty in terms of the full coefficient vector formula_39. Many other smoothing penalties can be written in the same way, and given the smoothing parameters the model fitting problem now becomes
formula_54,
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 formula_55 is the diagonal matrix of IRLS weights at convergence, and formula_56 is the GAM model matrix, then the model effective degrees of freedom is given by formula_57 where
formula_58,
is the effective degrees of freedom matrix. In fact summing just the diagonal elements of formula_59 corresponding to the coefficients of formula_17 gives the effective degrees of freedom for the estimate of formula_17.
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 
formula_62
(where formula_63 is the GLM scale parameter introduced only for later convenience), but we can immediately recognize this as a multivariate normal prior with mean formula_64 and precision matrix formula_65. Since the penalty allows some functions through unpenalized (straight lines, given the example penalties), formula_66 is rank deficient, and the prior is actually improper, with a covariance matrix given by the Moore–Penrose pseudoinverse of formula_66 (the impropriety corresponds to ascribing infinite variance to the unpenalized components of a smooth).
Now if this prior is combined with the GLM likelihood, we find that the posterior mode for formula_39 is exactly the formula_69 found above by penalized IRLS. Furthermore, we have the large sample result that 
formula_70
which can be used to produce confidence/credible intervals for the smooth components, formula_17.
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, formula_72, but these also need to be estimated. One approach is to take a fully Bayesian approach, defining priors on the (log) 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 (GCV) or the
Akaike information criterion (AIC). Finally we may choose to maximize the Marginal Likelihood (REML) obtained by integrating the model coefficients, formula_39 out of the joint density of formula_74, 
formula_75.
Since formula_76 is just the likelihood of formula_39, we can view this as choosing formula_78 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 (log) smoothing parameter vector requiring a penalized IRLS iteration to evaluate