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!