Testing if lme with corAR1() fits significantly better than without this correction for first-order autocorrelation

[Patrick_P9]

Dear everyone

I want to test if the correction for first order autocorrelation corAR1() improves significantly my model fit but don't know how to do this.

Accoring to Pinheiro & Bates (2000) p. 242 this should be done using the standard anova function to perform a Likelihood Ration Test. However I was told that from a theoretical point of view this would render inexact results because the chi-square-distribution-assumption is violated. Still I could not get an input on how to consider this in R.

Therefore I wanted to ask if anyone could lend me a hand on which alternatives could be used instead of the anova-LRT? It would be very appreciated.

dataset and if prefered r-script for the two models with and without corAR1() can be retrieved here

##reduced model-specific data-set:
datafclr <-read.csv("datafclr.csv", header = TRUE, sep = ",", dec = ".", fill = TRUE)

##without correction for autocorrelation:

tim1 <- lme(fixed=EERTmn ~ male + female + 
              (male:time7c)          + (female:time7c)          +                           
              (male:IERT_Cp)         + (female:IERT_Cp)         + 
              (male:IERT_Cp_Partner) + (female:IERT_Cp_Partner)-1,

            control=list(maxIter=100000), data=datafclr,                                             

            random=~male + female -1|dyade, na.action=na.omit) 

summary(tim1)

#with correction for autocorrelation of first order:

tim2 <- lme(fixed=EERTmn ~ male + female + 
              (male:time7c)          + (female:time7c)          +                           
              (male:IERT_Cp)         + (female:IERT_Cp)         + 
              (male:IERT_Cp_Partner) + (female:IERT_Cp_Partner)-1,

            control=list(maxIter=100000), data=datafclr,                                             

            random=~male + female -1|dyade, correlation=corAR1(), na.action=na.omit)

summary(tim2)

Many thanks in advance!

Best, Patrick

[alexpghayes]

This preprint includes a brief discussion of this exact problem.

The short version is you can use (ripping an example pretty much directly from ?corCAR1()):

library(nlme)
data("Ovary")

no_car <- lme(
  follicles ~ sin(2*pi*Time) + cos(2*pi*Time),
  data = Ovary,
  random = pdDiag(~sin(2*pi*Time))
)

with_car <- update(no_car, correlation = corCAR1(form = ~Time))

anova(no_car, with_car)
#>          Model df      AIC      BIC    logLik   Test  L.Ratio p-value
#> no_car       1  6 1638.082 1660.404 -813.0409                        
#> with_car     2  7 1565.535 1591.577 -775.7676 1 vs 2 74.54649  <.0001

^{Created on 2018-12-11 by the reprex package (v0.2.1)}

You can also use mgcv::gamm(). Some examples of that included in the paper above.