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:
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
ywhere a gaussian is present.
whas values of 1-20 where
xis between 11 and 30, and 51 to 70. Otherwise,
w + 1, so that there are no 0 values.
mgcv package makes it easy to specify a number of possible models for these data:
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
w as 1D smooths, producing a nice fit. Model 4 uses
w in a 2D smooth, also giving a nice fit. These two models are very similar, though not identical.
Model 5 models
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.
library(magrittr) library(tidyverse) library(mgcv) set.seed(1222) 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) print(plot.models)