# Multiple polynomial regression versus GAM

How does using GAM differ from using multiple polynomial regression? They seem to produce the same result. Below I run a polynomial regression using lm() with

$wage=\beta_0 + \beta_{age}\cdot age+\beta_{age^2}\cdot age^2 + \beta_{year}\cdot year+\beta_{year^2}\cdot year^2$

Then I do the same with gam() and compare the 2 models using anova(). The RSS for both is identical. Same using splines of $year$ and $age$.

If I am trying to explore the degree of nonlinearity between response and predictors, how would I do it differently using gam()? Is gam() meant to make it easier to compare a quadratic model to a cubic one? Thank you.

library(ISLR) # for Wage dataset
library(gam)

lm.wage.poly <- lm(wage ~ poly(age, 2, raw = TRUE) + poly(year, 2, raw = TRUE), data = Wage)
lm.wage.gam <- gam(wage ~ poly(age, 2, raw = TRUE) + poly(year, 2, raw = TRUE), data = Wage)
anova(lm.wage.poly, lm.wage.gam)
## Model 1: wage ~ poly(age, 2, raw = TRUE) + poly(year, 2, raw = TRUE)
## Model 2: wage ~ poly(age, 2, raw = TRUE) + poly(year, 2, raw = TRUE)

##   Res.Df     RSS Df Sum of Sq F Pr(>F)

## 1   2995 4771710
## 2   2995 4771710  0         0

library(splines)
# Using cubic splines for age and year
lm.wage.poly.spline <- lm(wage ~ bs(age) + bs(year), data = Wage)
lm.wage.gam.spline <- gam(wage ~ bs(age) + bs(year), data = Wage)
anova(lm.wage.poly.spline, lm.wage.gam.spline)

## Model 1: wage ~ bs(age) + bs(year)

## Model 2: wage ~ bs(age) + bs(year)

##   Res.Df     RSS Df Sum of Sq F Pr(>F)
## 1   2993 4754108
## 2   2993 4754108  0         0

• Not all generalized additive models have the form of a polynomial. Commented May 18, 2017 at 15:51

The two models you fitted are exactly equivalent (up to implementational details). The advantage of the GAM framework is that you are not limited to global basis expansions of your covariates. Instead you can use a range of penalized spline bases that may better adapt to the data rather than imposing a particular functional form as you did with the polynomial model.

I also wouldn't use the gam package. The penalized spline approach of Simon Wood, as implemented in the mgcv package that ships with R is much more useful that Hastie and Tibshirani's original.

Here I modify your example to use (thin plate) splines in the GAM via the s() function (for smooth) rather than the poly() function used in the lm(). I restrict the basis dimension for both terms but set them larger than the polynomial for technical reasons (if the true curves are polynomials of degree 2, that is on the edge of the set of functions the spline can represent and smoothness selection may not identify them — I come back to the issue of basis dimension size later). Here the splines use maximum 4 degrees of freedom, from k = 5 but we loose 1 df each due to identifiability constraints to accommodate a model intercept.

library("ISLR") # for Wage dataset
data(Wage)
library("mgcv")

lm.wage.poly <- lm(wage ~ poly(age, 2, raw = TRUE) + poly(year, 2, raw = TRUE), data = Wage)
lm.wage.gam <- gam(wage ~ s(age, k = 5) + s(year, k = 5), data = Wage)

anova(lm.wage.poly, lm.wage.gam)


That gives:

> anova(lm.wage.poly, lm.wage.gam)
Analysis of Variance Table

Model 1: wage ~ poly(age, 2, raw = TRUE) + poly(year, 2, raw = TRUE)
Model 2: wage ~ s(age, k = 5) + s(year, k = 5)
Res.Df     RSS      Df Sum of Sq      F    Pr(>F)
1 2995.0 4771710
2 2994.2 4745169 0.84022     26541 19.932 3.083e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


which indicates that the GAM uses slightly more degrees of freedom than the polynomial model, but has improved the fit to a degree.

If we look at the model summary

