# Incorporating auto-correlation structure into a negative binomial generalized additive mixed model using mgcv in R

I'm working with fish migration time series data, and I am modeling fish counts using environmental variables. I am using a generalized additive mixed model from the mgcv package in R. I am using a negative binomial distribution and year as a random effect. Additionally, the data has auto-correlation, and I have tried incorporating different auto-correlation structures with day nested within year, including AR1 and ARMA.

gamm(FishCount ~ s(Var1) + s(Var2) + Var3 + Var4,
family = negbin(0.82),
random = list(Year = ~ 1),
correlation = corARMA(form = ~ Day | 1, p = 1, q = 1))


Including the correlation structure gives error free output that improves the auto-correlation, substantially drops the AIC, and also improves the residual plots. However, I have not seen a correlation structure included in a negative binomial model before.

Is it valid to include AR1 or ARMA correlation structures into a negative binomial GAMM?

This is reasonable, in the same sense that working correlation matrices are included in GLMs to give the General Estimating Equations (GEE) model.

With gamm(), what you are getting is gamm() -> glmmPQL() -> lme() so you are really fitting a specially weighted linear mixed model. I won't ask how you found that $\theta = 0.82$...

I think your model specification is wrong for the correlations structure. If you want Day within Year, then the formula to use is: form = ~ Day | Year. At least, this way is explicit about the nesting.

I will add that this additive model plus correlation structure model isn't going to work in some cases, especially where the wiggliness of the smooth is of similar magnitude to the auto-correlation. Splines induce a correlation in the observations; this is most easily seen from the splines as random effects representation of the spline model. A random intercept induces a correlation between observations of the same group; similar things happen for splines.

If you have a smooth spline and a wiggly correlation (high autocorrelation) that is fine, and vice versa, but you can get into identifiability issues when trying to fit a splines and correlation structures in residuals (as here) as both terms increase in "complexity" (wiggly spline, high correlation).