# Interpretation of effect of factor-by-smooths in GAM

Given the following data and model, we'd like to assess if treatment affects the effect of conc on uptake.

library(gratia)
library(tidyverse)
library(mgcv)

# Data
data(CO2, package = "datasets")
plant <- CO2 |>
as_tibble() |>
rename(plant = Plant, type = Type, treatment = Treatment) |>
mutate(plant = factor(plant, ordered = FALSE),
uptake = ifelse(treatment == "chilled" & type == "Quebec", uptake+10, uptake))

# GAM
m_plant <- gam(uptake ~ treatment * type +
s(conc, by = treatment, k = 6) +
s(plant, bs = "re"),
data = plant,
method = "REML",
familly = Gamma(link = "log"))

# Fitted values
ds <- data_slice(m_plant,
treatment = levels(treatment),
type = levels(type),
conc = evenly(conc, n = 100))

fitted <-  fitted_values(m_plant,
data = ds,
scale = "response",
exclude = "s(plant)")

# Smooth estimates
se <- smooth_estimates(m_plant, "s(conc)",
partial_match = TRUE)


By predicting uptake, we can see that at given conc, uptake is higher for chilled plants for type = Mississippi and higher for nonchilled plants for type = Quebec. Can we still draw a general conclusion concerning the effect of treatment on the effect of conc on uptake, regardless of type, based on the smooth estimates? For instance, the plot of the smooth estimates shows that at e.g. conc = 500, the effect on uptake is stronger for nonchilled plants. Does such difference in intercept/position along the y-axis have any meaning when considering the global smooth effect? Alternatively, is it possible to compare the smooth effect for the types seperately (other than predicting values)?

Edit:

I think my question is related to this issue/discussion on github. Maybe I am just looking for something like smooth_estimates(..., group_means = TRUE). I tried to calculate that based on the difference smooth..

# Difference smooth with group means
ds_gm <- difference_smooths(m_plant, select = "s(conc)",
partial_match = TRUE,
group_means = TRUE) %>%
select(.diff, conc)

# Filter smooth estimate for nonchilled
se2 <- se %>%
filter(treatment == "nonchilled") %>%
select(.estimate, conc)

# Calculate smooth estimate for chilled from difference
left_join(ds_gm, se2) %>%
mutate(.estimate_chilled = .estimate - .diff) %>%
ggplot(.) +
geom_line(aes(x = conc, y = .estimate)) +
geom_line(aes(x = conc, y = .estimate_chilled), colour = "red")


Would this allow to conclude that the effect of conc on uptake is always stronger for chilled than nonchilled plants (confidence intervals left aside)?

It sounds like you are interested in the statistical estimation of differences between the two smooths formed by the s(conc, by = treatment, k = 6) term. If this is what you are interested in, your can compute differences between (pairs of) smooths using gratia::difference_smooths() but you can also do this yourself as it is based on a particular set of contrasts forming a specific linear combination of the model coefficients.

library("gratia")
library("mgcv")
library("tibble")
library("dplyr")

data(CO2, package = "datasets")
plant <- CO2 |>
as_tibble() |>
rename(plant = Plant, type = Type, treatment = Treatment) |>
mutate(plant = factor(plant, ordered = FALSE),
uptake = ifelse(treatment == "chilled" & type == "Quebec", uptake+10, uptake))

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

diffs <- difference_smooths(m_plant, select = "s(conc)")


With the dev version of {gratia} on GitHub this produces

> diffs
# A tibble: 100 × 9
.smooth .by       .level_1  .level_2   .diff    .se .lower_ci .upper_ci  conc
<chr>   <chr>     <chr>     <chr>      <dbl>  <dbl>     <dbl>     <dbl> <dbl>
1 s(conc) treatment nonchill… chilled  -0.262  0.0353    -0.331 -0.193      95
2 s(conc) treatment nonchill… chilled  -0.236  0.0310    -0.297 -0.175     104.
3 s(conc) treatment nonchill… chilled  -0.211  0.0273    -0.264 -0.157     113.
4 s(conc) treatment nonchill… chilled  -0.185  0.0246    -0.233 -0.137     122.
5 s(conc) treatment nonchill… chilled  -0.161  0.0229    -0.205 -0.116     132.
6 s(conc) treatment nonchill… chilled  -0.136  0.0223    -0.180 -0.0927    141.
7 s(conc) treatment nonchill… chilled  -0.113  0.0227    -0.158 -0.0688    150.
8 s(conc) treatment nonchill… chilled  -0.0911 0.0237    -0.137 -0.0447    159.
9 s(conc) treatment nonchill… chilled  -0.0700 0.0248    -0.119 -0.0213    168.
10 s(conc) treatment nonchill… chilled  -0.0503 0.0259    -0.101  0.000551  177.
# ℹ 90 more rows
# ℹ Use print(n = ...) to see more rows


