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

Question

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.

Answer 1

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 
#Link function: identity 
#
#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.