# How to interpret and report nonlinear effects from Generalized Additive Models

By Nicholas Clark in rstats mgcv

May 2, 2024

## What are Generalized Additive Models (GAMs)?

Generalized Additive Models (GAMs) are flexible tools that replace one or more predictors in a Generalized Linear Model (GLM) with *smooth functions* of predictors. These are helpful for learning arbitrarily complex, nonlinear relationships between predictors and conditional responses without needing *a priori* expectations about the shapes of these relationships. Rather, they are learned using
*penalized smoothing splines*.

How do these work? The secret is a *basis expansion*, which in lay terms means that the covariate (time, in this example) is evaluated at a smaller set of *basis functions* designed to cover the range of the observed covariate values. Below is one particular type of basis, called a cubic regression basis.

Each basis function acts on its own local neighbourhood of the covariate. Once we have constructed the basis, we can estimate weights for each function. The weights allow the basis functions to have different impacts on the shape of the spline. They also give us a target for learning splines from data, as the weights act as the regression coefficients.

There are many more types of basis expansions that can be used to form penalized smooths, including multidimensional smooths, spatial smooths or even monotonic smooths.

### When should you use a Generalized Additive Model?

GAMs are incredibly useful for several reasons. Some of these include:

- You suspect there are nonlinear relationships among your response(s) and your predictor(s)
- You have no prior knowledge about what kind of shapes these relationships might take
- You are aware that high-order polynomials can capture nonlinear effects, but you know that these are extremely unwieldy and dangerous
- You’d like to learn these shapes from your data without severe risks of overfitting

But unfortunately, GAMs are often considered difficult to interpret (at least compared to simpler linear models). For better or worse,
we as scientists are often trained to look first at summary outputs from regression models and try to interpret the coefficients, their standard errors and their approximate *p*-values directly. I don’t agree with this practice, but it is so ingrained in tradition that it isn’t going away anytime soon. But this type of direct (dare I say, *lazy*?) interpretation is even more challenging when it comes to GAMs. This is because we can get summary statistics for tens or even hundreds of basis coefficients, none of which can be interpreted individually.

### How do we interpret nonlinear effects from Generalized Additive Models (GAMs)?

I’ve taught several workshops on GAMs and given many more seminars, and consistently this is the most common question I receive. In this post, I illustrate the challenges of smoothing spline interpretation, and I provide 3 pointers that you can follow to start understanding, interpreting and reporting nonlinear effects from GAMs. I will make use of the widely popular
the {`mgcv`

} R package for fitting GAMs to data. And I will show how the powerful
{`marginaleffects`

} R package can be used to easily generate plots, summaries and hypothesis tests to help us interpret these models.

Note however that this post emphasizes how to *interpret* fitted nonlinear functions. I will not talk about how to address statistical significance or how to select and build models. There are many other useful resources to guide these steps, most notably the help documentation in the {`mgcv`

} package (starting with `?mgcv::gam.models`

).

## Environment setup and the initial GAM model

This tutorial relies on the following packages:

```
library(mgcv) # Fit and interrogate GAMs
library(tidyverse) # Tidy and flexible data manipulation
library(marginaleffects) # Compute conditional and marginal effects
library(ggplot2) # Flexible plotting
library(patchwork) # Combining ggplot objects
```

I’ll also define a custom `ggplot2`

theme

```
theme_set(theme_classic(base_size = 12, base_family = 'serif') +
theme(axis.line.x.bottom = element_line(colour = "black",
size = 1),
axis.line.y.left = element_line(colour = "black",
size = 1),
panel.spacing = unit(0, 'lines'),
legend.margin = margin(0, 0, 0, -15)))
```

Now I load the `CO2`

data from the built-in datasets that are installed with R. These data contain measurements from an experiment on the cold tolerance of the grass species *Echinochloa crus-galli*, which you can read more about using `?CO2`

```
data(CO2, package = "datasets")
```

Inspect the data structure and dimensions

```
dplyr::glimpse(CO2)
```

