I've been exploring a number of tools for forecasting, and have found Generalized Additive Models (GAMs) to have the most potential for this purpose. GAMs are great! They allow for complex models to be specified very succinctly. However, that same succinctness is causing me some confusion, specifically in regard to how GAMs conceive of interaction terms and covariates.

Consider an example data set (reproducible code at end of post) in which y is a monotonic function perturbed by a couple of gaussians, plus some noise:

enter image description here

The data set has a few predictor variables:

  • x: The index of the data (1-100).
  • w: A secondary feature that marks out the sections of y where a gaussian is present. w has values of 1-20 where x is between 11 and 30, and 51 to 70. Otherwise, w is 0.
  • w2: w + 1, so that there are no 0 values.

R's mgcv package makes it easy to specify a number of possible models for these data:

enter image description here

Models 1 and 2 are fairly intuitive. Predicting y only from the index value in x at default smoothness produces something vaguely correct, but too smooth. Predicting y only from w results in a model of the "average gaussian" present in y, and no "awareness" of the other data points, all of which have a w value of 0.

Model 3 uses both x and w as 1D smooths, producing a nice fit. Model 4 uses x and w in a 2D smooth, also giving a nice fit. These two models are very similar, though not identical.

Model 5 models x "by" w. Model 6 does the opposite. mgcv's documentation states that "the by argument ensures that the smooth function gets multiplied by [the covariate given in the 'by' argument]". So shouldn't Models 5 & 6 be equivalent?

Models 7 and 8 use one of the predictors as a linear term. These make intuitive sense to me, as they are simply doing what a GLM would do with these predictors, and then adding the effect to the rest of the model.

Lastly, Model 9 is the same as Model 5, except that x is smoothed "by" w2 (which is w + 1). What's strange to me here is that the absence of zeros in w2 produces a remarkably different effect in the "by" interaction.

So, my questions are these:

  • What is the difference between the specifications in Models 3 and 4? Is there some other example that would draw out the difference more clearly?
  • What, exactly, is "by" doing here? Much of what I've read in Wood's book and this website suggests that "by" produces a multiplicative effect, but I'm having trouble grasping the intuition of it.
  • Why would there be such a notable difference between Models 5 and 9?

Reprex follows, written in R.


