I am seeking advice on how to effectively eliminate autocorrelation from a linear mixed model. My experimental design and explanation of fixed and random factors can be found here from an earlier question I asked:
Crossed fixed effects model specification including nesting and repeated measures using glmm in R
I have treated day as numeric even though I only have four sampling time points (so I could treat it as a categorical predictor as well). Aside: Although four sample points is very few, I don’t think that this is the root of the problem as this same dataset is giving me this residual autocorrelation issues using a different response variable that has 24 time points.
My issue is that I have tried a number of different autocorrelation structures and can’t seem to achieve the random, non-significant residuals needed to confirm a lack of autocorrelation. I am using the function lme in the R package nlme to deal with autocorrelation.
I have tried the various autocorrelation classes with variations to form
1) corAR1 (autoregressive process of order 1).
2) corARMA (autoregressive moving average process)
3) corCAR1 (continuous autoregressive process)
4) corGaus (Gaussian spatial correlation)
With form varying in the following ways with these different autocorrelation classes:
form=~1
form=~1| TankNumb/RecruitID2
form=~Day| TankNumb/RecruitID2
If we look at a model without the time factor "Day" added, the ACF and PACF plots look like this.
lme4_lognormal_notime<-lmer(Arealog~Temperature*Culture+(1|TankNumb/RecruitID2), data=growthSR_noNA)
acf(residuals(lme4_lognormal_notime, retype="normalized"))
pacf(residuals(lme4_lognormal_notime, retype="normalized"))
Also, if I look at the residuals of the model without “Day” included, I do not see any strong pattern in the residuals that would make me think there is a temporal autocorrelation problem.
plot(residuals(lme4_lognormal_Ben_notime, retype="normalized")~growthSR_noNA$Day)
Now for two different models with autocorrelation structure to hopefully eliminate autocorrelation:
nlme_lognormal_mult_cor<-lme(Arealog~Temperature*Culture*Day, random=~1|TankNumb/RecruitID2,correlation=corAR1(form=~1), data=growthSR_noNA)
nlme_lognormal_mult_cortime<-lme(Arealog~Temperature*Culture*Day, random=~1|TankNumb/RecruitID2,correlation=corAR1(form=~Day|TankNumb/RecruitID2), data=growthSR_noNA)
ARMA_nlme_lognormal_mult_cor<-lme(Arealog~Temperature*Culture*Day, random=~1|TankNumb/RecruitID2,correlation=corARMA(form=~1, p=0, q=1), data=growthSR_noNA)
The AIC suggests that the simplest correlation structure is the best.
AIC(nlme_lognormal_mult,nlme_lognormal_mult_cor, nlme_lognormal_mult_cortime,ARMA_nlme_lognormal_Ben_mult_cor)
df AIC
nlme_lognormal_mult 15 1233.997
nlme_lognormal_mult_cor 16 1184.389
nlme_lognormal_mult_cortime 16 1235.997
ARMA_nlme_lognormal_Ben_mult_cor 16 1198.451
As I mentioned above, I have tried a number of different cor functions (the four listed above) and different form specifications. They all end up with ACF/PCF plots like the last two models with a first lag at below 0.2 in the ACF plot and a PCF plot with the first three lags around 0.10.
I have also read a number of sites describing how to specify corARMA models based on diagnosing the ACF plots and have tried a number of variations of p and q parameters.
Questions:
- Does anyone have some advice on which type of correlation structure that might elimate this autocorrelation problem based on the patterns in my ACF/PCF plots? Should I be diagnosing based on a model with or without Day included?
2.Is there ever an acceptable level of autocorrelation? This post (Do autocorrelated residual patterns remain even in models with appropriate correlation structures, & how to select the best models?) states that small amounts of autocorrelation probably won't impact the model coefficients very much. "The estimate is slightly larger than zero so will have negligible effect on the model fit and hence you might wish to leave it in the model if there is a strong a priori reason to assume residual autocorrelation." Potentially there is some autocorrelation that is not being caused by temporal autocorrelation, like outliers? Is there a cut-off, for example, autocorrelation below 0.1? I have extremely small 95% confidence intervals, so it doesn't take a lot of autocorrelation in my models to be significantly too much.
Any advice is appreciated!