```
## Rows: 84
## Columns: 5
## $ Plant <ord> Qn1, Qn1, Qn1, Qn1, Qn1, Qn1, Qn1, Qn2, Qn2, Qn2, Qn2, Qn2, …
## $ Type <fct> Quebec, Quebec, Quebec, Quebec, Quebec, Quebec, Quebec, Queb…
## $ Treatment <fct> nonchilled, nonchilled, nonchilled, nonchilled, nonchilled, …
## $ conc <dbl> 95, 175, 250, 350, 500, 675, 1000, 95, 175, 250, 350, 500, 6…
## $ uptake <dbl> 16.0, 30.4, 34.8, 37.2, 35.3, 39.2, 39.7, 13.6, 27.3, 37.1, …
```

Now some manipulations for modeling

```
plant <- CO2 |>
as_tibble() |>
rename(plant = Plant, type = Type, treatment = Treatment) |>
mutate(plant = factor(plant, ordered = FALSE))
```

I fit a Generalized Additive Model to these data, specifying a parametric interaction between `treatment`

and `type`

as well as random (hierarchical) intercepts for each `plant`

and using Carbon Dioxide uptake (`uptake`

) as the non-negative, real-valued response. The nonlinear part comes from a
*factor smooth interaction* that allows a nonlinear effect of Carbon Dioxide concentration (`conc`

) to vary for different levels of the cold treatment variable (`treatment`

). Given the nature of the response, a `Gamma()`

distribution is an appropriate observation model

```
model_1 <- gam(uptake ~ treatment * type +
s(plant, bs = "re") +
s(conc, by = treatment, k = 7),
data = plant,
method = "REML",
family = Gamma(link = "log"))
```

Look at the sumary of this model

```
summary(model_1)
```

```
##
## Family: Gamma
## Link function: log
##
## Formula:
## uptake ~ treatment * type + s(plant, bs = "re") + s(conc, by = treatment,
## k = 7)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.51609 0.06376 55.143 < 2e-16 ***
## treatmentchilled -0.11321 0.09018 -1.255 0.213987
## typeMississippi -0.31196 0.09018 -3.459 0.000978 ***
## treatmentchilled:typeMississippi -0.36041 0.12753 -2.826 0.006311 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(plant) 7.032 8.00 8.022 <2e-16 ***
## s(conc):treatmentnonchilled 5.193 5.68 83.888 <2e-16 ***
## s(conc):treatmentchilled 5.013 5.55 58.970 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.957 Deviance explained = 95.5%
## -REML = 239.08 Scale est. = 0.010327 n = 84
```

There seem to be lots of “*significant*” effects here, but how on Earth do we interpret these smooth terms? The summary gives no indication of the magnitude, direction or overall shape of these smooth effects. Nor do the coefficients, for that matter:

```
coef(model_1)
```

```
## (Intercept) treatmentchilled
## 3.516091294 -0.113205794
## typeMississippi treatmentchilled:typeMississippi
## -0.311956077 -0.360409395
## s(plant).1 s(plant).2
## -0.041226486 -0.015294663
## s(plant).3 s(plant).4
## 0.056521148 -0.040103411
## s(plant).5 s(plant).6
## 0.026931174 0.013172237
## s(plant).7 s(plant).8
## -0.055925539 0.049887951
## s(plant).9 s(plant).10
## 0.006037588 -0.215582180
## s(plant).11 s(plant).12
## 0.097118479 0.118463702
## s(conc):treatmentnonchilled.1 s(conc):treatmentnonchilled.2
## 5.232322869 1.133239546
## s(conc):treatmentnonchilled.3 s(conc):treatmentnonchilled.4
## -1.904830565 0.129187847
## s(conc):treatmentnonchilled.5 s(conc):treatmentnonchilled.6
## 0.236060128 1.202602610
## s(conc):treatmentchilled.1 s(conc):treatmentchilled.2
## 3.690591324 2.117499952
## s(conc):treatmentchilled.3 s(conc):treatmentchilled.4
## -2.367697749 0.063467378
## s(conc):treatmentchilled.5 s(conc):treatmentchilled.6
## 0.236541597 0.967747793
```

The coefficients labelled `s(conc):treatmentnonchilled`

and `s(conc):treatmentchilled`

are the basis function weights that control the shape of these splines. But without having a full understanding of what these basis functions look like and how they act collectively to form the spline, we are fairly lost for interpretation.

Now I will walk through 3 steps that you can use to help with interpreting and reporting nonlinear effects from GAMs. These steps are by no means exhaustive, but given the

## Step 1: Vizualising GAM outputs

