17.1 Generalized Additive Models

Generalized additive models are a kind of linear regression, but instead of finding coefficients of predictor variables (e.g., intercepts, slopes), the model finds a “smooth” response function for each predictor. Typically this means that a piecewise cubic function (spline) is used to approximate the relationship between two variables. We can compute predicted values, confindence and prediction intervals, and show the smooth response function that arose from the model. We don’t provide a simple list of coefficients like we did with linear regression, because the spline curve is defined by many numbers which are not usually informative on their own.

Generalized additive models are most commonly used when there is no theoretical motivation for a functional relationship between the variables being studied. We’ll look at the Mauna Loa atmospheric CO2 concentration. These data increase year over year and have well-established interannual oscillations, but there is no clear function for either of these patterns.

We will start by using the last decade of data from Lesson 1.

Here is a plot of atmospheric CO2 concentration over time.

my_theme = theme_linedraw() + theme(text = element_text(size = 14))
p1 <- co2 %>% ggplot(aes(decimal_date, co2_avg)) + geom_point() + my_theme +
  labs(x = "Year (decimal)", y = "Atmospheric CO2 (ppm)")
p1

The next plot shows each observation, minus the annual mean, as a function of the time of year but not the actual year. Each year’s observations are overlapped.

p2 <- co2 %>% ggplot(aes(year_fraction, co2_anomaly)) + geom_point() + my_theme +
  labs(x = "Year fraction (decimal)", y = "Atmospheric CO2 anomaly (ppm)")
p2

Here are GAM fits to both of these pairs of data using the gam function from the mgcv package. There is a plot function in the mgcv package, but I’m using a ggplot function called draw in the gratia package.

g1 <- gam(co2_anomaly ~ s(year_fraction, bs="cs"), data = co2)
draw(g1, residuals=TRUE, rug=FALSE)

g2 <- gam(co2_avg ~ s(decimal_date, bs="cs"),  data = co2)
draw(g2, residuals=TRUE, rug=FALSE) 

Depending on your goal, you may find that the second plot is too smooth – the interannual oscillations are smoothed out completely. We can increase the number of “knots” (by setting k=25) in the spline to capture the oscillation, but this is not the way the smoothing is normally used:

g3 <- gam(co2_avg ~ s(decimal_date, bs="cs", k = 25),  data = co2)
draw(g3, residuals=TRUE, rug=FALSE) 

There is a summary function for GAM fits and the broom functions glance and tidy give access to these data, but the interpretation of this output is beyond the scope of the course:

summary(g1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## co2_anomaly ~ s(year_fraction, bs = "cs")
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.364e-15  2.673e-02       0        1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df     F p-value    
## s(year_fraction) 8.296      9 530.7  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.976   Deviance explained = 97.8%
## GCV = 0.092244  Scale est. = 0.085038  n = 119

glance(g1) %>% kable()

tidy(g1) %>% kable()

It is often useful to compute residuals and extract predicted values.

co2 %>% add_residuals(g1) %>% add_fitted(g1) %>% kable() %>%  scroll_box(height = "200px")

If we generate new data (dates), we can plot predictions too. Since the model is a piecewise cubic function, extrapolations are often dramatically unreliable. The downward facing cubic in the last “bump” is simply continued, with comically bad results. Extrapolation of models, and especially smooths, is somewhere between risky and foolish!

new_data <- tibble(decimal_date = seq(2017, 2022, by = 0.05))
new_data %>% add_fitted(g3) %>%
  ggplot(aes(decimal_date, .value)) + 
  geom_line(color = "blue", size = 2) +
  geom_point(aes(y= co2_avg), data = co2) + my_theme + xlim(2017, 2022)

If you want confidence intervals on the fitted values, use the confint function together with the name of the smooth you are extracting. Be aware that this function does not include the intercept (or grand mean) from the model, so the values are all centred on zero.

confint(g1, "s(year_fraction)", level = 0.95) %>% kable() %>%  scroll_box(height = "200px")

confint(g1, "s(year_fraction)", level = 0.95) %>% 
  ggplot(aes(year_fraction)) + 
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.25) + 
  geom_line(aes(y = est)) + my_theme

You can also use lm to fit splines; these are similar to GAMs but there are some important differences. The mgcv package has a lot of features not available with lm.

s1 <- lm(co2_avg ~ splines::bs(decimal_date, df = 5), data = co2)
tidy(s1) %>% kable(digits = 2)

glance(s1) %>% kable(digits = 2)

We can use augment to generate data to plot and combine it with the original data.

