How to include an interaction term in GAM?

Question

The following code evaluates the similarity between two time series:

set.seed(10)
RandData <- rnorm(8760*2)
America <- rep(c('NewYork','Miami'),each=8760)

Date = seq(from=as.POSIXct("1991-01-01 00:00"), 
           to=as.POSIXct("1991-12-31 23:00"), length=8760)

DatNew <- data.frame(Loc = America,
                     Doy = as.numeric(format(Date,format = "%j")),
                     Tod = as.numeric(format(Date,format = "%H")),
                     Temp = RandData,
                     DecTime = rep(seq(1, length(RandData)/2) / (length(RandData)/2),
                                   2))
require(mgcv)
mod1 <- gam(Temp ~ Loc + s(Doy) + s(Doy,by = Loc) +
  s(Tod) + s(Tod,by = Loc),data = DatNew, method = "ML")

Here, gam is used to evaluate how the temperature at New York and Miami vary from the mean temperature (of both locations) at different times of the day. The problem that I now have is that I need to include an interaction term which shows how the temperature of each location varies throughout at the day for different days of the year. I eventually hope to display all of this information on one graph (for each location). So, for Miami I hope to have one graph that shows how the temperature varies from the mean during different times of the day and different times of the year (3d plot?)

Answer 1

The "a" in "gam" stands for "additive" which means no interactions, so if you fit interactions you are really not fitting a gam model any more.

That said, there are ways to get some interaction like terms within the additive terms in a gam, you are already using one of those by using the by argument to s. You could try extending this to having the argument by be a matrix with a function (sin, cos) of doy or tod. You could also just fit smoothing splines in a regular linear model that allows interactions (this does not give the backfitting that gam does, but could still be useful).

You might also look at projection pursuit regression as another fitting tool. Loess or more parametric models (with sin and/or cos) might also be useful.

Part of decision on what tool(s) to use is what question you are trying to answer. Are you just trying to find a model to predict future dates and times? are you trying to test to see if particular predictors are significant in the model? are you trying to understand the shape of the relationship between a predictor and the outcome? Something else?

Answer 2

For two continuous variables then you can do what you want (whether this is an interaction or not I'll leave others to discuss as per comments to @Greg's Answer) using:

mod1 <- gam(Temp ~ Loc + s(Doy, bs = "cc", k = 5) + 
                         s(Doy, bs = "cc", by = Loc, k = 5) + 
                         s(Tod, bs = "cc", k = 5) + 
                         s(Tod, bs = "cc", by = Loc, k = 5) +
                         te(ToD, Day, by = Loc, bs = rep("cc",2)),
            data = DatNew, method = "ML")

The simpler model should then be nested within the more complex model above. That simpler model is:

mod0 <- gam(Temp ~ Loc + s(Doy, bs = "cc", k = 5) + 
                         s(Doy, bs = "cc", by = Loc, k = 5) + 
                         s(Tod, bs = "cc", k = 5) + 
                         s(Tod, bs = "cc", by = Loc, k = 5),
            data = DatNew, method = "ML")

Note two things here:

1. The basis-type for each smoother is stated. In this case we would expect that there are no discontinuities in Temp between 23:59 and 00:00 for ToD nor between DoY == 1 and DoY == 365.25. Hence cyclic cubic splines are appropriate, indicated here via bs = "cc".
2. The basis dimension is stated explicitly (k = 5). This matches the default basis dimension for each smooth in a te() term.

Together these features ensure that the simpler model really is nested within the more complex model.

For more see ?gam.models in mgcv.