# How to include an interaction term in GAM?

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?)

-
You might find the answer to this question stats.stackexchange.com/questions/18937/… relevant. –  jbowman Jul 21 '12 at 14:42

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?

-
Suppose you have two predictors $x_1,x_2$ - isn't $$y = f_1(x_1) + f_2(x_2) + f_3(x_1x_2) + \varepsilon$$ still an gam? This could be thought of as an 'interaction' in some sense. Also, I think the gam package allows you to fit models like $$y = f_1(x_1) + f_2(x_2) + f_3(x_1)x_2 + f_4(x_2)x_1 + \varepsilon$$ (which I think is what the by argument is doing). Is that still a gam? –  Macro Jul 21 '12 at 17:43
@Macro, It depends on if you want to split hairs or not (well technically what you wrote would maybe be am's since you are not really using the g part). –  Greg Snow Jul 21 '12 at 18:02
@Macro, I'm sure more has been said on the subject, but see page 4 on GAMs of this Venable's article, Exegeses on Linear Models. It isn't clear exactly how non-additive main effects and interaction effects are simultaneously identified. –  Andy W Jul 22 '12 at 2:20
many thanks for your comments. My main aim here is to not predict future values but to simply see how each time series alters from the mean of both series. For example, the outcome of mod1 shows how the time series at each location varies from the mean firstly at the time scale of day of year (Doy) and then time of day (Tod). From this I would like to see how each series varies at a function of Doy and Tod which I expect the time series to differ greatly during the summer period. –  Kate Jul 24 '12 at 9:40

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.

-