[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).
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)
More information about the R-help
mailing list