if (FALSE) {
library(agridat)
data(pederson.lettuce.repeated)
dat <- pederson.lettuce.repeated
libs(lattice)
dat <- dat[order(dat$day),]
xyplot(weight ~ day|trt, dat, type='l', group=plant, layout=c(3,1),
main="pederson.lettuce.repeated")
# Pederson used this SAS MIXED model for unstructured covariance
# proc mixed data=Project.Spacingdata;
# class trt plant day;
# model weight=trt day trt*day;
# repeated day / subject=plant type=un r rcorr;
# This should give the same results as SAS, but does not.
libs(nlme)
dat <- transform(dat, plant=factor(plant), day=factor(day))
datg <- groupedData(weight ~ day|plant, data=dat)
un1 <- gls(weight ~ trt * day, data=datg,
correlation=corSymm(value=rep(.6,55), form = ~ 1 | plant),
control=lmeControl(opt="optim", msVerbose=TRUE,
maxIter=500, msMaxIter=500))
logLik(un1)*2 # nlme has 1955, SAS had 1898.6
# Comparing the SAS results in Pederson (page 16) and the nlme results, we notice
# the SAS correlations in table 5.2 are unusually low for the first
# column. The nlme results have a higher correlation in the first column
# and just "look" better
un1
}
Run the code above in your browser using DataLab