Does a GAM really work better in this case 🤔?

The answer is yes. Yes it does. Here’s a plot showing it:

ggplot(climate_data, aes(x = year, y = temp_anomaly)) +
  geom_point() +
  geom_smooth(aes(color = "GAM smooth"), method = "gam", se = FALSE, formula = y ~ s(x)) +
  geom_smooth(aes(color = "Linear fit"), method = "lm", se = FALSE,  formula = y ~ x) +
  scale_color_manual(
    name = NULL,
    values = c("GAM smooth" = "blue", "Linear fit" = "red")
  ) +
  labs(x = "Year", y = "Temp Anomalies (C)") +
  labs(title = "Temp Anomalies vs Time",
       x = "Year",
       y = "Temp Anomalies (C)",
       caption = "Sources: NOAA/GML (Co2)") +
  theme_minimal()

I’ll admit these data don’t look particularly wiggly at a glance, but we can still see that the lm (red line) really doesn’t fit the data perfectly. Don’t believe me? Well let’s look at it another way.

m_lm <- gam(temp_anomaly ~ year, data = climate_data, method = "REML")
m_gam <- gam(temp_anomaly ~ s(year, k = 10, bs = "cr"), data = climate_data, method = "REML")
# Build table
model_table <- data.frame(
  Model      = c("Linear", "GAM"),
  LogLik     = c(as.numeric(logLik(m_lm)),  as.numeric(logLik(m_gam))),
  Parameters = c(attr(logLik(m_lm), "df"),  attr(logLik(m_gam), "df")),
  AIC        = c(AIC(m_lm), AIC(m_gam))
)
# pretty table with kable
model_table |>
  kable(digits = 2, 
        col.names = c("Model", "Log-Likelihood", "Parameters", "AIC")) |>
  kable_styling(bootstrap_options = c("striped", "hover"))

So above I am comparing the models using Akaike’s Information Criterion (AIC). The long and short of AIC is that it compares the log likelihood of your model to the number of parameters it took to estimate that model (and penalizes you for using more parameters). We can see that the log likelihood is higher for our GAM and that AIC is much lower (which is better). This blog post is not about AIC (although I think that doing one on AIC would be useful), but it should be noted that AIC is not a perfect metric. The metric has caused some controversy due to it’s misue in Ecology specifically $ ^{4,5} $.

When we fit a GAM we do not actually get ‘parameters’ in the traditional sense. Instead we get something called effective degrees of freedom (edf). So it’s misleading for me to lump these both into a column called ‘parameters’. Unlike a linear model where the number of parameters is fixed and whole (for example in the parameters cell for linear model we have a estimated an intercept, slope and residual variance so we have 3 degrees of freedom used up), a GAM’s smooth is penalized toward simplicity during fitting. The EDF reflects how much flexibility (wiggliness) the smooth actually used — an EDF of 1 means the smooth collapsed to a straight line, while higher values indicate more ‘wiggliness’. It’s called “effective” because it’s not a count of literal parameters but a continuous measure of model complexity that the data and penalty jointly determine. Because GAMs do not have real ‘parameters’ calculating AIC actually uses EDF in their stead.

“Wiggly” data: what’s the basis?

I have mentioned that GAMs work by fitting wiggly lines to the data. I want to build some intuition on how this works. Underneath the hood, when we fit a GAM we are fitting what are called basis functions to the data. A basis function is essentially a simple “shape” that is fit to your data. The final “smooth” function is weighted by minimizing the residuals between the data and the basis function and the complexity (how wiggly that basis function is) by mgcv automatically. We can actually see this action with our own data by controlling a parameter called $ k $ (AKA the number of knots or basis functions used to fit the final smooth) in MGCV. Let’s see this in action by playing around with the $ k $ parameter.

# generate predictions for each knot (k)
preds <- map(c(2, 5, 10), function(k) {
  m <- gam(temp_anomaly ~ s(year, k = k, bs = "cr"), data = climate_data, method = "REML")
  tibble(
    year  = climate_data$year,
    temp   = fitted(m),
    k     = paste0("k = ", k)
  )
}) %>% 
  bind_rows() %>%
  mutate(k = factor(k, levels = c("k = 2", "k = 5", "k = 10")))
# Plot
ggplot(climate_data, aes(x = year, y = temp_anomaly), alpha = 0.4) +
  geom_point() +
  geom_line(data = preds, aes(x = year, y = temp, color = k), linewidth = 0.8) +
  scale_color_viridis_d(name = "Model") +
  labs(x = "Year", x = "Temp Anomaly (C)") +
  theme_bw() +
  facet_wrap(~k)