data.ex <- tibble(
  x = 1:100,
  w = c(rep(0, 10), 1:20, rep(0, 20), 1:20, rep(0, 30)),
  w2 = w + 1,
  y = dnorm(x, mean = rep(c(20, 60), each = 50), sd = 3) + (seq(0, 1, length = 100)^2) / 2 + rnorm(100, sd = 0.01)

models <- tibble(
  model = 1:9,
  formula = c('y ~ s(x)', 'y ~ s(w)', 'y ~ s(x) + s(w)', 'y ~ s(x, w)', 'y ~ s(x, by = w)', 'y ~ s(w, by = x)', 'y ~ x + s(w)', 'y ~ w + s(x)', 'y ~ s(x, by = w2)'),
  gam = map(formula, function(x) gam(as.formula(x), data = data.ex)),
  data.to.plot = map(gam, function(x) cbind(data.ex, predicted = predict(x)))

plot.models <- unnest(models, data.to.plot) %>%
  mutate(facet = sprintf('%i: %s', model, formula)) %>%
  ggplot(data = ., aes(x = x, y = y)) +
  geom_point() +
  geom_line(aes(y = predicted), color = 'red') +
  facet_wrap(facets = ~facet)
  • 1
    $\begingroup$ It's a little anti-social to folks here to use the tidyverse package as a dependency of the reprex; I use quite a few of those packages and still I required an install fest just to run your code. Minimal, i.e. listing just the packages needed, would have been more useful. That said, thanks for the reprex; I'm just running it now $\endgroup$ Commented Dec 11, 2017 at 19:52

1 Answer 1


Q1 What's the difference between models 3 and 4?

Model 3 is a purely additive model

$$y = \alpha + f_1(x) + f_2(w) + \varepsilon$$

so we have a constant $\alpha$ plus the smooth effect of $x$ plus the smooth effect of $w$.

Model 4 is a smooth interaction of two continuous variables

$$y = \alpha + f_1(x, w) + \varepsilon$$

In a practical sense, Model 3 says that no matter what the effect of $w$, the effect of $x$ on the response is the same; if we fix $x$ at some known value and vary $w$ over some range, the contributions from $f_1(x)$ to the fitted model remains the same. Verify this if you want by predict()-ing from model 3 for a fixed value of x and a few different values of w and use the type = 'terms' argument of the predict() method. You'll see a constant contribution to the fitted/predicted values for s(x).

This is not the case of model 4; this model says that the smooth effect of $x$ varies smoothly with the value of $w$ and vice-versa.

Note that unless $x$ and $w$ are in the same units or we expect the same wiggliness in both variables, you should be using te() to fit the interaction.

m4a <- gam(y ~ te(x, w), data = data.ex, method = 'REML')

pdata <- mutate(data.ex, Fittedm4a = predict(m4a))
ggplot(pdata, aes(x = x, y = y)) +
  geom_point() +
  geom_line(aes(y = Fittedm4a), col = 'red')

enter image description here

In one sense, model 4 is fitting

$$y = \alpha + f_1(x) + f_2(w) + f_3(x, w) + \varepsilon$$

where $f_3$ is the pure smooth interaction of the "main" smooth effects of $x$ and $w$, and which have been removed, for sake of identifiability, from the basis of $f_3$. You can get this model via

m4b <- gam(y ~ ti(x) + ti(w) + ti(x, w), data = data.ex, method = 'REML')

but note this estimates 4 smoothness parameters:

  1. the one associated with the main smooth effect of $x$
  2. the one associated with the main smooth effect of $w$
  3. the one associated with the marginal smooth of $x$ in the interaction tensor product smooth
  4. the one associated with the marginal smooth of $w$ in the interaction tensor product smooth

The te() model contains just two smoothness parameters, one per marginal basis.

A fundamental problem with all these models is that the effect of $w$ is not strictly smooth; there's a discontinuity where the effect of $w$ falls to 0 (or 1 in w2). This is showing up in your plots (and the one I show in detail here).

Q2 What, exactly, is "by" doing here?

by variable smooths can do a number of different things depending on what you pass to the by argument. In your examples the by variable, $w$ is continuous. In that case what you get a varying coefficient model. This is a model in which the linear effect of $w$ varies smoothly with $x$. In equation terms this is what your, say model 5, is doing

$$y = \alpha + f_1(x)w + \varepsilon$$

If this isn't immediately clear (it wasn't to me when I first looked at these models) for some given values of $x$ we evaluate the smooth function at this value and this then become the equivalent of $\beta_1 w$; in other words, it is the linear effect of $w$ at the given value of $x$, and those linear effects vary smoothly with $x$. See section 7.5.3 in the second edition of Simon's book for a concrete example where the linear effect of a covariate varies a smooth function of space (lat and long).

Q3 Why would there be such a notable difference between Models 5 and 9?

The difference between models 5 and 9 I think is simply due to multiplying by 0 or multiplying by 1. With the former, the effect of the only term in the model $f_1(x)w$ is 0 because $f_1(x) \times 0 = 0$. In model 9 you have $f_1(x) \times 1 = f_1(x)$ in those areas where there is no contribution from the gaussians of $w$. As $f_1(x)$ is an ~ exponential function, you get this superimposed upon the overall effect of $w$.

In other words, model 5 contains zero trend everywhere $w$ is 0, but model 9 includes an ~ exponential trend everywhere $w$ is 0 (1), upon which the varying coefficient effect of $w$ is superimposed.

  • $\begingroup$ That's a helpful answer to Q1, thanks! The choice of sums of 1D smooths or a single 2D smooth seems roughly (again, roughly) analogous to main effects vs. interactions in standard linear modeling. But this makes the existence of the by parameter all the more perplexing. $\endgroup$
    – jdobres
    Commented Dec 11, 2017 at 20:19
  • $\begingroup$ I've added something on Q2 now, which hopefully explains what those models are doing. I'll look at Q3 now. $\endgroup$ Commented Dec 11, 2017 at 20:25
  • $\begingroup$ I think answer to Q3 is just a simple issue of arithmetic interacting with what the varying-coefficient part of the models in 5 and 9. $\endgroup$ Commented Dec 11, 2017 at 20:33
  • $\begingroup$ Very helpful! To clarify Q2, are you saying that the predictor provided in the "by" argument essentially becomes an extra coefficient for the result of the smoothed predictor? I suspect that my intuition is wrong, as it should lead to Model 5 looking similar to Model 2. $\endgroup$
    – jdobres
    Commented Dec 11, 2017 at 20:41
  • $\begingroup$ In M2 you have $\alpha + f_1(w)$ and in M5 you have $\alpha + f_1(x)w$, the former is a smooth of $w$, the latter says the linear effect of $w$ varies as a smooth function of $x$. These are quite different; notice that in 5 the effect of $w$ is linear (there's no decrease when $w$ > 10), but that linearity is modified smoothly via $x$. In M2 you have a smooth effect of $w$. $\endgroup$ Commented Dec 11, 2017 at 20:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.