[R] ML vs. REML with gls()

[Manuel A. Morales]

Hello listmembers,

I've been thinking of using gls in the nlme package to test for serial 
correlation in my data set. I've simulated a sample data set and have 
found a large discrepancy in the results I get when using the default 
method REML vs. ML.

The data set involves a response that is measured twice a day (once for 
each level of a treatment factor). In my simulated data set, I have a 
mean "day effect" but no serial correlation between days, otherwise. 
However, if I include an AR1 structure in my model, I find significant 
evidence of serial correlation when I fit the model using maximum 
likelihood! I don't see the same result when I use restricted maximum 
likelihood. I'm not a statistician, but I'm wondering if this is 
expected. I also realize that this particular example would work well 
with day as a random effect, but then I'm not sure how I would test for 
serial correlation across days. I have Pinheiro and Bate's extremely 
useful book, but I can't seem to find anything that addresses this. I'm 
using R 1.9.1 with version 3.1-50 of nlme.

Thanks very much for any help. I've included my sample code below.

Manuel Morales

#Below I generate a day variable (note two trials per day).
day<-rep(1:25,each=2)
#Below I generate a random "day effect"
day.effect<-c(); for(i in 1:25) {
	tmp<-rnorm(1);
	day.effect<-append(rep(tmp,2),day.effect)}
#Below I randomize the application of the treatment effect w/i days.
trt<-c(); for(i in 1:25) {
	if (rnorm(1)>0) trt<-append(c(0,1),trt)
	else trt<-append(c(1,0),trt)}
#Below I generate the response variable. Trt increases the response.
resp<-trt+day.effect+rnorm(50)

#Below I analyze the data using REML or ML w/ and w/o AR1 errors.
library(nlme)
base.gls<-gls(resp~factor(trt)+factor(day))
ar.gls<-gls(resp~factor(trt)+factor(day),corr=corAR1())
base.gls.ml<-gls(resp~factor(trt)+factor(day),method="ML")
ar.gls.ml<-gls(resp~factor(trt)+factor(day),corr=corAR1(),method="ML")

#Below I compare models using REML or ML
anova(base.gls,ar.gls); anova(base.gls.ml,ar.gls.ml)