r time-series gam cyclic stl-decomposition

can I estimate a time varying seasonal effect in R with GAMM?

I would like to use a generalized additive model to investigate time-series data in R. My data are monthly and I would like to estimate a seasonal effect and a longer run trend effect. I have followed some helpful posts by Gavin Simpson here and here:

My data look like this:

I have the full data set available on my github page:

I have attempted to specify a generalized additive model with smooth seasonal and trend terms as follows:

    df <- read.csv('trips.csv')
    head(df)
    # A tibble: 276 × 2
     date ntrips
   <date>  <int>
    1  1994-01-01    157
    2  1994-02-01    169
    3  1994-03-01    195
    4  1994-04-01    124
    5  1994-05-01    169

    #add a time column
    trips <- tbl_df(trips) %>% mutate(time=as.numeric(date))

    mod1 <- gamm(ntrips~s(month,bs="cc",k=12) + s(time),data=trips)

I extracted the estimate of the seasonal effect as follows:

    pred <- predict(mod1$gam,newdata=trips,type="terms")
    seas <- data.frame(s=pred[,1],date=trips$date)
    ggplot(seas,aes(x=date,y=s)) + geom_line()

This plot is included below:

My question is: in the original data the seasonal peaks move around a little from year to year. In the embarassingly simple GAM I have specified the seasonal effect is constant. Is there a way to accommodate time varying seasonality with a GAM?

I have analyzed these data using the STL approach of Cleveland et al.:

Using the STL paradigm, how wiggly or smooth one allows the seasonal effects to be seems to be a matter of preference or choice. I would prefer if I could allow the data to tell me the difference between random error and a shifting seasonal peak. GAMS seem better suited to this goal as they lend themselves more readily to statistical model fitting-type exercises...but I would like to know if there is a parameter in the R package for fitting gams that allows time varying seasonal effects.

Solution

The answer is: yes, a GAM model can be formulated for the question you are interested in. If we assume that the trend and seasonal components of the model interact smoothly we have a the smooth equivalent of a continuous-continuous interaction. Such interaction can be fitted in a GAM using a tensor product of the two marginal smooths:

the seasonal cyclic smooth, and
the longer-term trend smooth

Incidentally, I have further blog posts on these:

Read those for more detail, but the essential aspects are to fit the following model:

## fix the knots are just passed the ends of the month numbers
## this allows Dec and Jan to have their own estimates
knots <- list(month = c(0.5, 12.5))

## original model, fixed seasonal component
m1 <- gam(ntrips ~ s(month, bs="cc", k=12) + s(time), data = trips,
          knots = knots)

## modified model with
m2a <- gam(ntrips ~ te(month, time, bs = c("cc","tp"), k = c(12, 10)), data = trips,
          knots = knots))

An alternative to the second model is an ANOVA-like decomposition of the two main effects plus interaction. In the modified model above, all three components are bound up in the single tensor product smooth, the te() part of the model.

The ANOVA-like decomposition variant would be fitted using

m2b <- gam(ntrips ~ ti(month, bs = 'cc', k = 12) +
             ti(time, bs = 'tp', k = 10) +
             ti(month, time, bs = c("cc","tp")), data = trips,
           knots = knots)

The third ti() then is the smooth interaction, separated from the main smooth effects of the seasonal and long-term trend.

I've shown these fitted using gam() but they can be used with gamm() also if you need to include an ARMA process for the model residuals.