Visuals are by-far our best go-to tool for understanding GAMs. This is why the {`mgcv`

} package has some simple but powerful plotting routines that we can apply directly to fitted `gam`

objects. See `?mgcv::plot.gam`

for more details about default plotting.

### Plotting GAM effects on the link scale

The most common visualisattion that users make when working with GAMs are
*partial effects* plots on the link scale. Basically, these plots show the individual component effect of a smooth function on the link scale, conditional on all other terms in the model being set to zero. Because most smooth terms in {`mgcv`

} are zero-centred to maintain identifiability, they are convenient to look at because we can immediately make some interpretations. For example, we can see the partial effect of the `conc`

smooth for the `nonchilled treatment`

by selecting the second smooth term in the model:

```
plot(model_1, select = 2, shade = TRUE)
abline(h = 0, lty = 'dashed')
```

This plot shows that, for the `nonchilled treatment`

, the effect of `conc`

on the linear predictor is mostly negative (below zero) for small values of `conc`

, but it quickly becomes positive (above zero) for intermediate values of `conc`

. It then begins to plateau at larger values of `conc`

. We can also look at the effect for the other level of `treatment`

```
plot(model_1, select = 3, shade = TRUE)
abline(h = 0, lty = 'dashed')
```

This plot looks broadly similar, but it is not identical to the `nonchilled treatment`

effect.

These plots are a good first step to interrogate GAMs, but they have several limitations:

- They only show what the model expects to see if
*all other terms are dropped*. These are not helpful if we want to aggregate over multiple smooths (i.e.*What is the average smooth effect of*) or in situations where covariates occur in multiple terms. For example, a formula such as:`conc`

?`y ~ s(rain) + s(temp) + s(year) + ti(rain, temp) + ti(rain, year) + ti(temp, year)`

is perfectly allowed in`gam()`

