I am modelling the occurrence of a species at 5 different sites on an hourly basis (presence/absence), based on a range of temporal predictors (e.g. time of the year, day/night cycle, tides ...). Covariates are indicated by x1, x2, ... in the code below. For information, I have ~70 000 data points.

I am using a HGAM structure, as introduced in Pedersen et al.2019. For each covariate, I first investigated different specification options (global smoother or not, shared penalty or not ...), and selected the best one based on AIC and my research question.

When putting all of the terms together in the model, I end up with a structure like this:

model <- bam(response ~ offset(log(offset)) + s(year, bs="re") + 
      Site + s(x1, m=2, bs="cc", k=8) + 
      s(x1, Site, bs = "fs", xt=list(bs="cc"), m=2) + 
      s(x2, bs = "cc", by=Site, m=2, k=8) + 
      s(x3, m=2, bs="cc", k=10) + 
      s(x3, by = Site, bs = "cc", m=1, k=10) + 
      s(x4, bs = "cc", by=Site, m=2, k=8) + 
      s(x5, bs = "tp", by=Site, m=2, k=8), family = "binomial", 
      data = data.all, method="REML", cluster=cl, select=TRUE)

The explained deviance of the model is only about 13%. Given the high temporal resolution (hourly), I expect some temporal autocorrelation in the residuals. I read in the bam() documentation, and in other posts (here and there), that this could be specified using the rho argument in bam() following the order of the dataset.

  • However, I am not sure this applies when the model as a hierarchical structure ? I know it is feasible in gamm(), but the problem then is that I cannot specify multiple factor-smoother interaction terms as currently written in the model ...
  • Is it something that could be tackled in brms ? Or any other way ?

For information, here is the output of the acf and pcaf functions plot:

enter image description here enter image description here

  • $\begingroup$ Is that the ACF of your data or the normalised/standardised residuals? $\endgroup$ Apr 6, 2023 at 7:50
  • $\begingroup$ This is the ACF for the model residuals plotted with the following: acf(resid(model), lag.max = 100, main = "ACF") $\endgroup$
    – Timelate
    Apr 6, 2023 at 7:57
  • $\begingroup$ I edited the question to also include the output of pacf(resid(model), lag.max = 100, main = "pACF") $\endgroup$
    – Timelate
    Apr 6, 2023 at 8:03

1 Answer 1

  1. The hierarchical structure is mainly for the "fixed" effects, but you do need to consider it when it comes to an autocorrelation structure. With both correlation in gamm() and rho in bam(), unless told otherwise, the model will assume a single long time series.

    This behaviour may be what you want, but you may, for example in your case, want to nest the correlation structure within day rather than just the temporal ordering of the samples. Such a nesting would say that the autocorrelation structure operates at the lowest grouping level (within day) — you can't opt to nest it within year when you have sub-daily level data unless you create a variable that orders the samples at the sub-daily level within year.

    With correlation you would use corAR1(form = ~ t | year) to indicate you want an AR(1) nested within year, where t orders the observations within year. To get this nested within site you just need to augment the right-hand side of the formula. One way to do this is to create a variable year_site = interaction(year, site, drop = TRUE) in your data and then modify the formula to be corAR1(form = ~ t | year_site). There may be more direct ways to do this within the formula but all the examples I have seen use a single grouping factor. With rho you need to create a logical vector that is FALSE everywhere, except at the first observation of each "series" that you want the correlation structure to be nested within.

    Importantly, you are estimating the same single parameter $\rho$ regardless of how you specify the AR(1), that operates within each level of the grouping variable (or equivalently for other ARMA terms). In that sense, the parameters are global in the way they describe a common autocorrelation structure within each level of the hierarchy.

  2. brms works essentially the same way. See https://paul-buerkner.github.io/brms/reference/autocor-terms.html and then the individual correlation functions linked from that page which document arguments such as time (for ordering the observations) and gr for the grouping factor.

  • $\begingroup$ Hi @Gavin Simpson, thank you very much for the comprehensive answer. I did try after asking my question to specify a rho term in bam(), just as you suggested. The ACF plot appeared the same but I think I may have wrongly defined my 'events' when assigning the TRUE/FALSE values. I went by default for a new event for each new independent time series (so each new deployment at a monitoring site), but that might be too coarse as deployments lasted for months. I could maybe indeed choose daily blocks. Looking at the ACF plots there seems to be some kind of 24h component in the correlation... $\endgroup$
    – Timelate
    Apr 6, 2023 at 8:01
  • $\begingroup$ I also specified year as a random effect mostly because of one predictor (Julian day) for which I have at best two replicates (two years of data), but most others predictors are cyclic over much shorter time periods (tides, day/night) so I don't expect year to have a huge effect on the estimates for these... $\endgroup$
    – Timelate
    Apr 6, 2023 at 8:23
  • 1
    $\begingroup$ For bam() you need to plot the ACF/PACF of the standardised residuals which are in m$rsd (IIRC) for model object m. If you used the code as cited in your replies to my comment on your question, you are getting residuals that don't include the autocorrelation structure/process you specified. $\endgroup$ Apr 6, 2023 at 8:44
  • 1
    $\begingroup$ The main value-add of fitting a single HGAM when there are no global smooths is that you have only a single model to check and validate. Another advantage is that you have a way to statistically compare curves fitted to different groups (for example via gratia::difference_smooths()) which would not be possible if you had separate models for each group $\endgroup$ Apr 6, 2023 at 10:49
  • 1
    $\begingroup$ That's great ! I'll keep this specification then ! Thank you :) $\endgroup$
    – Timelate
    Apr 6, 2023 at 11:10

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.