This question uses R, but is about general statistical issues.
I'm analysing the effects of mortality factors (% mortality due to disease and parasitism) on moth population growth rate over time, where larval populations were sampled from 12 sites once a year for 8 years. The population growth rate data displays a clear but irregular cyclical trend over time.
The residuals from a simple generalised linear model (growth rate ~ %disease + %parasitism + year) displayed a similarly clear but irregular cyclical trend over time. Therefore, generalised least squares models of the same form were also fitted to the data with appropriate correlation structures to deal with the temporal autocorrelation, e.g. compound symmetry, autoregressive process order 1 and autoregressive moving average correlation structures.
Models all contained the same fixed effects, were compared using AIC, and were fitted by REML (to allow comparison of different correlation structures by AIC). I'm using the R package nlme and the gls function.
The GLS models' residuals still display almost identical cyclical patterns when plotted against time. Will such patterns always remain, even in models that accurately account for the autocorrelation structure?
I have simulated some simplified but similar data in R below my second question, which shows the issue based on my current understanding of the methods needed to assess temporally autocorrelated patterns in model residuals, which I now know are wrong (see answer).
I have fitted GLS models with all possible plausible correlation structures to my data, but none are actually substantially better fitting than the GLM without any correlation structure: just one GLS model is marginally better (AIC score = 1.8 lower), whilst all the rest have higher AIC values. However, this is only the case when all models are fitted by REML, not ML where GLS models are clearly much better, but I understand from stats books you must only use REML to compare models with different correlation structures and the same fixed effects for reasons I won't detail here.
Given the clearly temporally autocorrelated nature of the data, if no model is even moderately better than the simple GLM what is the most appropriate way to decide which model to use for inference, assuming I'm using an appropriate method (I eventually want to use AIC to compare different variable combinations)?
Q1 'simulation' exploring residual patterns in models with and without appropriate correlation structures
Generate simulated response variable with a cyclical effect of 'time', and a positive linear effect of 'x':
time <- 1:50 x <- sample(rep(1:25,each=2),50) y <- rnorm(50,5,5) + (5 + 15*sin(2*pi*time/25)) + (x/1)
y should display a cyclical trend over 'time' with random variation:
And a positive linear relationship with 'x' with random variation:
Create a simple linear additive model of "y ~ time + x":
require(nlme) m1 <- gls(y ~ time + x, method="REML")
The model displays clear cyclical patterns in the residuals when plotted against 'time', as would be expected:
And what should be a nice, clear lack of any pattern or trend in the residuals when plotted against 'x':
A simple model of "y ~ time + x" that includes an autoregressive correlation structure of order 1 should fit the data much better than the previous model because of the autocorrelation structure, when assessed using AIC:
m2 <- gls(y ~ time + x, correlation = corAR1(form=~time), method="REML") AIC(m1,m2)
However, the model should still display almost identically 'temporally' autocorrelated residuals:
Thank you very much for any advice.