Today's topics

• What are GAMs?

• How to fit GAMs in R with mgcv

• Examples

Logistics

• HTML Slide deck bit.ly/gam-efi-22 © Simpson (2022)

• RMarkdown bit.ly/gam-efi-git

HadCRUT4 time series

Linear Models

$$y_i\sim N({\mu }_i,{\sigma }^2)$$

$${\mu }_i={\beta }_0+{\beta }_1{x}_{1i}+{\beta }_2{x}_{2i}+\cdots +{\beta }_j{x}_{ji}$$

Assumptions

1. linear effects of covariates are good approximation of the true effects
2. conditional on the values of covariates, $y_i|X\sim N(0,{\sigma }^2)$
3. this implies all observations have the same variance
4. $y_i|X$ are independent

An additive model address the first of these

Is this linear?

Polynomials perhaps…

Polynomials perhaps…

We can keep on adding ever more powers of $x$ to the model — model selection problem

Runge phenomenon — oscillations at the edges of an interval — means simply moving to higher-order polynomials doesn't always improve accuracy

HadCRUT data set

## Load Data
tmpf <- tempfile()
curl_download("https://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/HadCRUT.4.6.0.0.annual_nh.txt", tmpf)
gtemp <- read.table(tmpf, colClasses = rep("numeric", 12))[, 1:2] # only want some of the variables
names(gtemp) <- c("Year", "Temperature")
gtemp <- as_tibble(gtemp)

File format

HadCRUT data set

gtemp

## # A tibble: 172 × 2
##     Year Temperature
##    <dbl>       <dbl>
##  1  1850          NA
##  2  1851          NA
##  3  1852          NA
##  4  1853          NA
##  5  1854          NA
##  6  1855          NA
##  7  1856          NA
##  8  1857          NA
##  9  1858          NA
## 10  1859          NA
## # … with 162 more rows

Fitting a GAM

library('mgcv')
m <- gam(Temperature ~ s(Year), data = gtemp, method = 'REML')
summary(m)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Temperature ~ s(Year)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.020523   0.009733  -2.109   0.0365 *
## ---
## 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) 7.838   8.65 145  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.88   Deviance explained = 88.5%
## -REML = -91.196  Scale est. = 0.016294  n = 172

Fitted GAM

Generalized Additive Models

Source: GAMs in R by Noam Ross

How is a GAM different?

In LM we model the mean of data as a sum of linear terms:

$$y_i={\beta }_0+\sum _j{\beta }_j{x}_{ji}+{ϵ}_i$$

A GAM is a sum of smooth functions or smooths

$$y_i={\beta }_0+\sum _j{s}_j\left(x_{ji}\right)+{ϵ}_i$$

where ${ϵ}_i\sim N(0,{\sigma }^2)$, $y_i\sim \text{Normal}$ (for now)

Call the above equation the linear predictor in both cases

Fitting a GAM in R

model <- gam(y ~ s(x1) + s(x2) + te(x3, x4), # formuala describing model
             data = my_data_frame,           # your data
             method = 'REML',                # or 'ML'
             family = gaussian)              # or something more exotic

s() terms are smooths of one or more variables

te() terms are the smooth equivalent of main effects + interactions

How did gam() know?

Wiggly things

Basis expansions

In the polynomial models we used a polynomial basis expansion of $x$

• $x^0=1$ — the model constant term
• $x^1=x$ — linear term
• $x^2$
• $x^3$
• …

Splines

Splines are functions composed of simpler functions

Simpler functions are basis functions & the set of basis functions is a basis

When we model using splines, each basis function $b_k$ has a coefficient ${\beta }_k$

Resultant spline is a the sum of these weighted basis functions, evaluated at the values of $x$

$$s\left(x\right)=\sum _{k=1}^K{\beta }_k{b}_k\left(x\right)$$

Splines formed from basis functions

Weight basis functions ⇨ spline

How do GAMs learn from data?

Maximise penalised log-likelihood ⇨ β

Wiggliness

$${\int }_R[f^{\prime \prime }{]}^{2}dx={\beta }^{T}S\beta =W$$

(Wiggliness is 100% the right mathy word)