Note: By default $ k = 10 $. Also notice that I have introduced the parameter bs. This stands for basis and is what we control to change the type of basis function that mgcv uses (essentially the type of basis determines the shape). The most important distinction is that some basis functions are ‘local’ (e.g. cr which stands for cubic regression spline used above) whereas others are ‘global’ (e.g., the default bs parameter tp which stands for thinplate). A local basis function takes k to be literal. It finds the ‘densest’ regions of data and assigns one of the k availble knots to those dense parts and selects the locations of k that minimize the distance between the function and the data. A global spline (again like the tp option) fits k basis functions to the data. Each of those basis functions span the entire dataset. While both methods end up fitting k basis functions the tp option (for example) will fit an $ n x n $ matrix (where $ n $ is the number of datapoints) and select the top k of those functions. So computationally the global option is much more expensive for large datasets.

m7 <- gam(temp_anomaly ~ s(year, k = 7, bs = "cr"), data = climate_data, method = "REML")

# let's start by looking at the basic summary output
summary(m7)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## temp_anomaly ~ s(year, k = 7, bs = "cr")
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.48909    0.01316   37.17   <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(year) 2.861  3.515 173.8  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.919   Deviance explained = 92.3%
## -REML = -43.402  Scale est. = 0.009524  n = 55

The key quantities to notice here are edf (again the effective degrees of freedom) and it’s associated p-value (note that this p-value does not have the same interpretation or formula as a p-value in a standard parameteric statistical test, but it does give some notion of the signficance of the shape of the relationship nonetheless). Remember that for this model we manually set $ k = 7 $. The fact that edf is not close to the $ k $ ceiling that we enforced (but still greater than 2) tells us two things. First, it is statistically true that the relationship between temp_anomaly and year is nonlinear. Second, the fact that edf is not close to \(k\) tells us that mgcv may have run out of ‘wiggly-ness’ to assign to the relationship. One thing you may notice is that unlike a typical model fit with lm (or some other parametric statistical model) we actually don’t get effect sizes here, meaning we get some information about the shape of the relationship between year and temp_anomaly but not any information about whether or not that relationship is positively or negatively related. To see whether the relationship is positive or negative we need to look at something called the partial effect. Also notice that the intercept in this case 0.49 is the average temp anomaly across the study years.

A partial effect is just the partial contribution of a predictor (i.e., year) to our response (i.e., temp_anomaly). Thankfully there are handy, programmable ways to see these partial effects clearly with a companion package to mgcv called gratia.

draw(m7)

So this plotstill looks pretty smooth and linear (and also the confidence bounds are pretty small), but regardless we can interpret it. The partial effect of year on temperature anomalies is early in the time period and positive later in the time period. This is consistent with the general trend of global warming where global temperatures have been rising (relatively steadily) year over year since the dawn of the industrial revolution.

Plot with gratia

gratia::draw(m_multi)

Notice that we get two side-by-side plots (one per smoothed term)…but these plots no longer seem so sensible. In fact, let’s check to see that this is the case by running a model with co2_mean and temp_anomaly.

m_co2_mean <- gam(temp_anomaly ~ s(co2_mean, bs = "cr"), data = climate_data)
gratia::draw(m_co2_mean)

OK so something very funky is going on with the multi-variate GAM partial effects above…Thankfully, I have a pretty strong hunch about what’s going on. To confirm that hunch, let’s take a look at one more function.

mgcv::concurvity(m_multi)

##                  para   s(year) s(co2_mean)
## worst    1.843496e-28 0.9999336   0.9999336
## observed 1.843496e-28 0.9998354   0.9997859
## estimate 1.843496e-28 0.9978072   0.7275050

The mgcv function concurvity measures the…well concurvity between two smoothed predictors. This is analogous to the idea of multicollinearity in parametric space, and what the output of the concurvity test is saying is…telling. Concurvity is a more generalized than that. It actually is a measure of how well s(year) can predict s(co2_mean) (and vice versa). So the first column is how well s(year) is predicted by s(co2_mean) and the second is how well s(co2_mean) is predicted by s(year). The rows correspond to the worst case (i.e., an upper limit on the concurvity), the observed case (based on the actual fitted values), and the estimate case (generally the most realistic and relevant). What we cansay is that all of these terms are highly concurve with one another and making it difficult to suss out true variation. In this case it would likely be wise to drop one of these covariates entirely.