2
$\begingroup$

I analyzed the effect of temperature (4 different areas) on laying date: LDT ~ Aa3+Bb+Cc+Dd. Because of autocorrelation in residuals I used gls models with method="RELM" to find the correlation structure of a model set in which I included different combinations of fixed effects. I used acf plots to see if the issue was gone. I also tried Breusch-Godfrey and Durbin-Watson tests, but I saw in a post that these would not work well when autocorrelation is at high lags like in my case. Thus according to acf plots only some models needed a correlation structure. For model selection comparing nested models I used the method="ML". Please help me to understand this, it's getting more and more confusing. What's the right/suggested criteria for model selection when I need a correlation structure and test different effects?

  1. Should I use the correlation structure found with REML in ML models? Or do I have to check again for autocorrelation in ML models? (which I guess doesn't make sense because that's why we use REML). Am I right only checking for other model assumptions ignoring residuals autocorrelation? In two cases the autocorrelation structure from REML used with ML shows autocorrelation in the acf plot. Is it correct to use it?

  2. Does residual autocorrelation at high lag need a correlation structure? Visually there's not much signs of heterogeneity and in a few cases I wasn't able to find a correlation structure (models m8 and m11 method=REML).

Here's my data:

structure(list(year = 1991:2015, Aa3 = c(11.7967532509252, 11.4555638363487, 
13.1706280624658, 9.83853834721082, 11.0892774079976, 9.28558456688599, 
8.37008645576346, 16.685201223273, 14.9208415784334, 12.5005737304688, 
14.9006742778577, 9.53653243215462, 10.6907120554071, 13.3834204958197, 
9.8151839406867, 13.2178794459293, 11.6116308379592, 11.5517068427906, 
12.2147025125069, 11.5025536252741, 11.4172066671807, 12.803326094778, 
15.5187240868284, 12.148594585218, 15.0073361045436), Bb = c(-3.18143519666869, 
-4.39866028668509, 1.78033518645253, -1.09087811744935, 0.781504864298432, 
-3.35258413204158, -0.147817393686718, 2.20858102713945, 2.78942101113926, 
0.133688841044735, 7.54352335295978, -3.01163687673316, 4.76122905649298, 
-0.067047827360323, 6.75503612812823, 4.59545579459801, -3.10992282189532, 
5.19488285181458, -1.75941821116188, 3.67164104630711, -1.37833278441656, 
0.889547512007415, -1.47006917494296, 2.96885387269044, -0.599857256084644
), Cc = c(-3.14888958162661, -5.09421300863669, 0.859087760125301, 
-3.77697950592429, -2.35219581389364, -3.95249382467829, -3.33432632797512, 
-1.20699309912482, -2.85922008621721, -3.713746771849, -1.33312677641961, 
-2.46973133916438, -1.72655447645258, -3.11445919729557, -1.25910057594707, 
-1.57561781314935, -4.75200348290639, 0.160156234390008, -0.626073828499617, 
-1.68243433179146, -3.08619266129513, -2.68649754048611, -4.31190628871281, 
-1.472913440597, -3.08940582665639), Dd = c(-1.37625014361211, 
-5.05573658662682, -1.18101124482994, -2.52954173368564, -0.261439424402551, 
-3.97714305204502, -1.98631268669575, 1.59101167566638, -3.05814029469207, 
-5.33150060317094, -2.69995009478398, -2.56073069852939, -1.92038394703582, 
-2.70067713120402, -0.0833998736213011, -3.31705501780788, -3.06175213982075, 
-1.19264921300549, -2.02601964613968, -2.31430484547332, -1.1192687988281, 
-1.61056051815255, -1.08186142865347, -2.30663057215071, -2.99631168141082
), medini = c(162L, 172L, 157L, 161L, 161L, 166L, 161L, 158L, 
168L, 168L, 164L, 167L, 160L, 163L, 163L, 165L, 167L, 162L, 164L, 
164L, 163L, 164L, 165L, 162L, 163L)), .Names = c("year", "Aa3", 
"Bb", "Cc", "Dd", "medini"), class = "data.frame", row.names = c(NA, 
-25L))

Models 0 - 7, 9, 12 didn't need a correlation structure. But as I mentioned, there were two models I wasnt able to find one.

require(nlme)
m8 = gls(medini~Bb+Cc ,method="REML", control = list(maxIter=100),data=LDT )
plot(ACF(m8,resType="normalized"),alpha=0.05)
m8.1 <- update(m8, correlation=corARMA(form = ~ 1,p=1))
plot(ACF(m8.1,resType="normalized"),alpha=0.05)
m8.1.1 <- update(m8, correlation=corARMA(form = ~ 1,p=1,q=1))
plot(ACF(m8.1.1,resType="normalized"),alpha=0.05)
m8.2 <- update(m8, correlation=corARMA(form = ~ 1,p=2))
plot(ACF(m8.2,resType="normalized"),alpha=0.05)
m8.2.1 <- update(m8, correlation=corARMA(form = ~ 1,p=2, q=1))
plot(ACF(m8.2.1,resType="normalized"),alpha=0.05)
m8.2.2 <- update(m8, correlation=corARMA(form = ~ 1,p=2, q=2),control=glsControl(tolerance=1e-4,msVerbose=F))
plot(ACF(m8.2.2,resType="normalized"),alpha=0.05)
m8.3 <- update(m8, correlation=corARMA(form = ~ 1,p=3))