> summary(lm.wage.gam)

Family: gaussian

Formula:
wage ~ s(age, k = 5) + s(year, k = 5)

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 111.7036     0.7268   153.7   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df     F  p-value
s(age)  3.84  3.984 70.97  < 2e-16 ***
s(year) 1.00  1.000 13.58 0.000232 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.0899   Deviance explained = 9.13%
GCV = 1587.9  Scale est. = 1584.8    n = 3000


We note that gam() selected a linear function of year (it uses 1 effective degree of freedom == EDF). An ~ 4 EDF term was fitted to age. This is in contrast to the polynomial model, which used 2 DFs per covariate regardless.

We should check that the basis dimensions specified (k = 5) were sufficiently large:

> gam.check(lm.wage.gam)

Method: GCV   Optimizer: magic
Smoothing parameter selection converged after 9 iterations.
The RMS GCV score gradient at convergence was 0.0001096254 .
The Hessian was positive definite.
Model rank =  9 / 9

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

k'  edf k-index p-value
s(age)  4.00 3.84    0.98    0.09 .
s(year) 4.00 1.00    1.01    0.72
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


The test for s(age) is marginal and we would be advised to increase k and refit. Here I increase both terms basis dimension, but year only has a few unique values so we can't go beyond k = 7 for that term.

m2 <- gam(wage ~ s(age, k = 10) + s(year, k = 7), data = Wage)


Now smoothness selection suggests age should be more strongly non-linear that before:

> summary(m2)

Family: gaussian

Formula:
wage ~ s(age, k = 10) + s(year, k = 7)

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 111.7036     0.7266   153.7   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df     F  p-value
s(age)  5.462  6.568 43.51  < 2e-16 ***
s(year) 1.000  1.000 14.18 0.000169 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.0905   Deviance explained = 9.24%
GCV = 1587.7  Scale est. = 1583.7    n = 3000


The fitted smooths are now:

The conceptual difference between the polynomial basis expansion and the spline basis expansion used in the GAM is that the former is global; data at age 20 contribute to the fit of the quadratic function of age at age 80. Because of this global nature, the polynomial can be too inflexible for many nonlinear relationships. As the spline basis is local, the fit at a particular age for example is largely dominated by the values of the response around the same age.

If we use lo and s basis expansion in gam, we will see the difference. Here is the demo:

Using no basis expansion and polynomial, gam and lm are the same.

> library(gam)
> anova(lm(mpg~wt,data=mtcars),gam(mpg~wt,data=mtcars))
Analysis of Variance Table

Model 1: mpg ~ wt
Model 2: mpg ~ wt
Res.Df    RSS Df  Sum of Sq F Pr(>F)
1     30 278.32
2     30 278.32  0 5.6843e-14
> anova(lm(mpg~poly(wt,2,raw = T),data=mtcars),gam(mpg~poly(wt,2,raw = T),data=mtcars))
Analysis of Variance Table

Model 1: mpg ~ poly(wt, 2, raw = T)
Model 2: mpg ~ poly(wt, 2, raw = T)
Res.Df    RSS Df Sum of Sq F Pr(>F)
1     29 203.75
2     29 203.75  0         0


Using lo and s expansion, gam and lm are different.

> anova(lm(mpg~lo(wt),data=mtcars), gam(mpg~lo(wt),data=mtcars))
Analysis of Variance Table

Model 1: mpg ~ lo(wt)
Model 2: mpg ~ lo(wt)
Res.Df    RSS  Df Sum of Sq      F   Pr(>F)
1   30.0 278.32
2   26.6 177.47 3.4    100.86 4.4462 0.009422 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(lm(mpg~s(wt),data=mtcars), gam(mpg~s(wt),data=mtcars))
Analysis of Variance Table

Model 1: mpg ~ s(wt)
Model 2: mpg ~ s(wt)
Res.Df    RSS Df Sum of Sq      F  Pr(>F)
1     30 278.32
2     27 190.00  3    88.321 4.1837 0.01483 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1