Autocorrelation in GAM (R)

Question

I am fitting a GAM to a time series, to check for any trend and seasonality. I expect my data to be autocorrelated and as I am fitting some models I need to check if there is autocorrelation in the residuals, as part of model diagnostics.

I've read that it's recommended to apply the ACF() function to normalized residuals to accurately assess autocorrelation. Is this correct? If so, how can I get normalized residuals from a gam()-fitted (from MGCV package) model in R?

My model is something like this:

    gam(abundance ~ species + s(year_factor, bs= "re") + 
    s(julian_day_number, by=species)+s(month_as_numeric, bs="cc", by=species), 
    data=data, method="REML", family=tw())

I have 5 species and I want to see if there are significant trends and seasonality in each.

• Trend is Julian day number
• Seasonality is month as numeric
• Random effect is year as factor

I have tried to plot ACF on Pearson residuals,

residuals(model, type="Pearson")

as they seemed the closest type to normalized residuals after reading some help documents on R Documentation.

Answer 1

You only need normalized or standardized residuals when you have fitted a correlation function in the model. As you haven't done that here, you can just use the standard residuals for the ACF diagnostic. The Pearson residuals would work too; they are standardized by the estimated variance of the observations.

Answer 2

I agree with Gavin here that the type of residual you extract is probaly less crucial here. But given that you're fitting the model to several different time series (one for each level of species), you'll want to be sure that you inspect possible autocorrelation for each series separately. If you do find unmodelled temporal autocorrelation left, you can impose an AR(1) model for the working residuals in bam(), where you can also supply a vector of class logical to the AR.start argument (i.e. arrange the data by species and then by julian_day_number and add a new variable that sets to TRUE if that row represents the first observation in a time series and FALSE otherwise). The downside here of course is that you have to specify the AR1 coefficient, it cannot be estimated. I suspect there are ways to do this with gamm() as well but I'll have to defer to others to explain the correct syntax.

If you prefer more complex dynamics or would like to make future predictions however, I must shamefully recommend my package {mvgam}. It can handle all of the nonlinear effects you might be interested in but can also incorporate latent temporal dynamics such as autoregressive, Gaussian Process or multivariate temporal processes. It is also designed specifically for working with multiple time series.