We penalize wiggliness to avoid overfitting

Making wiggliness matter

$W$ measures wiggliness

(log) likelihood measures closeness to the data

We use a smoothing parameter $\lambda$ to define the trade-off, to find the spline coefficients $B_k$ that maximize the penalized log-likelihood

$$L_p=\log \left(\text{Likelihood}\right)-\lambda W$$

HadCRUT4 time series

Picking the right wiggliness

Maximum allowed wiggliness

We set basis complexity or "size" $k$

This is maximum wigglyness, can be thought of as number of small functions that make up a curve

Once smoothing is applied, curves have fewer effective degrees of freedom (EDF)

EDF < $k$

Maximum allowed wiggliness

$k$ must be large enough, the $\lambda$ penalty does the rest

Large enough — space of functions representable by the basis includes the true function or a close approximation to the tru function

Bigger $k$ increases computational cost

In mgcv, default $k$ values are arbitrary — after choosing the model terms, this is the key user choice

Must be checked! — gam.check()

GAM summary so far

1. GAMs give us a framework to model flexible nonlinear relationships

2. Use little functions (basis functions) to make big functions (smooths)

3. Use a penalty to trade off wiggliness/generality

4. Need to make sure your smooths are wiggly enough

A cornucopia of smooths

The type of smoother is controlled by the bs argument (think basis)

The default is a low-rank thin plate spline bs = 'tp'

Many others available

A bestiary of conditional distributions

A GAM is just a fancy GLM

Simon Wood & colleagues (2016) have extended the mgcv methods to some non-exponential family distributions

Smooth interactions

Two ways to fit smooth interactions

1. Bivariate (or higher order) thin plate splines
   • s(x, z, bs = 'tp')
   • Isotropic; single smoothness parameter for the smooth
   • Sensitive to scales of x and z
2. Tensor product smooths
   • Separate marginal basis for each smooth, separate smoothness parameters
   • Invariant to scales of x and z
   • Use for interactions when variables are in different units
   • te(x, z)

Tensor product smooths

There are multiple ways to build tensor products in mgcv

1. te(x, z)
2. t2(x, z)
3. s(x) + s(z) + ti(x, z)

te() is the most general form but not usable in gamm4::gamm4() or brms

t2() is an alternative implementation that does work in gamm4::gamm4() or brms

ti() fits pure smooth interactions; where the main effects of x and z have been removed from the basis

Factor smooth interactions

Two ways for factor smooth interactions

1. by variable smooths
   • entirely separate smooth function for each level of the factor
   • each has it's own smoothness parameter
   • centred (no group means) so include factor as a fixed effect
   • y ~ f + s(x, by = f)
2. bs = 'fs' basis
   • smooth function for each level of the function
   • share a common smoothness parameter
   • fully penalized; include group means
   • closer to random effects
   • y ~ s(x, f, bs = 'fs')

Random effects

When fitted with REML or ML, smooths can be viewed as just fancy random effects

Inverse is true too; random effects can be viewed as smooths

If you have simple random effects you can fit those in gam() and bam() without needing the more complex GAMM functions gamm() or gamm4::gamm4()

These two models are equivalent

m_nlme <- lme(travel ~ 1, data = Rail, ~ 1 | Rail, method = "REML") 
m_gam  <- gam(travel ~ s(Rail, bs = "re"), data = Rail, method = "REML")

Random effects

The random effect basis bs = 're' is not as computationally efficient as nlme or lme4 for fitting

• complex random effects terms, or
• random effects with many levels

Instead see gamm() and gamm4::gamm4()

• gamm() fits using lme()
• gamm4::gamm4() fits using lmer() or glmer()

For non Gaussian models use gamm4::gamm4()

Next steps

Read Simon Wood's book!

Lots more material on our ESA GAM Workshop site

https://noamross.github.io/mgcv-esa-workshop/

Noam Ross' free GAM Course

https://noamross.github.io/gams-in-r-course/

A couple of papers:

Also see my blog: www.fromthebottomoftheheap.net

Acknowledgments

Slides

• HTML Slide deck bit.ly/gam-efi-22 © Simpson (2022)
• RMarkdown bit.ly/gam-efi-git