How to account for temporal autocorrelation in a hierarchical generalized additive model (HGAM/GAMM) with a negative binomial distribution?

Question

I am using a harvest data set (count data) that I am trying to infer population abundance from. This data set contains annual harvest estimates from 20 Zones of varying areas, spanning 34 years. While catch-per-unit-effort is often used as a response variable in similar analyses, I would like to model Effort as a predictor variable and keep my harvest data (Catch) as counts. I would also like to include annual Price as a predictor variable in the model.

Based on diagnostic plots, it is clear that my data fit a negative binomial distribution better than a poisson. There are not many zeros in my data either (~2%).

As of right now, I have separated my data into training (first 27 yrs) and testing (last 7 yrs) data sets. To account for the different areas within each zone, I divided Effort by the area (i.e., Effort_km2) of each Zone. I then scaled and centered the Effort_km2 and Price variables (mean=0, sd=1), and only centered the year data around 0. This is my general model structure using the gam function in the mgcv package:

log_area <- log(area)

m1 <- gam(Catch ~ s(Effort_km2, bs="tp") + 
                  s(Price, bs="tp") + 
                  s(Year, bs="tp") +
                  s(Year, by=Zone, m=1, bs="tp") +
                  s(Zone, bs="re", k=20) + 
                  offset(log_area),   #control for different sized Zone areas
                data = dat,
                family = nb(link="log"),
                method = "REML")

Here, I fit a global smoother for year, as well as a group smoother for each zone that will allow for the wiggliness to vary by zone.

When I examine the normalized residuals from this model using acf function, I see evidence of temporal autocorrelation unless I increase the basis size (k values) of the Price and the two Year covariates up to their near-maximum values (e.g., s(Year, bs="tp", k=25)). (Similarly, the k.check function indicated the k-values need to be that high.) I saw mentioned in a tutorial (by Gavin Simpson) that this would indicate that the model is basically modeling the autocorrelation in the data, which should be accounted for in a different way. I should also note that there are a few NAs in my data.

So my question is how can I develop a GAMM with a negative binomial distribution that can account for temporal autocorrelation? I tried the following model using the gamm function:

th1 <- m1$family$getTheta # estimate theta

m2 <- gamm(Catch ~ s(Effort_km2, bs="tp") + 
                  s(Price, bs="tp") + 
                  s(Year, bs="tp") +
                  s(Year, by=Zone, m=1, bs="tp") +
                  offset(log_area),   #control for different sized Zone areas
                data = dat,
                correlation = corCAR1(form = ~ Year | Zone)
                random=list(Zone=~1),
                family = negbin(theta=th1),
                method = "REML",
                niterPQL = 1000)

However, I get an error message about singularity, which I think comes from having the two Year covariates. Is there a way to specify the year terms to get both 'global' and 'group-level' smoothers? And is there a better way to specify the autocorrelation structure? Any suggestions about how to best proceed would be really wonderful. Thanks!

Answer 1

Autocorrelation can be seen as a trend in a time series, albeit a wiggly one. Hence you can just crank up the k if you want to get the residuals as independent as possible.

Effort should be included as an offset term using + offset(log(Effort_km2)) in the model formula. This turns the model into one for the response on the catch per unit effort scale. Making Effort_km2 a smooth term in the model is saying that the effort itself affects the response in a non-linear way. Using the offset gets a model for $\frac{\text{catch}}{\text{effort}}$, without you having to modify the response counts themselves. You can combine multiple variables in the offset term, so you can include Effort and Area by summing the log effort and log area terms.

If the fit is returned, you're just getting a singularity warning, look where the singularity os coming from. Often in these models it's the correlation term that can't be identified with a wiggly trend also in the model, for reasons due to the fact described above in my first paragraph. If you get a returned model, you can check this by looking at the estimate CAR(1) term and it's confidence interval with the intervals() function applied to the $lme component of the model fit.

You can also try simpler models and build up from there; it is likely too that the multiple Year smooths are causing you fitting problems. Instead, try a model with just one Year smooth: s(Year, by = zone, bs = "tp") so that you just get a separate smooth per zone without the global smooth.