Gavin Simpson and colleagues have produced a lot of useful teaching material on GAMs that I’ve made use of here. For example:

Gavin’s blog

Webinar 1

Webinar 2

• Understand some diagnostics

• Understand more about model selection techniques with GAMs

• Become familiar with some of the flexibility of GAMs

A generalized additive model is an extension of the GLM with a predictor that is a sum of smoothed (and possibly unsmoothed) terms:

\[ y_i = \beta_0 + \sum_j s_j(x_{j,i}) + \sum_k \beta_k x_{k,i} + \epsilon_i \]

And with splines, our smoothed terms are a sum of weighted basis functions:

\[ s(x)=\sum_k \beta_k b_k(x) \]

Ask yourself:

• data/error distribution?

• predictors?

• type of smoother(s)?

Check candidate models via diagnostics

• Effect of choice of \(k\) (number of knots)?

Chicago air quality and death rate (Peng and Welty (2004); {NMMAPSdata})

death: total deaths per day

pm25median: median particles < 2.5 mg per cubic m

o3median: ozone in parts per billion

tmpd: temperature in Fahrenheit

## get data
data(chicago, package = "gamair")

## fit gam
fit_gam_air <- mgcv::gam(death ~ s(pm25median) + 
                           s(o3median) + 
                           s(tmpd), 
                         family = "poisson", 
                         data = chicago, method = "REML")

Smoothed Terms

##                    edf   Ref.df     Chi.sq      p-value
## s(pm25median) 3.518573 4.440907  20.860574 0.0006254809
## s(o3median)   3.430259 4.305665   6.373377 0.2230510318
## s(tmpd)       4.801666 5.878771 244.698414 0.0000000000

• Effective degrees of freedom (edf) - degree of non-linearity, but with upper bound set by \(k\)

• Reference degrees of freedom (Ref.df) - used for calculating test statistic

• p-value - null hypothesis is a linear effect (this can include a flat linear effect, i.e., no effect)

  • Biased low because they don’t account for need to estimate \(\lambda\)

  • Works well asymptotically

  • Doesn’t work for unpenalized smooths (random effects)

• Check whether \(k\) is adequate - can set it very large to start, but this increases computation time

• Do some iterative checking with larger values of \(k\) for multiple parameters, effective degrees of freedom should stabilize

## check model fit
gam.check(fit_gam_air) 

## fit with k = 20 for all
fit_gam_air <- mgcv::gam(death ~ s(pm25median, k = 20) + 
                           s(o3median, k = 20) + 
                           s(tmpd, k = 20), 
                         family = "poisson", 
                         data = chicago, method = "REML")

## check p-values
gam.check(fit_gam_air) 

## 
## Method: REML   Optimizer: outer newton
## full convergence after 3 iterations.
## Gradient range [-1.827205e-06,5.041726e-07]
## (score 2867.819 & scale 1).
## Hessian positive definite, eigenvalue range [0.1852984,1.155305].
## Model rank =  58 / 58 
## 
## 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(pm25median) 19.00  3.57    0.99    0.32  
## s(o3median)   19.00  3.59    0.97    0.24  
## s(tmpd)       19.00  5.01    0.94    0.06 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• When we select \(\lambda\) we are only penalizing wiggliness

• However, this only penalizes functions of \(X\) with a second derivative > 0

• Thus, linear relationships can’t be removed from the model this way

Use select = TRUE to apply penalties to the linear terms as well

## apply penalties to linear terms 
fit_gam_air <- mgcv::gam(death ~ s(pm25median, k = 20) + 
                           s(o3median, k = 20) + 
                           s(tmpd, k = 20), 
                         select = TRUE,
                         family = "poisson",
                         data = chicago, method = "REML")

Check to see if edf < 1

## check edf values
summary(fit_gam_air)$s.table

##                    edf   Ref.df     Chi.sq     p-value
## s(pm25median) 3.567278 4.543784  20.830454 0.000704636
## s(o3median)   3.591463 4.540236   7.006517 0.201616686
## s(tmpd)       5.006420 6.297518 243.866273 0.000000000

We have lots of options when fitting GAMs

so we have to give very careful thought to our goal(s)

We have lots of options when fitting GAMs

so we have to give very careful thought to our goal(s)

Just because we can fit a GAM doesn’t mean it will be interpretable or useful!

Cubic regression basis because have a lot of wiggliness in the data

Let’s make \(k\) really large then

## fit GAM with cubic regression basis
fit_co2_cubic <- gam(co2 ~ s(c.month, k = 300, bs = "cr"),
                     method = 'REML',
                     data = co2s)

Let’s now distinguish within-year variation from the long-term trend

• cyclic spline on month (1:12)

• cubic regression spline on c.month (1:507)

## fit GAM with 
b2 <- gam(co2 ~ s(month, bs = "cc") + 
            s(c.month, bs = "cr", k = 300), 
          knots = list(month = c(0.5, 12.5)),
          method = 'REML',
          data = co2s)

Building GAMs requires more choices than building GLMs

It requires thinking hard about the system and the question

main effects & interactions te(x1, x2)

• separate smoothness parameter and separate marginal basis

main effects & interactions te(x1, x2)

• separate smoothness parameter and separate marginal basis

interactions only ti(x1, x2)

• when the main effect are specified elsewhere

main effects & interactions te(x1, x2)

• separate smoothness parameter and separate marginal basis

interactions only ti(x1, x2)

• when the main effect are specified elsewhere

bivariate thin plate splines s(x1, x2, bs = 'tp')

• single smoothness parameter, only use when units are the same (eg, over space)

main effects & interactions te(x1, x2)

• separate smoothness parameter and separate marginal basis

interactions only ti(x1, x2)

• when the main effect are specified elsewhere

bivariate thin plate splines s(x1, x2, bs = 'tp')

• single smoothness parameter, only use when units are the same (eg, over space)

smooth by different levels of a factor f + s(x1, by = f)

• be sure to include the factor as a fixed effect

• akin to random regression slopes

We can also fit random effects

• mgcv::gam(...,bs = “re”)`

• mgcv::gamm()

• gamm4::gamm4()

Thus, these models are equivalent:

• lme(y ~ 1 + (1|f), method = "REML")

• gam(y ~ s(f, bs = "re"), method = "REML")

• GAMs using splines offer tremendous flexibility

• There is some tradeoff between flexibility and ease of understanding

• Many of our relationships in ecology aren’t linear

• What are some examples?