# LMM model with nested design and random slopes: how to fit the autocorrelation term?

I am currently running a linear mixed model with a nested design and random slopes.

For example, let's imagine some monthly captures of wild rabbits in kilograms in 5 sites during 21 years:

site<- rep(rep(c("Golden Cave","Ringo's place","Damned Dam","Knockampton","Easy Fuzzy"),each=12),21)
year <- rep(2000:2020, each=12*5)
month <- rep(seq(1,12),21*5)
rabbit_captures <- rnorm(12*21*5, 50, 10)
dataset <- as.data.frame(cbind(site,year,month,rabbit_captures))
dataset$$rabbit_captures <- as.numeric(dataset$$rabbit_captures)


And the corresponding model:

library(nlme)
library(MASS)

model_lme <- lme(fixed = log(rabbit_captures) ~ site,
random = ~ site|year/month,
data = dataset, method = "ML",
control = lmeControl(opt = 'optim'))


With this simulated distribution of rabbit_captures, the model does not converge, sorry I tried, but it still corresponds to the design of my experiment. In other words, with a more complex dataset, this model runs properly. However, I am uncomfortable doing spatial modeling without considering time autocorrelation...

So here is my question:

How to consider time autocorrelation within the model?

Knowing that:

• correlation = corAR1(value = 0.9, form = ~ site|year/month) is an incorrect formulation. Variable "Site" is made of factors and "A covariate for this correlation structure must be integer valued."

• correlation = corAR1(value = 0.9, form = ~ 1|year/month) should work (with my example, it won't) but I really don't know if it makes sense statistically...

Typically, you'd use something like corCAR1(form = ~ date | site) or corAR1(form = ~ months_since_start_of_timeseries | site). This specifies time as an auto-correlation co-variate and groups by site. If your time series are perfectly regular and the data sorted by time, you could use corAR1(form = ~ 1 | site).

However, I suggest that you do not use a mixed model and include time as a fixed effect. Since this would not be a linear effect, you could use a GAM.

#never use as.data.frame(cbind(...)), do this instead:
dataset <- data.frame(site,year,month,rabbit_captures)

library(mgcv)
dataset$$site <- factor(dataset$$site)
dataset$$yearmon <- (dataset$$year - min(dataset$$year)) * 12 + dataset$$month

fit <- gam(log(rabbit_captures) ~  site + s(yearmon, by = site), data = dataset)
gam.check(fit)
abline(0, 1)


summary(fit)
#Family: gaussian
#
#Formula:
#log(rabbit_captures) ~ site + s(yearmon, by = site)
#
#Parametric coefficients:
#                   Estimate Std. Error t value Pr(>|t|)
#(Intercept)        3.886624   0.013556 286.714   <2e-16 ***
#siteEasy Fuzzy    -0.005637   0.019171  -0.294    0.769
#siteGolden Cave    0.008913   0.019171   0.465    0.642
#siteKnockampton   -0.023536   0.019171  -1.228    0.220
#siteRingo's place  0.001926   0.019171   0.100    0.920
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Approximate significance of smooth terms:
#                               edf Ref.df     F p-value
#s(yearmon):siteDamned Dam    5.826  6.984 1.116   0.338
#s(yearmon):siteEasy Fuzzy    1.000  1.000 2.710   0.100 .
#s(yearmon):siteGolden Cave   2.114  2.638 1.390   0.316
#s(yearmon):siteKnockampton   5.523  6.667 1.597   0.138
#s(yearmon):siteRingo's place 1.000  1.000 0.814   0.367
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#R-sq.(adj) =  0.0132   Deviance explained = 2.84%
#GCV = 0.047072  Scale est. = 0.046307  n = 1260

plot(fit, pages = 1)


The time smoothers should take care of auto-correlation. You could also use an additional time smoother independent of site to model a general, site-independent trend.