Interpretation of contour plots (mgcv)

When we plot a GAM model using the mgcv package with isotropic smoothers, we have a contour plot that looks something like this:

x axis for one predictor,
y axis for another predictor,
the main is a function s(x1, x2) (isotropic smother).

Suppose that in this model we have many other isotropic smoothers like:

y ~ s(x1, x2) + s(x3, x4) + s(x5, x6)

My doubts are: when interpreting the contour plot for s(x1, x2), what happens to the others isotropic smoothers? Are they "fixed at their medians"? Can we interpret a s(x1, x2) plot separately?

Solution

Because this model is additive in the functions you can interpret the functions (the separate s() terms) separately, but not necessarily as separate effects of covariates on the response. In your case there is no overlap between the covariates in each of the bivariate smooths, so you can also interpret them as the effects of the covariates on the response separately from the other smoothers.

All of the smooth functions are typically subject to a sum to zero constraint to allow the model constant term (the intercept) to be an identifiable parameter. As such, the 0 line in each plot is the value of the model constant term (on the scale of the link function or linear predictor).

The plots shown in the output from plot.gam(model) are partial effects plots or partial plots. You can essentially ignore the other terms if you are interested in understanding the effect of that term on the response as a function of the covariates for the term.

If you have other terms in the model that might include one or more covariates in another terms, and you want to look at how the response changes as you vary that term or coavriate, then you should predict from the model over the range of the variables you are interested in, whilst holding the other variables at some representation values, say their means or medians.

For example if you had

model <- gam(y ~ s(x, z) + s(x, v), data = foo, method = 'REML')

and you want to know how the response varied as a function of x only, you would fix z and v at representative values and then predict over a range of values for x:

newdf <- with(foo, expand.grid(x = seq(min(x), max(x), length = 100),
                               z = median(z)
                               v = median(v)))
newdf <- cbind(newdf, fit = predict(model, newdata = newdf, type = 'response'))
plot(fit ~ x, data = newdf, type = 'l')

Also, see ?vis.gam in the mgcv package as a means of preparing plots like this but where it does the hard work.