[R] ML vs. REML with gls()

Manuel A. Morales Manuel.A.Morales at williams.edu
Fri Sep 3 23:02:01 CEST 2004

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).
#Below I generate a random "day effect"
day.effect<-c(); for(i in 1:25) {
#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.

#Below I analyze the data using REML or ML w/ and w/o AR1 errors.

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

More information about the R-help mailing list