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)
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
