2
$\begingroup$

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...

Thanks in advance!

$\endgroup$
2
$\begingroup$

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)

plot of measured vs fitted values

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)

plot of time smoothers

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.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.