Presenting 2D smoother of GAM

Question

I had a look around and couldn't find the answer to my problem, so hopefully this is a new question.

I tried to fit a GAM with two continuous explanatory variables, one of them is Day of the year and therefore uses a cyclic spline smoother. The data is normal distributed.
It would be great to present the result as a 3D surface with x and y as explanatory and z as response variables. I fitted two models using R library mgcv:

Model1<-gam(y ~ s(a, bs = "cc") + s(b, bs = "cs"))

Model2<-gam(y ~ te(a, b, bs = c("cc", "cs")))

Model2 contains a tensor, as the normal s() function can't account for the cyclicity. As this is a proper interaction, I am sure you can present the model as a 3D surface, e.g.

data<-expand.grid(a,b)
prediction<-predict(Model2, data)
ggplot(data, aes(x=a, y=b))+geom_raster(aes(fill = prediction))

The problem is, in my model output I don't get single terms for a and b but only for their tensor product. But would it be incorrect to present Model1 in the same way as Model2? And would it be incorrect or rather overfitting the model to use something like this:

Model3<-gam(y ~ te(a, b, bs = c("cc", "cs")) + s(a, "cc") + s(b, "cs"))

I know this might be hard to answer without the actual data. I did model selection and model comparisons based on AIC, Model 3 has the lowest AIC. But I am sometimes not sure about basic rules of model fitting, what is allowed and what isn't.

Maybe there is no single answer, but some tips and recommendations would be appreciated. Many thanks :)

Answer 1

Your Model 3 is asking for numerical trouble; basis functions used to represent the "main effects" of $f(a)$ and $f(b)$ are also found in the tensor product smooth (not tensor), which can lead to identifiability problems. {mgcv} does it's best to avoid these problems and detect repeated basis functions and remove them from the model, but this process isn't infallible and even if it works well, you can get inflated credible intervals on the estimated smooths.

A better approach to Model 3 is to us the ti() smooth, for a tensor product interaction

gam(y ~ s(a, "cc") + s(b, "cs") + ti(a, b, bs = c("cc", "cs")))

The basis for the ti() smooth has had these "main effects basis functions" removed from the basis so there is no identifiability issue. The only problem with this model is that it uses more smoothing parameters than the equivalent

gam(y ~ te(a, b, bs = c("cc", "cs")))

and hence the two models often show some subtle differences because the former (ti()-containing) is actually a more complex one (te()-containing).

I'm not sure what you mean by "in my model output I don't get single terms for a and b but only for their tensor product", but if you mean something like the situation where you have multiple terms that include a or b:

gam(y ~ s(a) + s(b) + te(a, c) + te(b, d) + te(a, b)) # wrong

then from the above discussion about ti() we should ideally proceed by fitting models with all the main effects specified via s() and then use ti() terms to fit the interactions

m <- gam(y ~ s(a) + s(b) + s(c) + s(d) +
             ti(a, c) + ti(b, d) + ti(a, b))

Now say you want to plot a 2d surface of the main effects plus interaction for a and b, you need to predict() from the model but exclude the effects of the smooth terms that aren't involved in the specific interaction**

# note, if you have `NA` in any of the covariates you'll need to handle that
new_df <- with(data,
               expand.grid(a = seq(min(a), max(a), length = 50),
                           b = seq(min(b), max(b), length = 50),
                           c = mean(c),
                           d = mean(d)))
p <- predict(m, newdata = new_df,
             exclude = c("s(c)", "s(d)", "ti(a,c)", "ti(b,d)"),
             se.fit = TRUE)
p <- cbind(new_df, as.data.frame(p))

Note that smooth must be named exactly as {mgcv} names them internally and shows them the output from summary(m).

These predictions will include the intercept (constant), but you can use type = "terms" and modify the code above to get predicted values without the intercept/constant by summing up the rows of matrix of predicted values it returns (standard errors are a little more tricky to get this way but if you are plotting a 2d surface you aren't showing an interval for that surface at the same time anyway - but it can be done).

The interpretation here is that we have estimated the marginal effect of $f(a,b)$ correctly given the effects of the other covariates on the response, but then we ignore the additive effects of those other terms for the purposes of visualizing the partial effect of $f(a,b)$. This is the same principle as we use to get partial effects plots of smooths produced by plot(m) say, it's just that we've got the intercept included as well.

**You could also condition the predictions on some values of the other covariates, which would require you to specify some values for c and d that you want to condition on, rather than the arbitrary values that were plugged in above. To do this, using the values I plugged in above, would result in outputs that are conditional upon c and d being at their men values in the data used to fit the model