a1 <- augment(s1, data = co2, se_fit = TRUE, interval="prediction") 
a1 %>%  ggplot(aes(decimal_date)) + 
  geom_ribbon(aes(ymin = .lower, ymax = .upper), alpha = 0.2) + 
  geom_line(aes(y= .fitted)) +
  geom_point(aes(y = co2_avg)) + 
  my_theme

For a spline that closely traces the data, increase df to 26 or more. This increases the numbers of knots or separate cubics used to approximate the data.

17.2 Locally Estimated Scatterplot Smoothing (LOESS)

LOESS smooths are constructed by making a large number of quadratic (or possibly linear) regression lines as a window moves along the x-axis. The degree of the fits and the width of the window (and other details) can be adjusted. The predictions from these many local regressions are then computed and plotted as a “smooth” of the data. We usually just imagine the line is drawn through the data in a way that allows it to track fluctuations without specifying a model.

l1 <- loess(co2_anomaly ~ year_fraction, data = co2)

As with the GAMs, the summary function is not particularly easy to interpret and we will not explore the details, but you can see that this is a quadratic model with a “span” of 0.75, meaning that 75% of the data are used for each regression. (The weighting of each point in the regression varies as a function of x, resulting a continuous change in the predictions.)

summary(l1)

## Call:
## loess(formula = co2_anomaly ~ year_fraction, data = co2)
## 
## Number of Observations: 119 
## Equivalent Number of Parameters: 4.58 
## Residual Standard Error: 0.4647 
## Trace of smoother matrix: 5  (exact)
## 
## Control settings:
##   span     :  0.75 
##   degree   :  2 
##   family   :  gaussian
##   surface  :  interpolate      cell = 0.2
##   normalize:  TRUE
##  parametric:  FALSE
## drop.square:  FALSE

We can use predict, residuals and augment on these model objects as we did for lm fits in the previous lesson.

augment(l1, se_fit = TRUE) %>% kable() %>%  scroll_box(height = "200px")

augment(l1, se_fit = TRUE) %>%
  ggplot(aes(x = year_fraction, y = .fitted)) + 
  geom_line() + 
  geom_ribbon(aes(ymin = .fitted - 2*.se.fit, ymax = .fitted + 2*.se.fit), alpha = 0.20) + 
  my_theme

This line is only drawn using 12 points, so we might want to generate a smoother prediction by creating a new set of dates. We also add the original data to the plot.

new_data = tibble(year_fraction = seq(min(co2$year_fraction), 
                                      max(co2$year_fraction),
                                      length = 100))
augment(l1, newdata = new_data, se_fit = TRUE) %>%
  ggplot(aes(x = year_fraction, y = .fitted)) + 
  geom_line() + 
  geom_ribbon(aes(ymin = .fitted - 2*.se.fit,
                  ymax = .fitted + 2*.se.fit), 
              alpha = 0.20) + 
  geom_point(aes(y = co2_anomaly), data = co2) +
  my_theme

The LOESS function will not predict outside of the range of data provided, so I had to select the range of new_data above carefully to get the line to start and end near the first and last points.

Now we make two plots that contrast the results with different sized windows (span) and degree of the local regression models.

l2 <- loess(co2_avg ~ decimal_date, data = co2, 
            degree = 1, span = 0.05)
l3 <- loess(co2_avg ~ decimal_date, data = co2, 
            degree = 1, span = 0.25)
new_data = tibble(decimal_date = seq(min(co2$decimal_date), 
                                      max(co2$decimal_date),
                                     length = 300))
a2 <- augment(l2, newdata = new_data, se_fit = TRUE) 
a3 <- augment(l3, newdata = new_data, se_fit = TRUE) 
ggplot(a2, aes(x = decimal_date, y = .fitted)) + 
  geom_line(col="green") + 
  geom_ribbon(aes(ymin = .fitted - 2*.se.fit, 
                  ymax = .fitted + 2*.se.fit), 
              alpha = 0.20, fill="green") + 
  geom_point(aes(y = co2_avg), data = co2) +
  geom_line(data = a3, 
            col="blue") + 
  geom_ribbon(data = a3, 
              aes(ymin = .fitted - 2*.se.fit, 
                  ymax = .fitted + 2*.se.fit), 
              alpha = 0.20, fill="blue") + 
  my_theme

LOESS fits are slow and take a lot of memory compared to other methods (both time and storage requirements increase like the square of the number of points), so they are usually only used for small data sets.

17.3 Further reading

• LOESS (Wikipedia)
• Generalized additive models (Wikipedia)