(The CRAN version will produce something similar but the select argument is named something else and the variables in the returned object have different names.)

We can visualise this with the draw() method

diffs |> draw()


How do we interpret this? We can see what's happening via a few steps (that difference_smooths() doesn't actually do directly. First we create some values of conc and treatment to compute predicted values for:

ds <- data_slice(m_plant, conc = evenly(conc),
treatment = evenly(treatment))


Note this also creates values for plant and type (constant) but we'll remove those effects from the fitted values in the next step, using terms to choose just the effects of the two smooths created by the factor by smooth:

fv <- fitted_values(m_plant, data = ds, scale = "link",
terms = smooths(m_plant)[1:2])


Note we are doing this on the link (log) scale as that is what is shown on the plot you were discussing.

Now we can compute a difference using some dplyr foo (but we don't get any uncertainty estimates doing it this waY)

fv |>
group_by(conc) |>
summarise(.diff = -diff(.fitted))


This produces

# A tibble: 100 × 2
conc   .diff
<dbl>   <dbl>
1   95  -0.262
2  104. -0.236
3  113. -0.211
4  122. -0.185
5  132. -0.161
6  141. -0.136
7  150. -0.113
8  159. -0.0911
9  168. -0.0700
10  177. -0.0503
# ℹ 90 more rows
# ℹ Use print(n = ...) to see more rows


From this we see that chilled starts off having higher uptake than nonchilled, but that difference decreases quickly up to conc ~ 200, after which nonchilled leads to larger uptake values than chilled up to about conc 600, after which the difference decreases towards 0. Given this, we see that the effect of conc on uptake changes more rapidly as we increase conc up to about conc = 300, and after that we see similar changes in uptake as we increase conc in the two levels of the treatment, with nonchilled leading to larger uptake at a given conc. All of this is really in terms of $$\log(\text{uptake})$$ of course.

What we just looked at is the difference in the effects of conc on uptake between treatments excluding the treatment mean effects. So we are focusing on how the effect on uptake varies with conc alone, regardless of the constant offset between the two treatment levels. In other words we focus on differences in the shape of the smooth effects on the log scale.

If you're interested in the group means being included in the comparison, you can do that with group_means = TRUE:

difference_smooths(m_plant, select = "s(conc)", group_means = TRUE)


Does such difference in intercept/position along the y-axis have any meaning when considering the global smooth effect?

Yes, it is the magnitude of the effect of conc on the response centred on the model constant term (so centred about the mean for treatment == "nonchilled" & type == "Quebec")

However, your language suggests that you shouldn't be doing any of this on the log scale. The way you describe the problem it seems like you are interested in the actual response of uptake as a function of conc for the two treatment levels. In which case you should be using the response scale for this, excluding the effects of model terms you are not interested in:

fv2 <- fitted_values(m_plant, data = ds,
terms = c("treatmentchilled", smooths(m_plant)[1:2]))

fv2 |>
ggplot(aes(x = conc, y = .fitted, group = treatment)) +
geom_ribbon(aes(ymin = .lower_ci, ymax = .upper_ci, fill = treatment),
alpha = 0.2) +
geom_line(aes(colour = treatment))


This plot is showing the estimated effects of conc on uptake plus the difference between chilled and nonchilled levels. We've excluded the type effects here, but we are missing the average uptake to scale the y axis in terms of observed uptake values.

If you want a difference on this response scale, then the easiest way to get it is to use {marginaleffects} as Nick showed. Use comparisons() instead of avg_comparisons() if you want actual comparisons for specified data values, and you'll need to figure out what datagrid to use if you want evenly-spaced values over conc for plotting.

But I actually think you should be comparing treatment effects within type, as what does the average across types mean?

You can approach this with {gratia} using fitted_samples() to do posterior sampling, and then difference the fitted values for the comparisons of interests for each posterior draw. ({gratia} doesn't do the delta method to get confidence intervals on the response scale for difference_smooths-like operations. Instead we can do it with brute force via posterior sampling.)

First, set up the data we want to generate fitted values for

ds <- data_slice(m_plant, conc = evenly(conc), type = evenly(type),
treatment = evenly(treatment)) |>
mutate(.row = row_number())


Now we take a large set of posterior draws

fs <- fitted_samples(m_plant, data = ds, n = 50000, seed = 42)


add on the variables from ds

fs <- fs |> left_join(ds, by = join_by(.row == .row))


Now we compute differences between the two treatment levels and summarise the posterior distribution of those differences with the 0.25, 0.5, and 0.975 probability quantiles

library("tidyr")

quantile_fun <- function(x, probs = c(0.025, 0.5, 0.975), ...) {
tibble::tibble(
.value = quantile(x, probs = probs, ...),
.q = probs * 100
)
}

diffs <- fs |>
pivot_wider(id_cols = c("type", "conc", ".draw"),
names_from = "treatment", values_from = ".fitted") |>
mutate(.diff = nonchilled - chilled) |>
group_by(type, conc) |>
reframe(quantile_fun(.diff)) |>
pivot_wider(
id_cols = c(type, conc), names_from = .q, values_from = .value,
names_prefix = ".q"
)


Finally we plot:

diffs |>
ggplot(aes(x = conc, y = .q50, group = type)) +
geom_ribbon(aes(ymin = .q2.5, ymax = .q97.5, fill = type), alpha = 0.2) +
geom_line(aes(colour = type)) +
labs(y = "Estimated difference",
title = "Difference between treatment levels",
subtitle = "Comparison: nonchilled - chilled")


yielding

which matches what is shown in your faceted fitted values plots in the original post (just differencing the two curves per facet).

What I think this illustrates is that we have to be very clear about what we want to compare when working with GAMs, especially when we are working with non-identity link functions. Are we interested in differences of shape only? Differences among which factor levels? On the response scale or the link scale?

To close, how does my brute force method compare with the elegance of marginaleffects::comparisons()?

meff <- comparisons(m_plant, newdata = ds, variables = "treatment",
by = c("conc", "type"))

meff |> # note I'm negating the values because the comparison is opposite mine
ggplot(aes(x = conc, y = -estimate, group = type)) +
geom_ribbon(aes(ymin = -conf.low, ymax = -conf.high, fill = type),
alpha = 0.2) +
geom_line(aes(colour = type))


Which is pretty similar to the one I produced via brute force sampling, and that's reassuring.

Hopefully this gives you examples of the different ways you can compare your smooths in the GAM, and among those is the answer you wanted.

• Thanks, that is very helpful and I think I am slowly getting convinced that I have to work with fitted values for my purpose. Yet, something like smooth_estimates(..., group_means = TRUE), would that actually be useful at all (to be able to compare the effect of treatment at the scale of the effect of conc)? Commented Mar 20 at 13:03
• Adding group means would actually require some careful rethinking because of how it would interact with overall_uncertainty, and besides, these are partial effect plots and meant for comparison of the estimated smooth functions, not linear combinations of the model terms. With fitted-values() or the marginaleffects package you can achieve what you want (by specifying the terms or exclude arguments (only one of) to choose which terms to include or exclude respectively from the model predictions. Just ask for things on the link scale. Commented Mar 20 at 13:54

As you can see from Gavin's reply, this can be done with your {gratia} outputs. Hwer I briefly show how you can tackle this question with {marginaleffects}. I first compute the average (marginal) contrast for the two levels of treatment on the link scale (aggregated across values of conc). I then show how you can do this for a sequence of conc values to effectively ask whether one treatment consistently gives larger linear predictor values than another. Hopefully that gets at what your're interested in, but if not there are many other types of comparisons you can make, including specific hypothesis tests. See the {'marginaleffects'} hypothesis vignette for details.

library(tidyverse)
library(mgcv)
library(marginaleffects)
library(ggplot2)

# Data
data(CO2, package = "datasets")
plant <- CO2 |>
as_tibble() |>
rename(plant = Plant, type = Type, treatment = Treatment) |>
mutate(plant = factor(plant, ordered = FALSE),
uptake = ifelse(treatment == "chilled" & type == "Quebec",
uptake+10, uptake))

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

# Conditional smooths on the link scale (these are what we want to difference)
plot_predictions(m_plant,
condition = c('conc', 'treatment',
'type'),
type = 'link') +
theme_bw()



# Average difference (again on the link scale) between
# treatments (aggregating over 'type' and 'conc' values)
avg_comparisons(m_plant,
variables = "treatment",
by = "type", type = 'link')
#>
#>       Term                         Contrast        type Estimate Std. Error
#>  treatment mean(chilled) - mean(nonchilled) Quebec         0.186     0.0939
#>  treatment mean(chilled) - mean(nonchilled) Mississippi   -0.477     0.0939
#>      z Pr(>|z|)    S    2.5 % 97.5 %
#>   1.98   0.0479  4.4  0.00168  0.370
#>  -5.07   <0.001 21.3 -0.66071 -0.292
#>
#> Columns: term, contrast, type, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted

# Our estimand of interest: average differencea (again on the link scale) between
# treatments, aggregated over 'type' values for a sequence
# of 'conc' values
avg_comparisons(m_plant,
variables = "treatment",
by = c("conc", "type"),
#>
#>       Term                         Contrast conc        type Estimate
#>  treatment mean(chilled) - mean(nonchilled)   95 Quebec         0.448
#>  treatment mean(chilled) - mean(nonchilled)   95 Mississippi   -0.215
#>  treatment mean(chilled) - mean(nonchilled)  175 Quebec         0.241
#>  treatment mean(chilled) - mean(nonchilled)  175 Mississippi   -0.422
#>  treatment mean(chilled) - mean(nonchilled)  250 Quebec         0.131
#>  treatment mean(chilled) - mean(nonchilled)  250 Mississippi   -0.532
#>  treatment mean(chilled) - mean(nonchilled)  350 Quebec         0.102
#>  treatment mean(chilled) - mean(nonchilled)  350 Mississippi   -0.560
#>  treatment mean(chilled) - mean(nonchilled)  500 Quebec         0.101
#>  treatment mean(chilled) - mean(nonchilled)  500 Mississippi   -0.562
#>  treatment mean(chilled) - mean(nonchilled)  675 Quebec         0.134
#>  treatment mean(chilled) - mean(nonchilled)  675 Mississippi   -0.529
#>  treatment mean(chilled) - mean(nonchilled) 1000 Quebec         0.145
#>  treatment mean(chilled) - mean(nonchilled) 1000 Mississippi   -0.517
#>  Std. Error     z Pr(>|z|)    S   2.5 % 97.5 %
#>      0.1004  4.46   <0.001 16.9  0.2509  0.644
#>      0.1004 -2.14   0.0324  4.9 -0.4115 -0.018
#>      0.0974  2.47   0.0134  6.2  0.0500  0.432
#>      0.0974 -4.33   <0.001 16.0 -0.6124 -0.231
#>      0.0971  1.35   0.1785  2.5 -0.0597  0.321
#>      0.0971 -5.47   <0.001 24.4 -0.7221 -0.341
#>      0.0996  1.02   0.3062  1.7 -0.0933  0.297
#>      0.0996 -5.62   <0.001 25.7 -0.7557 -0.365
#>      0.1005  1.00   0.3169  1.7 -0.0964  0.297
#>      0.1005 -5.59   <0.001 25.4 -0.7588 -0.365
#>      0.1007  1.33   0.1847  2.4 -0.0638  0.331
#>      0.1007 -5.25   <0.001 22.7 -0.7262 -0.331
#>      0.1008  1.44   0.1492  2.7 -0.0522  0.343
#>      0.1008 -5.13   <0.001 21.7 -0.7146 -0.319
#>
#> Columns: term, contrast, conc, type, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted

# Plot these differences and add some labels
plot_comparisons(m_plant,
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") +
theme_bw()


Created on 2024-03-22 with reprex v2.1.0

• Thanks, unfortunately that does not answer my question. I am looking for how to interpret the global smooth effect of treatment (unconditional of type) OR comparing the smooth per type (but with the partial effect on the y-axis, not based on fitted values of uptake). Commented Mar 20 at 6:19
• What I provided gives you the average (global) effect of treatment, irrespective of type. But if you want to make comparisons per level of type, and you want this on the partial effect (i.e. link) scale, use the by argument: avg_comparisons(m_plant, variables = "treatment", by = "type", type = 'link'). I didn't originally catch the typo in your original call to the family, where you use: familly = Gamma(link = "log"). This is being ignored so you are actually getting a Gaussian model with link = 'identity'. Fixing this will produce more meaningful comparisons. Commented Mar 20 at 8:36
• Thanks! Still not sure if we are talking about the same aim though. I added some thoughts to my question, maybe that clarifies what I intend to do. Commented Mar 20 at 9:37
• Ah ok, I think I get it now. I've edited my response to show how you can target these questions (differences among smooths, with uncertainties accounted for). Hope that clarifies, but if not I'd suggest you read over the immense doc for the {marginaleffects} Commented Mar 20 at 10:24
• Thanks for all your effort, I'll definitely look into the marginaleffects package. Commented Mar 20 at 13:52