Generalized Additive Models

Abstract

Regression models play an important role in many applied settings, by enabling predictive analysis, revealing classification rules, and providing data-analytic tools for understanding the interactions between different variables. Although attractively simple, traditional linear regression often fails in practical situations because real-life effects are usually nonlinear. This article describes statistical methods that may be used to identify and characterize general nonlinear regressions, without requiring the analyst to prespecify the form of the nonlinear relationship. These methods form the basis of the generalized additive models approach to data analysis. The choice of the link function g permits the analyst to incorporate many of the benefits of the methodology of generalized linear models. The logit link is widely viewed as highly appropriate for the analysis of binary data, and other standard link functions are similarly useful. Besides the identity link g(μ) = μ, which yields the additive model corresponding to ordinary multiple regression, people often use g(μ) = lnμ to give a log-additive model for Poisson count data (See Poisson Distribution). These link functions arise from exponential family sampling models, which also include the gamma and negative‐binomial distributions. These families generate the well-known class of generalized linear models 12, which may all be extended in the same way to generalized additive models by replacing linear terms with nonparametric functions. In this section we describe a modular algorithm for fitting additive models and their generalizations. The basic component is the scatterplot smoother for fitting nonlinear curves in a way that permits a great variety of functional forms. Suppose that we have a scatterplot of points (x i , y i ) like that shown in Fig. 1. Here y is a response or outcome variable, and x is a predictor. We wish to fit a smooth curve f(x) that summarizes the dependence of y on x (see Curve Fitting). More formally, we want to fit the model y = f(x) + ε, where f(x) is specified in a flexible way, yet required to be reasonably smooth. If we were to find the curve that simply minimized ∑[y i − f(x i )]2, the result would be an interpolating curve that would not be smooth at all. Notice that ∫f″(x)2 dx measures the “wiggliness” of the function f: linear functions f have ∫f″(x)2 dx = 0, while nonlinear functions produce values larger than zero. Here λ is a nonnegative smoothing parameter that must be chosen by the data analyst. It governs the trade-off between the goodness of fit* to the data (as measured by ∑[y i − f(x i )]2) and the wiggliness of the function. Larger values of λ force f to be smoother. In fact the interpolating curve corresponds to λ = 0 at one extreme, and the straight line fit is the limit as λ → ∞. For any value of λ, the solution to 1 is a cubic spline, i.e., a piecewise cubic polynomial with pieces joined smoothly at the values of x in the data set. Fast and stable numerical procedures are available for computation of the fitted curve. The right panel of Fig. 1 shows a cubic-spline fit to the data. What value of λ did we use in Fig. 1 ? In fact it is not convenient to express the desired smoothness of f in terms of λ, since the meaning of λ depends on the units of the explanatory variable x. Instead, it is possible to define an “effective number of parameters” or “degrees of freedom” for a cubic spline smoother, and then use a numerical search to determine the value of λ that yields this number. In Fig. 1 we chose the degrees of freedom to be five. Roughly speaking, this means that the complexity of the curve is about the same as that of a polynomial regression of degree four. However, the cubic spline smoother “spreads out” its parameters in a more even manner, and hence allows more diverse curves than a polynomial regression. Note that the degrees of freedom of a smoother need not be an integer. The preceding discussion sought to minimize the least squares fit to the data, subject to a penalty on the roughness of the fitted function. For the logistic regression model and other generalized additive models, the appropriate criterion is penalized maximum likelihood or penalized maximum log partial likelihood. To maximize these, the backfitting procedure is used in conjunction with a maximum likelihood or maximum partial likelihood algorithm. The usual Newton–Raphson method for maximizing log likelihoods in these models can be cast in an iteratively reweighted least squares form. This involves repeated use of weighted linear regression of a constructed response variable on the covariates: each regression yields a new value of the parameter estimates which give a new constructed variable, and the process is iterated. In the generalized additive model, the weighted linear regression is simply replaced by a weighted backfitting algorithm. We describe the algorithm in more detail for logistic regression below, and in more generality in Chap. 6 of 8. Iterate: Construct an adjusted dependent variable Construct weights w i = p i old ( 1 − p i old ) . Compute η new = A w z , the weighted additive model fit to z , where A w is the matrix whose rows represent estimates of β0 and the values of the unknown functions f j at the corresponding covariate values. If the warming procedure had been given in a randomized manner, we could simply focus on the post-February 1988 data and compare the survival rate of those who received the new procedure with those who did not. However, allocation was not random, so we can only try to assess the effectiveness of the warming procedure as it was applied. For this analysis, we use all of the data (1983–1988). To adjust for changes that might have occurred over the five-year period, we include the date of the operation as a covariate. But operation date is strongly confounded with the use of the warming operation, and thus a general nonparametric fit for date of operation might unduly remove some of the effect attributable to the warming procedure. To avoid this, we allow only a linear effect (increasing or decreasing) for operation date. Hence we assume that any time trend is linear. We fit a generalized additive logistic model to the binary response death, with smooth terms for age and weight, a linear term for operation date, a categorical variable for diagnosis, and a binary variable for the warming operation. To start the algorithm, all the smooth terms are fitted with four degrees of freedom. For interpretation, one wants to understand the influence each covariate has on the logit transform of the probability of death. The curves for age and weight are shown in Fig. 2. The tick marks on the x-axis indicate the points at which data were available. As one would expect, the highest risk is for the lighter babies, with a decreasing risk over 3 kg. Somewhat surprisingly, there seems to be a low-risk age around 200 days, with higher risk for younger and older children. In Table 1 each line summarizes the fit for the factor listed in the column on the left; here diag1–diag5 are the five indicator variables for the six diagnosis categories, and df is the number of degrees of freedom used for that variable. Note that the numerical algorithm is unable to achieve exactly 4 degrees of freedom for the age and weight effects, but rather uses 3.80 and 3.86 degrees of freedom, respectively. For ease of interpretation, the estimated curves for age and weight are decomposed into a linear component and the remaining nonlinear component (the linear component is essentially obtained as a weighted least-squares linear fit of the estimated curve on the predictor, while the nonlinear part is the residual). Then the next three columns, headed coefficient, standard error, and z-score, report the estimated coefficient, the standard error of the coefficient, and the normal score, respectively, for the linear component of the corresponding factor. The last column, headed Nonlinear p-Value, gives the p-value for a test of nonlinearity of the effect. Note however that the z-scores of the other factors (e.g., treatment) are fully adjusted for the shapes of all fitted curves, not just for their linear parts. We find that the warming procedure is strongly significant, with an estimated coefficient of 1.43 and a standard error of 0.45, indicating a survival benefit. There are strong differences in the diagnosis categories, while the estimated effect of the operation date is not large. Since a logistic regression is additive on the logit scale but not on the probability scale, a plot of the fitted probabilities is often informative. Figure 3 shows the fitted probabilities broken down by age and diagnosis, and is a concise summary of the findings of this study. The beneficial effect of the treatment at the lower weights is evident. As with all nonrandomized studies, the results here should be interpreted with caution. In particular, one must ensure that the children were not chosen for the warming operation according to their prognosis. To investigate this, we performed a second analysis in which a dummy variable, say period, corresponding to before versus after February 1988, was inserted in place of the dummy variable for the warming operation treatment. The purpose of this was to determine whether the overall treatment strategy improved after February 1988. If this turned out not to be the case, it would imply that warming was used only for patients with a good prognosis, who would have survived anyway. A linear adjustment for the operation date was included as before. The results were qualitatively very similar to the first analysis: age and weight were significant, with effects similar to those in Fig. 2; the diagnosis category was significant, while the operation data (a linear effect) was not. The period was highly significant, with a coefficient of −1.12 and a standard error of 0.33. Hence there seems to be a significant overall improvement in survival after February 1988 15. Software for fitting generalized additive models is available as part of the S/S-PLUS statistical language 1, 3, in a Fortran program called gamfit available from Statlib (in general/gamfit at ftp site lib.stat.cmu.edu), and also in the GAIM package for MS-DOS computers (information on this is available from the authors). Computer‐Intensive Statistical Methods; Nonparametric Regression.