(and very useful in many contexts). But the effects of our predictors of interest are spread across multiple smooth functions and we cannot interpret these individually - These plots are shown on the link scale, which can often translate nonlinearly into effects on the scale of interest (the response scale. For example we have a log link function in this
`Gamma()`

observation model, so the partial effects are not the same as what we’d expect on the actual scale of`uptake`

- We cannot ask more targeted questions using these plots, such as
*How does the effect of*or`conc`

differ among treatment groups?*How much more uptake would we expect if*`conc`

increased from 150 to 300 mL/L?

For the remainder of this post, I’ll switch over to making plots and summaries using the elegant {`marginaleffects`

} package. It turns out that the functions available in this package make answering *all* of the above questions easy and seamless. Moreover, we could easily fit different kinds of models (i.e. GLMs, location-scale models, survival models etc…) to these same data and use *the exact same functions* for making model comparisons. To begin illustrating, I’ll tackle the first point above by computing and plotting the average smooth effect of `conc`

across all `treatment`

levels. Again this plot is made on the link scale, but it is no longer zero-centred because the intercepts are taken into account

```
plot_predictions(model_1, condition = 'conc',
type = 'link') +
labs(y = "Linear predictor (link scale)",
title = "Average smooth effect of concentration",
subtitle = "Aggregated across treatments and types")
```

Of course we can very readily split these plots by `treatment`

and `type`

:

```
plot_predictions(model_1, condition = c('conc', 'treatment', 'type'),
type = 'link') +
labs(y = "Linear predictor (link scale)",
title = "Average smooth effect of concentration",
subtitle = "Per treatment, per type")
```

There are many more possibilities for making simple, one-line-of-code targeted plots using `plot_predictions()`

. See `?marginaleffects::plot_predictions`

for guidance. But we can do much more than plot conditional predictions using this framework. For example,
we often want to understand how the *slope* of a fitted function changes across the range of a covariate. We can calculate the slopes (i.e. the first derivatives) of our smooth functions using the `plot_slopes()`

function, which conveniently has very similar arguments and structure to `plot_predictions()`

:

```
plot_slopes(model_1, variables = 'conc',
condition = c('conc', 'treatment'),
type = 'link') +
geom_hline(yintercept = 0, linetype = 'dashed') +
labs(y = "1st derivative of the linear predictor",
title = "Conditional slopes of the concentration effect",
subtitle = "Per treatment, per type")
```

This plot makes it more clear that the slope is positive up until we hit `conc`

values near 260, beyond which the function plateaus. Now we will tackle the second point above, which relates to probably the second-most common question users ask when working with GAMs in {`mgcv`

}:

**How do I plot smooth effects on the outcome scale?**

### Plotting GAM effects on the response scale

In fact this is dead-simple with {`marginaleffects`

}. We just change `type = 'response'`

and away we go! Of course we now have options to show the observed values on these plots, which is something we might often want to do

```
plot_predictions(model_1, condition = 'conc',
type = 'response', points = 0.5,
rug = TRUE) +
labs(y = "Expected response",
title = "Average smooth effect of concentration",
subtitle = "Aggregated across treatments and types")
```

Here is the average smooth of `conc`

once again, but this time on the scale of the response variable. The points and rugs show observed values, which makes it clear that a single effect of `conc`

would not be enough to capture the heterogeneity in the data. Splitting by `treatment`

and `type`

again makes things more clear

```
plot_predictions(model_1, condition = c('conc', 'treatment', 'type'),
type = 'response', points = 0.5,
rug = TRUE) +
labs(y = "Expected response",
title = "Average smooth effect of concentration",
subtitle = "Per treatment, per type")
```

These types of plots give us more ability to interrogate and criticize our fitted models. Some questions you can start to ask when trying to interpret or report functions from GAMs include:

- Does the function reach an aysmptote anywhere over its range?
- Is the function highly wiggly, or does it appear smooth?
- Does the function appear to have multiple modes, and does this make sense?
- Are there large regions over which you have few data points, and does the function’s uncertainty grow appropriately?
- Are there any points that look like outliers, which the function may be responding strongly to?

And in case you were wondering, yes you can also plot slopes on the scale of the response:

```
plot_slopes(model_1, variables = 'conc',
condition = c('conc', 'treatment', 'type'),
type = 'response') +
geom_hline(yintercept = 0, linetype = 'dashed') +
labs(y = "1st derivative of the expected response",
title = "Conditional slopes of the concentration effect",
subtitle = "Per treatment, per type")
```

By now you may be wondering *where* these predictions are coming from. This brings me to the next tool in your arsenal for interpreting GAMs, simulation

## Step 2: Simulating from fitted GAM models

In a related blog post,
I show how simulating data in R can be essential to understand models and their limitations. This is especially applicable to GAMs because this simulation involves producing the
*linear predictor (i.e. design) matrix*. This matrix is what we get when our covariates are evaluated at their respective basis functions.

### Using `type = 'lpmatrix'`

to understand what coefficients mean in a GAM model

Generating the linear predictor matrix, which I’ll refer to as `\(X_{lp}\)`

, is simple in {`mgcv`

}. We just use `predict.gam()`

and set `type = 'lpmatrix'`

. This can be done for *new* data or for the original data that were used to fit the model. Here is the latter:

```
Xp <- predict(model_1, type = 'lpmatrix')
dim(Xp)
```

```
## [1] 84 28
```

You can see there are 84 rows in the matrix, corresponding to the 84 observations in our data. But we have 28 columns, many of which represent the basis functions for the two smooth terms in the model

```
colnames(Xp)
```

```
## [1] "(Intercept)" "treatmentchilled"
## [3] "typeMississippi" "treatmentchilled:typeMississippi"
## [5] "s(plant).1" "s(plant).2"
## [7] "s(plant).3" "s(plant).4"
## [9] "s(plant).5" "s(plant).6"
## [11] "s(plant).7" "s(plant).8"
## [13] "s(plant).9" "s(plant).10"
## [15] "s(plant).11" "s(plant).12"
## [17] "s(conc):treatmentnonchilled.1" "s(conc):treatmentnonchilled.2"
## [19] "s(conc):treatmentnonchilled.3" "s(conc):treatmentnonchilled.4"
## [21] "s(conc):treatmentnonchilled.5" "s(conc):treatmentnonchilled.6"
## [23] "s(conc):treatmentchilled.1" "s(conc):treatmentchilled.2"
## [25] "s(conc):treatmentchilled.3" "s(conc):treatmentchilled.4"
## [27] "s(conc):treatmentchilled.5" "s(conc):treatmentchilled.6"
```

These correspond to the coefficients, `\(\beta\)`

, that we extracted earlier from the fitted model

```
beta <- coef(model_1)
all.equal(names(beta), colnames(Xp))
```

```
## [1] TRUE
```

If we use linear algebra to cross-multiply these coefficients with the design matrix `\((X_{lp}\beta)\)`

we get the predicted values on the link scale:

```
preds <- as.vector(Xp %*% beta)
```

Running these through the inverse link function (which is `exp()`

in the case of our log link) gives us the fitted values from the model

```
all.equal(fitted(model_1), exp(preds))
```

```
## [1] TRUE
```

As with the examples
in my other blog post about simulation in R, you can investigate uncertainties in these predictions by drawing `\(\beta\)`

coefficients from the model’s implied multivariate normal posterior distribution. But I won’t cover that again here.

This is all well and good, but it becomes much more useful when we supply *new* data to the `predict.gam()`

function. For example, what might we predict if we had a new measurement for plant `Qn1`

measured on `treatment = nonchilled`

in `type = Mississippi`

and with CO2 concentration `conc = 278`

? This could be a targeted scenario that we are interested in:

```
newXp <- predict(model_1, type = 'lpmatrix',
newdata = data.frame(plant = 'Qn1',
treatment = 'nonchilled',
type = 'Mississippi',
conc = 278))
dim(newXp)
```

```
## [1] 1 28
```

Here we have taken our relatively simple scenario and run it through all necessary basis evaluations to produce the design matrix. Generating expected `uptake`

from this scenario is then as simple as:

```
exp(newXp %*% beta)
```

```
## [,1]
## 1 27.56059
```

It might be difficult to appreciate the power of this approach, but the possibilities to explore targeted scenarios are endless once you master these steps. For example, we can use these steps to see what the underlying basis functions for the effect of `conc`

look like. Simply simulate a sequence of fake `conc`

values to cover the range of the covariate and run those through the `predict.gam()`

function (keeping the other covariates fixed to particular values):

```
conc_seq <- seq(from = min(plant$conc), max(plant$conc), length.out = 500)
newdat <- data.frame(conc = conc_seq,
plant = 'Qn1',
treatment = 'nonchilled',
type = 'Mississippi')
newXp <- predict(model_1, type = 'lpmatrix', newdata = newdat)
```

Find which coefficients belong to the first smooth of `conc`

```
conc_coefs <- model_1$smooth[[2]]$first.para:model_1$smooth[[2]]$last.para
conc_coefs
```

```
## [1] 17 18 19 20 21 22
```

And now set all cells in the `\(X_{lp}\)`

matrix that don’t correspond to these coefficients to zero

```
newXp_adj <- matrix(0, ncol = NCOL(newXp), nrow = NROW(newXp))
newXp_adj[,conc_coefs] <- newXp[, conc_coefs]
```

Generate the predictions on the link scale and plot the functions

```
plot_dat = do.call(rbind, lapply(seq_along(conc_coefs), function(basis){
data.frame(basis_func = paste0('bs', basis),
conc = conc_seq,
value = newXp_adj[, conc_coefs[basis]] * beta[conc_coefs[basis]])
}))
ggplot(plot_dat, aes(x = conc, y = value, col = basis_func)) +
geom_line(linewidth = 0.75) + theme(legend.position = 'none') +
labs(y = 'Basis function value')
```

Essentially these steps are what we need to interrogate GAMs. But as you can see, it can be tedious to create the `newdata`

object because we need to specify values for all of the covariates. Even the ones we aren’t interested in (such as `plant`

or `type`

in the above example). This is where the automation from {`marginaleffects`

} becomes so extraordinarily valuable. By using sophisticated rules to automatically set values for missing predictors, creating these `newdata`

scenarios is effortless. The magic happens using the `datagrid()`

function:

```
newdat <- datagrid(conc = conc_seq, model = model_1)
head(newdat)
```

```
## uptake treatment type plant conc
## 1 27.2131 nonchilled Quebec Qn1 95.00000
## 2 27.2131 nonchilled Quebec Qn1 96.81363
## 3 27.2131 nonchilled Quebec Qn1 98.62725
## 4 27.2131 nonchilled Quebec Qn1 100.44088
## 5 27.2131 nonchilled Quebec Qn1 102.25451
## 6 27.2131 nonchilled Quebec Qn1 104.06814
```

Here you can see that only the `conc`

variable changes across rows, while all other variables are fixed. You can use `?marginaleffects::datagrid`

to see what the default rules are for fixing these (and how you can change them if you wish). Combining {`mgcv`

}’s `predict.gam(..., type = 'lpmatrix')`

and {`marginaleffects`

}’s `datagrid(...)`

function gives you all the flexibility you need to ask targeted questions using fake data simulation. In fact, when you call `plot_predictions()`

or `plot_slopes()`

, behind the scenes is a clever run of calls to `datagrid()`

to set up the design matrices.

### Sidenote: `marginaleffects`

works for *many* model classes

Obviously {`marginaleffects`

} is incredibly useful for interrogating GAMs as it provides a clean, simple interface for understanding nonlinear effects. But in fact the power of this framework is much broader, as it can support a *huge* number of model classes in R (over 100 at the time of writing). This makes model comparison an absolute breeze. For example, let’s fit a GLM with some wacky polynomial interaction effects to the same data (Please note I’m not advocating for this model. No. It is a *bad* model)

```
model_2 <- glm(uptake ~ treatment * type +
poly(conc, 4) * treatment + plant,
data = plant,
family = Gamma(link = "log"))
```

And we can remake one of the plots of the concentration effect from the GAM that we made earlier

```
plot_predictions(model_1, condition = c('conc', 'treatment', 'type'),
type = 'response', points = 0.5,
rug = TRUE) +
labs(y = "Expected response",
title = "Average smooth effect of concentration",
subtitle = "Per treatment, per type")
```

We can produce the exact same plot, but now with the GLM. Notice how not a single character in the call to `plot_predictions()`

has to change apart from the model argument

```
plot_predictions(model_2, condition = c('conc', 'treatment', 'type'),
type = 'response', points = 0.5,
rug = TRUE) +
labs(y = "Expected response",
title = "Average smooth effect of concentration",
subtitle = "Per treatment, per type")
```

Ok but how do we go beyond these relatively “simple” plots to ask more direct questions from our models?

## Step 3: Making useful comparisons with GAMs

In their recent book
Regression and Other Stories, authors Gelman, Hill and Vehtari advocate that regression outputs should be interpreted as *comparisons*. This is helpful because most of our models are descriptive rather than emulations of some causal assumptions we might have about the data-generating process. But it also helps to give us ideas about how we can use regression models to compare different scenarios and ask the model what it thinks might be plausible.

### Difference smooths: how does a smooth effect vary among groups?

A recent post on Cross Validated on computing differences among smooths in factor-smooth interactions (which partly motivated this post) provides useful context to help showcase the power of {`marginaleffects`

} for making comparisons in GAMs. The user was interested in fitting the same type of GAM that I fit above (including a factor-smooth interaction) to understand how nonlinear effects vary among factor levels. The question is straightforward enough, but it turns out this can be challenging to calculate by hand. But fortunately we have the family of `comparisons()`

functions from {`marginaleffects`

} to help. These are designed to *compare* predictions among different scenarios to compute differences, calculate ratios or run any user-specified function (i.e. changes in log odds or elasticities, if we so wish). See `?marginaleffects::comparisons`

for more guidance.

As an example, we can calculate the average expected difference in `uptake`

values between treatments by aggregating over `type`

and `conc`

. This gives us the “global” effect of `treatment`

, irrespective of `type`

or `conc`

values and can be a useful summary measure

```
avg_comparisons(model_1, newdata = datagrid(conc = conc_seq,
treatment = unique,
type = unique),
variables = "treatment",
by = "type")
```

```
##
## Term Contrast type Estimate Std. Error
## treatment mean(chilled) - mean(nonchilled) Quebec -4.76 3.04
## treatment mean(chilled) - mean(nonchilled) Mississippi -10.71 2.20
## z Pr(>|z|) S 2.5 % 97.5 %
## -1.57 0.117 3.1 -10.7 1.19
## -4.87 <0.001 19.8 -15.0 -6.40
##
## Columns: term, contrast, type, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
## Type: response
```

Here we can see that, on the scale of the response, the average difference among treatments is stronger for `type = Mississippi`

than for `type = Quebec`

. This contrast comes with an estimated effect size and *p*-value, which you can read more about in `?marginaleffects::comparisons`

Alternatively, we can target the original poster’s estimand of interest: what is the average difference between treatments across a sequence of `conc`

values (i.e. what is the difference among the smooths?)

```
avg_comparisons(model_1, newdata = datagrid(conc = conc_seq,
treatment = unique,
type = unique),
variables = "treatment",
by = c("conc", "type"))
```

```
##
## Term Contrast conc type Estimate
## treatment mean(chilled) - mean(nonchilled) 95.0 Quebec 0.413
## treatment mean(chilled) - mean(nonchilled) 95.0 Mississippi -3.114
## treatment mean(chilled) - mean(nonchilled) 96.8 Quebec 0.373
## treatment mean(chilled) - mean(nonchilled) 96.8 Mississippi -3.183
## treatment mean(chilled) - mean(nonchilled) 98.6 Quebec 0.332
## Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 1.60 0.258 0.79644 0.3 -2.73 3.55
## 1.02 -3.057 0.00224 8.8 -5.11 -1.12
## 1.61 0.232 0.81662 0.3 -2.78 3.53
## 1.03 -3.102 0.00192 9.0 -5.19 -1.17
## 1.62 0.205 0.83731 0.3 -2.84 3.51
## --- 990 rows omitted. See ?print.marginaleffects ---
## treatment mean(chilled) - mean(nonchilled) 996.4 Mississippi -11.426
## treatment mean(chilled) - mean(nonchilled) 998.2 Quebec -4.436
## treatment mean(chilled) - mean(nonchilled) 998.2 Mississippi -11.427
## treatment mean(chilled) - mean(nonchilled) 1000.0 Quebec -4.435
## treatment mean(chilled) - mean(nonchilled) 1000.0 Mississippi -11.428
## Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2.75 -4.156 < 0.001 14.9 -16.81 -6.04
## 3.97 -1.119 0.26334 1.9 -12.21 3.34
## 2.76 -4.147 < 0.001 14.9 -16.83 -6.03
## 3.98 -1.115 0.26467 1.9 -12.23 3.36
## 2.76 -4.138 < 0.001 14.8 -16.84 -6.02
## Columns: term, contrast, conc, type, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
## Type: response
```

This gives us the expected difference among the two smooths for a broad sequence of `conc`

values. It is extremely useful because it already takes into account any variation in intercepts or other effects where `treatment`

or `conc`

might show up in the model. Of course looking at a returned `tibble`

with 1,000 rows is not very helpful, so we can plot these differences instead:

```
plot_comparisons(model_1, newdata = datagrid(conc = conc_seq,
treatment = unique,
type = unique),
variables = "treatment",
by = c("conc", "type"),
type = 'link') +
geom_hline(yintercept = 0, linetype = 'dashed') +
labs(y = "Estimated difference",
title = "Difference between treatment levels",
subtitle = "Chilled - nonchilled, per type")
```

We can also ask questions such as *At what values of conc does the nonlinear slope grow the fastest?*. This kind of question is extremely useful in time series modeling, for example if we wanted to know when a population trend was increasing (or decreasing) the fastest. Supplying a custom function to the

`comparisons()`

function achieves this for us:```
max_growth = function(hi, lo, x) {
dydx <- (hi - lo) / 1e-6
dydx_max <- max(dydx)
x[dydx == dydx_max][1]
}
comparisons(model_1,
newdata = datagrid(conc = conc_seq,
treatment = unique,
type = unique),
variables = list("conc" = 1e-6),
vcov = FALSE,
by = "treatment",
comparison = max_growth)
```

```
##
## Term Contrast treatment Estimate
## conc +1e-06 nonchilled 157
## conc +1e-06 chilled 151
##
## Columns: term, contrast, treatment, estimate
## Type: response
```

Here we can see that, for both treatments, the slope of the nonlinear effect appears to be increasing most rapidly near values of `conc = 155`

.

Or we can ask *At what values of conc does the slope differ significantly from zero?*. By once again aggregating over

`type`

and computing first derivatives across a range of `conc`

values, we can use the `hypotheses()`

function which, by default, asks if a quantity is significantly different from zero:```
hypotheses(slopes(model_1,
newdata = datagrid(conc = conc_seq,
treatment = unique,
type = unique),
variables = "conc",
by = c("conc", "treatment"),
type = 'link')) %>%
dplyr::filter(p.value <= 0.05)
```

```
##
## Term Contrast conc treatment Estimate Std. Error z Pr(>|z|) S
## conc mean(dY/dX) 95.0 nonchilled 0.00809 0.000876 9.23 <0.001 65.0
## conc mean(dY/dX) 95.0 chilled 0.00647 0.000829 7.80 <0.001 47.2
## conc mean(dY/dX) 96.8 nonchilled 0.00809 0.000875 9.24 <0.001 65.1
## conc mean(dY/dX) 96.8 chilled 0.00647 0.000829 7.81 <0.001 47.3
## conc mean(dY/dX) 98.6 nonchilled 0.00808 0.000874 9.25 <0.001 65.2
## 2.5 % 97.5 %
## 6.37e-03 0.00980
## 4.84e-03 0.00809
## 6.37e-03 0.00980
## 4.84e-03 0.00809
## 6.37e-03 0.00980
## --- 196 rows omitted. See ?print.marginaleffects ---
## conc mean(dY/dX) 278.2 nonchilled 0.00127 0.000503 2.53 0.0114 6.5
## conc mean(dY/dX) 280.0 nonchilled 0.00122 0.000511 2.39 0.0168 5.9
## conc mean(dY/dX) 281.8 nonchilled 0.00117 0.000521 2.25 0.0242 5.4
## conc mean(dY/dX) 283.6 nonchilled 0.00113 0.000535 2.10 0.0355 4.8
## conc mean(dY/dX) 285.4 nonchilled 0.00108 0.000547 1.97 0.0484 4.4
## 2.5 % 97.5 %
## 2.87e-04 0.00226
## 2.21e-04 0.00222
## 1.53e-04 0.00219
## 7.67e-05 0.00218
## 7.52e-06 0.00215
## Columns: term, contrast, conc, treatment, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
## Type: link
```

Sometimes we might want more specific hypotheses, and these are often readily achievable using {`marginaleffects`

}. See for example
this post on Cross Validated looking into biases of slope estimates from GAMs.

## How do I report effects from GAMs in scientific publications?

Hopefully that rundown gives you some helpful tools and tips to interrogate your GAMs. I haven’t dived into multidimensional smooths in this post, but in principle the workflow doesn’t need to change much. You just need to be specific about the types of comparisons you want to make and let {`marginaleffects`

} turn the inference crank. I’ll sum up by attempting to tackle another one of the most common questions I get asked about GAMs: how do I report these things? Here are a few points that I think are valuable for reporting nonlinear effects of GAMs in scientific communications:

- Avoid reliance on
*p*-values and statistical significance. Discretizing these effects into some arbitrary binary classification does nothing but throw out valuable information - Focus on how these effects influence expected change
*on the scale of the outcome*. We use GAMs and GLMs because we want to model real things in the world that aren’t Gaussian, so we use link functions that are often nonlinear. Interpreting effects on the scale of the linear predictor is hard for everyone - Use simulations to understand what your model thinks is plausible, and fully describe (and justify) your simulation strategy
- Make comparisons against different models (i.e. polynomial, linear) to get a sense for how sensitive your conclusions are to the choice of functional form

Thanks very much for your time!

## Further reading

The following papers and resources offer useful material for working with and interrogating GAMs

Christensen, E. M., Simpson, G. L., & Ernest, S. M. (2019).
Established rodent community delays recovery of dominant competitor following experimental disturbance. *Proceedings of the Royal Society B*, 286(1917), 20192269.

Clark, N. J., Ernest, S. K. M., Senyondo, H., Simonis, J. L., White, E. P., Yennis, G. M., & Karunarathna, K. A. N. K. (2023).
Multi-species dependencies improve forecasts of population dynamics in a long-term monitoring study. *Biorxiv*, doi:
https://doi.org/10.32942/X2TS34.

Karunarathna, K. A. N. K., Wells, K., & Clark, N. J. (2024).
Modelling nonlinear responses of a desert rodent species to environmental change with hierarchical dynamic generalized additive models. *Ecological Modelling*, 490, 110648.

Pedersen, E. J., Miller, D. L., Simpson, G. L., & Ross, N. (2019).
Hierarchical generalized additive models in ecology: an introduction with mgcv. *PeerJ*, 7, e6876.

Ross, N. (2019). GAMs in R: a free interactive course using mgcv

Simpson, G. L. (2018).
Modelling palaeoecological time series using generalised additive models. *Frontiers in Ecology and Evolution*, 6, 149.