I am currently modellingusing the 'mgcv' package in R. My response variable is called log.tr, representing the log of residence time. My data looks a little bit like this:

logtr <- seq(0.67, 4.29, length.out = 100)
day_ <- sample(1:365, 100, replace = TRUE)  # Random day values
year_ <- sample(2000:2020, 100, replace = TRUE)  # Random year values
TEMPERATURE <- rnorm(100, mean = 25, sd = 5)  # Random temperature values
gam_tres_df <- data.frame(log.tr = logtr, day_ = day_, year_ = year_, TEMPERATURE = TEMPERATURE)

I am attempting to fit a generalized additive model (GAM) using the 'gam' function from 'mgcv'.

The model formula I am using is:

gam(log.tr ~ s(day_, k = 40, bs = 'cc') + s(year_, k = 12) + s(TEMPERATURE, k = 40) + s(day_, year_), method = 'REML', data = gam_tres_df)

In this formula, I have included smooth terms for the variables day_, year_, TEMPERATURE, and an interaction term between day_ and year_. I have chosen specific degrees of freedom (k) for each smooth term.

However, I suspect that the variables day_ and year_ may be dependant. I am considering adding a correlation structure to the model to account for this potential correlation. Specifically, I am thinking of using the corARMA function with the form ~ 1|Year and an autoregressive moving average (ARMA) model with a lag of 3 (p = 3).

My question is whether I should include this correlation structure (corARMA(form = ~ 1|Year, p = 3) or maybe (form = ~ day | year)) in my GAM model or if the current model specification without the correlation structure is appropriate. Also I wanted to know how to define this correlation structure for seasonal daily data; since i want to consider residual variation from year to year.

Please let me know if you need further information or have any additional questions!


1 Answer 1

  1. You're likely better off modelling the residence time itself, not it's log transform, using a non-Gaussian family, but that could well be tricky if you do actually need the autocorrelation structure.

  2. You can't use the correlation with gam(); you can only use this with gamm(), FYI; gam() will silently ignore this argument if you provide it.

  3. There shouldn't really be a relationship between day_ and year_, their effects are unlikely to be similar; one is a within-year effect, the other a between-year effect. There may be an interaction between these two variables, if say the seasonal pattern has changed over time due to say climate change. These variables are not dependent however.

  4. The correlation structure available via gamm() (provided by the nlme package) don't account for correlation or dependence among covariates. They account for un-modelled correlation among observations.

  5. The correlation structures available do not cover seasonal ARMA models, only ARMA models.

  6. corARMA(form = ~ 1|year_, p = 3) and corARMA(form = ~ day_ | year_), p = 3) should give the same results if the data are in time order within year_. Both will describe an AR(3) nested within year, using the same autocorrelation parameter $\rho$

  7. As mentioned in 4. you only need it if there is unmodelled temporal autocorrelation. With 40 df for the day_ term, you can model quite a complex seasonal effect so maybe you won't need the AR(3)? You might want to consider using te(day_, year_, ...) to allow the seasonal signal to vary with the trend, and then you might not need the autocorrelation.

  8. Whether you need the autocorrelation structure will depend on how wiggly you want the smooths of day_ and year_ to be; all else equal, you can model the autocorrelation with wiggly smooth functions and then you won't need the autocorrelation structure. But if you are trying to estimate the seasonal and between-year trends without this autocorrelation (so you want simpler, smoother functions) then it is likely that you'll need the autocorrelation structure if you have data recorded at sufficient time resolution to separate the larger scales of temporal variation (seasonal and long-term trend) from the short-scale variation of the autocorrelation.

  9. You can fit models with no AR(p), and AR(1), an AR(2), and then an AR(3) and use tools to select among these 4 models. I have some example code here showing how to do this (for models with AR(1) vs AR(2) for example): https://fromthebottomoftheheap.net/2016/03/25/additive-modeling-global-temperature-series-revisited/

  • $\begingroup$ Thank u Gavin u the goat!! Really appreciate the effort you put into helping others. The only point that is not clear to me is 8. I mean, now I know that I don't need to introduce the autocorrelation, but I would like to study how the seasonal variation changes over time (I have 25 years of daily data). At first I was using Loess, but I read your post and I decided to dive into this GAMs stuff. But when i introduce the tensor product te(day_, year_) it seems like a lot of the trend goes to this tensor product (seasonal component variation). Is there a way to adjust this? $\endgroup$
    – Pablo
    May 26, 2023 at 10:05
  • $\begingroup$ Point 8. arises because there is an equivalency between wiggly smooth functions and autocorrelation. The latter is the degree to which you have short-run behaviour where samples are higher (lower) than expected under an assumption of i.i.d noise, in other words the observations close in time around any observed value are likely to be more similar to one another than expected under an independence model for the data. Neighbours in time are more similar than expected. That similarity among neighboring points is what a smooth models. Hence wiggly smooths and autocorrelation can be the same thing $\endgroup$ May 26, 2023 at 10:20
  • $\begingroup$ So in some cases you can't have both in the model (a complex trend and an autocorrelation process) as they both model the same thing. You can only do what you want if the scale at which the trend works is different to the scale of the autocorrelation. For example, in the global climate example from the blog, the trend is at a larger temporal scale (it varies over many years not over a few samples) while the autocorrelation is over over a few samples (smaller temporal scales). If they operate at similiar scales you likely can't include both as they are the same thing. $\endgroup$ May 26, 2023 at 10:22
  • $\begingroup$ I don't understand what you mean by the te() term including all the trend - it's supposed to. You have to predict for the same day of year over all years to get the trend for that day of year because the trend is different for every day of year. If you want a model which still has the average change over years then you can do y ~s(day_) + s(year_) + ti(day_, year_) instead of y ~ te(day_, year_) (I have excluded all the bs, k bits in those formulas), but the te() model is simpler as it involves fewer smoothing parameters. $\endgroup$ May 26, 2023 at 10:25
  • $\begingroup$ Sorry for not explaining myself well. What I was referring to is that after applying the stl() function the trend was positive (which is what is happening because it is a political decision) while the GAM results say that this trend is negative and attribute these changes in the last two years to the seasonal variation component 'te(day_, year_)' $\endgroup$
    – Pablo
    May 26, 2023 at 10:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.