anova(m10,m10.1,m10.1.1,m10.2,m10.2.1,m10.2.2)

I show details for models m13 and m15 because these correlation structures are different when I use ML method:

m13 = gls(medini~Bb+Cc+Dd  ,method="REML", data=LDT)
plot(ACF(m13,resType="normalized"),alpha=0.05)
m13.1 <- update(m13, correlation=corARMA(form = ~ 1,p=1))
plot(ACF(m13.1,resType="normalized"),alpha=0.05)
m13.1.1 <- update(m13, correlation=corARMA(form = ~ 1,p=1,q=1))
plot(ACF(m13.1.1,resType="normalized"),alpha=0.05)
m13.2 <- update(m13, correlation=corARMA(form = ~ 1,p=2))
plot(ACF(m13.2,resType="normalized"),alpha=0.05)
m13.2.1 <- update(m13, correlation=corARMA(form = ~ 1,p=2,q=1),control=lmeControl(opt="optim"))
plot(ACF(m13.2.1,resType="normalized"),alpha=0.05)
anova(m13,m13.1,m13.1.1,m13.2,m13.2.1)

m14.2.1 <- gls(medini~Aa3+Cc+Dd  ,method="REML", data=LDT, correlation=corARMA(form = ~ 1,p=2,q=1))
plot(ACF(m14.2.1,resType="normalized"),alpha=0.05)

m15 = gls(medini~Aa3+Bb+Cc+Dd ,method="REML", data=LDT)
plot(ACF(m15,resType="normalized"),alpha=0.05)
m15.1 <- update(m15, correlation=corARMA(form = ~ 1,p=1))
plot(ACF(m15.1,resType="normalized"),alpha=0.05)
m15.1.1 <- update(m15, correlation=corARMA(form = ~ 1,p=1,q=1))
plot(ACF(m15.1.1,resType="normalized"),alpha=0.05)
m15.2 <- update(m15, correlation=corARMA(form = ~ 1,p=2))
plot(ACF(m15.2,resType="normalized"),alpha=0.05)
m15.2.1 <- update(m15, correlation=corARMA(form = ~ 1,p=2,q=1))
plot(ACF(m15.2.1,resType="normalized"),alpha=0.05)
anova(m15,m15.1,m15.1.1,m15.2,m15.2.1)

Now these are all the models with ML to compare fixed effects in nested models using selected corARMA from RELM

m0 = gls(medini~1,method="ML", data=LDT)
m1 = gls(medini~Aa3,method="ML", data=LDT)
m2 = gls(medini~Bb,method="ML", data=LDT)
m3 = gls(medini~Cc,method="ML", data=LDT)
m4 = gls(medini~Dd,method="ML", data=LDT)
m5 = gls(medini~Aa3 + Bb ,method="ML", data=LDT)
m6 = gls(medini~Aa3 + Cc ,method="ML", data=LDT )
m7 = gls(medini~Aa3 + Dd ,method="ML",data=LDT )
m8 = gls(medini~Bb+Cc ,method="ML", control = list(maxIter=100),data=LDT )
m9 = gls(medini~Bb+Dd ,method="ML",data=LDT  )
m10.2.1 = gls(medini~Cc+Dd  ,method="ML",data=LDT,correlation=corARMA(form = ~ 1,p=2,q=1))
m11 = gls(medini~Aa3+Bb+Cc  ,method="ML",data=LDT)
m12 = gls(medini~Aa3+Bb+Dd  ,method="ML",data=LDT )

m13.2.1 and m15.2.1 show autocorrelation at high lag, should I ignore this because it was already defined with REML? or should I take m13.3 and m15.2.2 instead?

m13 = gls(medini~Bb+Cc+Dd  ,method="ML", data=LDT)
m13.2.1 <- update(m13, correlation=corARMA(form = ~ 1,p=2,q=1))
plot(ACF(m13.2.1,resType="normalized"),alpha=0.05)
m13.2.2 <- update(m13, correlation=corARMA(form = ~ 1,p=2,q=2),control=lmeControl(opt="optim"))
plot(ACF(m13.2.2,resType="normalized"),alpha=0.05)
m13.3 <- update(m13, correlation=corARMA(form = ~ 1,p=3))
plot(ACF(m13.3,resType="normalized"),alpha=0.05)

m14.2.1 <- gls(medini~Aa3+Cc+Dd,method="ML", data=LDT, correlation=corARMA(form = ~ 1,p=2,q=1),control=lmeControl(opt="optim"))

m15.2.1 <- gls(medini~Aa3+Bb+Cc+Dd ,method="ML", data=LDT, correlation=corARMA(form = ~ 1,p=2,q=1))
plot(ACF(m15.2.1,resType="normalized"),alpha=0.05)
m15.2.2 <- update(m15, correlation=corARMA(form = ~ 1,p=2,q=2),control=lmeControl(opt="optim"))
plot(ACF(m15.2.2,resType="normalized"),alpha=0.05)

anova(m0,mAa3,mBb,mCc,mDd,m5,m6,m7, m8,m9,m10.2.1,m11,m12,m13.2.2,m14.2.1, m15.2.2)
$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.