0
$\begingroup$

I'm trying to identify seasonal and diel patterns in multi-year data from different locations using a mgcv::gam. I have continuous data at one site and recording gaps during the summer months for the others.

I fit the gam like this:

mod1 <- gam(y ~ s(julian, by=site) + s(Hour, 
    bs='cc', by=site) + site + year, 
    data=dat, method='REML', select=TRUE)

GAM1

In the effect plots for julian day (seasonal pattern), I obviously have times with no data for some of the sites, but do see an 'assumed effect'. I'm not completely sure why the rug claims there was data. Is there any way to suppress this to not give false visual interpretations? Especially since the assumed effect is quite large with -40/-60, while in reality it probably is closer to what is seen at site 3.

Also, can I actually take the 'site' factor out of the smoother for julian despite having the gaps for some of the sites? I'm not expecting underlying drastically different seasonal patterns across sites, but lack the summer data, so wasn't sure whether that would be statistically sound.

mod <- gam(y ~ s(julian, bs='cc') + s(Hour, 
                 bs= 'cc', by=site) + site + 
               year, data=dat, method='REML', 
               select=TRUE)

GAM2

The data can be found here: https://www.dropbox.com/s/c3ypj3zj2na7pjo/dat.Rda?dl=0

$\endgroup$

1 Answer 1

1
$\begingroup$

Q1

In the effect plots for julian day (seasonal pattern), I obviously have times with no data for some of the sites, but do see an 'assumed effect'. I'm not completely sure why the rug claims there was data. Is there any way to suppress this to not give false visual interpretations?

The rug is showing the support of the data for the basis used to create the smooths. Even though in one sense the model includes three separate smooths, each with their own smoothness parameter(s), $\boldsymbol{\lambda}_j$, the basis functions themselves are shared across each of the smooths.

The only option to suppress this is to set rug = FALSE in the plot.gam() call. If you want to use other packages to visualise GAM fits, the {gratia} package has a draw() method for gam() fits which does separate out the observations into the panel for the appropriate factor level/smooth combination.

Something looks very off with this model; you haven't set the values in those missing months to some value have you? I can't recall that I've ever seen fits that are so poor for smooths of season.

Why isn't s(julian) a cyclic smooth here but it is in the second model with just a single seasonal smoother below (although the figure you show suggests that it didn't arise from the claimed model as the smooth is nowhere near joining up smoothly at the end points and your model code has errors/typos, suggesting you wrote this out rather than copy/pasting from your actual script.)

Q2

Also, can I actually take the 'site' factor out of the smoother for julian despite having the gaps for some of the sites? I'm not expecting underlaying drastically different seasonal patterns across sites, but lack the summer data, so wasn't sure whether that would be statistically sound.

mod <- gam(y ~ s(julian, bs = 'cc') + s(Hour, bs = 'cc', by = site) +
             site + year,
           data = dat, method = 'REML', select = TRUE)

(edited code for style and to fit over multiple lines)

You can; the statistical Gods won't smite you for fitting this model. They might get a bit tetchy if you then did any inference about the estimated model terms if you didn't check the independence assumption of this model. You should plot the residuals from mod against julian separately for each site; if the assumptions of the model hold, there should be no effect of julian in the residuals across the three sites. The residuals should be evenly and equally scattered about the 0 line throughout the range of julian.

An alternative might be to fit a model with a global or common seasonal smooth plus site-specific smooths. There are multiple ways to do this, each with its own interpretation, but two that jump to mind are:

  1. use a first-order derivative penalty on the site-specific smooths so that the factor by smooths are interpreted as differences from the common seasonal smooth. This model would be fitted as

     mod <- gam(y ~ site + year +
                  s(julian, bs = "cc") +
                  s(julian, bs = "cc", by = site, m = 1) +
                  s(Hour, by = site, bs = "cc"),
                data = dat, method = "REML", select = TRUE)
    

    The m = 1 bit is what tells {mgcv} to use a first derivative penalty for the smooth, so that it is interpreted as a difference from a flat function, which means difference from the common smooth.

  2. a second option would be to set the site with full data as the reference site (so relevel() the site factor so that the reference is "site3"), then create a new ordered factor for site, say site_o. Then fit the model as

     mod <- gam(y ~ site_o + year +
                  s(julian, bs = "cc") +
                  s(julian, bs = "cc", by = site_o) +
                  s(Hour, by = site, bs = "cc"),
                data = dat, method = "REML", select = TRUE)
    

    The ordered factor instructs gam() to create a more ANOVA-like parameterisation of the model, where the first smooth is for the reference level (hence the relevel()) and then the by smooths are for smooth difference between each of the other levels and the reference level, like the default contrasts in a lm() produces a parameterisation wherein the intercept term is for the reference group and the other terms for the factor are differences between the reference level and the stated level of the factor

Also, fyi it is select = TRUE not Select = T note the argument name is lowercase: none of your models in the post are doing any feature selection because of this. Again, I presume this is just a typo here, but worth double checking in your actual code as this sort of error will pass by quietly because of the S3 class system and the use of ... which gobbles up arguments that don't match any of the formal arguments to gam()

$\endgroup$
5
  • $\begingroup$ Thank you! I hope my various responses don't get lost because I can't format a comment nicely and stackexchange wants me to reply in a comment rather than an answer. Q1) I know I can supress the rug with rug=F but that still plots the predicted smooth in the effects plot, which is just wrong and major guess work from the model due to the lack of data. I worry about the suggestive interpretation in a publication for example. $\endgroup$
    – Anke
    Mar 17, 2022 at 19:24
  • $\begingroup$ Q1) I have not set any values for the missing month and I don't know why the ends don't fit nicely. This is exactly from my code, I only changed variable names for some anonymity. I didn't make julian cyclic in the first model because I thought it wasn't cyclical data since I had the recording gaps. I did run it with bs='cc' by site and it looked pretty much the same, though. $\endgroup$
    – Anke
    Mar 17, 2022 at 19:26
  • $\begingroup$ Q2) Thank you, I'll test those. And thanks for catching the typo in select, this too was actually directly from my code and of course is a major oversight! $\endgroup$
    – Anke
    Mar 17, 2022 at 19:26
  • $\begingroup$ You second model has an obvious typo (missing quote after the "bs = 'cc` in the julian day smooth), which makes me wonder how that code produced anything... Re Q1, it wasn't clear that you meant to suppress plotting the fitted smooth in the gap rather than the rug. I think hiding this terrible fit would be the wrong way to go; there's clearly something very wrong with this model and hiding it isn't going to change that. What is wrong is hard to say without the data and code (but I understand you can't share that). You shouldn't be publishing this model! $\endgroup$ Mar 17, 2022 at 20:56
  • $\begingroup$ The missing quote was indeed a typo that didn't exist in my original code. I uploaded the data to dropbox and added the link in my original post. Of course this isn't going to be published like this, as you rightly pointed out there seems to be quite a few things wrong with the model(s), still. $\endgroup$
    – Anke
    Mar 17, 2022 at 21:22

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.