Dear Mr. Graves: Thank you. You're right -- I should have said method="ML".
Putting sex into the model as a fixed effect enables us to see if the relationship between time and y differs between the sexes. But it doesn't do anything to the covariance of the random effects. Jake On Mon, 12 Jul 2004, Spencer Graves wrote: > Have you tried something like the following: > > mylmeMF0 <- lme( y ~ time, data=DAT, random=~time | id, method="ML") > > > > mylmeMF1 <- lme( y ~ time*sex, data=DAT, random=~time | id, method="ML") > > > anova(mylemFM0, mylmeMF1) > > Please check Pinheiro and Bates regarding "method". They explain > when "ML" and "REML" are appropriate. I don't have the book handy, and > I can't remember for sure, but I believe that "REML" is would only be > apropriate here if "sex" were NOT in the "fixed" model. > > hope this helps. spencer graves > > > Jacob Wegelin wrote: > > >How does one implement a likelihood-ratio test, to test whether the > >variances of the random effects differ between two groups of subjects? > > > >Suppose your data consist of repeated measures on subjects belonging to > >two groups, say boys and girls, and you are fitting a linear mixed-effects > >model for the response as a function of time. The within-subject errors > >(residuals) have the same variance in both groups. But the dispersion of > >the random effects differs between the groups. The boys' random effects > >-- say, the intercepts -- have greater variance than the girls'. One can > >see this by partitioning the data by sex and fitting two separate models. > > > >The model for the girls, > > > > library("nlme") > > mylmeF0 <- lme( y ~ time, data=DAT, random=~time | id, subset=sex=="F") > > > >yields a variance of about one for the random intercepts: > > > > StdDev Corr > >(Intercept) 0.9765052 (Intr) > >time 0.1121913 -0.254 > >Residual 0.1806528 > > > >whereas in the model for the boys, the corresponding variance is ten times > >that amount: > > > > mylmeM0 <- lme( y ~ time, data=DAT, random=~time | id, subset=sex=="M") > > > > StdDev Corr > >(Intercept) 10.1537946 (Intr) > >time 0.1230063 -0.744 > >Residual 0.1298910 > > > >I would like to use a likelihood ratio to test this difference. The > >smaller ("null") model would be > > > > mylme0 <- lme( y ~ time, data=DAT, random=~time | id ) . > > > >This model forces the random intercepts for both boys and girls to come > >from a single normal distribution. > > > >The larger model would allow the boys' and girls' random intercepts (or > >more generally their random effects) to come from separate normal > >distributions with possibly unequal variances. > > > >There must be some straightforward obvious way to fit the larger model, > >but I do not see it. > > > >Pinheiro and Bates, chapter 5.2, show how to model unequal *residual* > >("within-group") variances for the two groups using varIdent. They also > >tantalizingly say, "The single-level variance function model (5.10) can be > >generalized to multilevel models" (page 206), which seems to suggest that > >a solution to the current problem might exist. > > > >The pdMat classes provide a way to *constrain* the dispersion matrix of > >the random effects, not make it more general. > > > >Of course, one way to test for unequal variances is to apply an F-test for > >equal variances to the random intercepts. If the data are as shown at the > >bottom of this email, the test can be implemented as follows: > > > > stuff<-as.data.frame(summary(mylme0)$coefficients$random$id) > > stuff$sex<-factor(substring(row.names(stuff), 1,1)) > > mysplit<-split(stuff[,"(Intercept)"], stuff[,"sex"]) > > ns<-sapply(mysplit, length) > > vars<-sapply(mysplit, var) > > p<- 1-pf( vars["M"]/vars["F"], ns["M"]-1, ns["F"]-1) > > > >Alternatively, one could implement a permutation test for the ratio of the > >variances of the random intercepts--these variances derived from the two > >halves of the partitioned data. > > > >But surely there's a direct, model-based way to do this? > > > >Thanks for any suggestions > > > >Jake > > > >P.S. Here is the code by which the "data" were generated. > > > > nb<-1 > > ntimepts<-3 > > girls<-data.frame( > > y= rep(-nb:nb , each=ntimepts) > > , > > id=rep( paste("F", 1:(2*nb+1), sep=""), each=ntimepts) > > , > > time=rep(1:(2*nb+1), length=ntimepts) > > ) > > boys <-data.frame( > > y= rep(10*(-nb:nb) , each=ntimepts) > > , > > id=rep( paste("M", 1:(2*nb+1), sep=""), each=ntimepts) > > , > > time=rep(1:(2*nb+1), length=ntimepts) > > ) > > DAT<-rbind(girls,boys) > > DAT$y<-DAT$y + rnorm(nrow(DAT))/5 > > DAT$sex<-factor(substring( as.character(DAT[,"id"]), 1,1)) > > row.names(DAT)<-paste( DAT[,"id"], DAT[,"time"], sep=".") > > > >Jacob A. Wegelin > >Assistant Professor > >Division of Biostatistics, School of Medicine > >University of California, Davis > >One Shields Ave, TB-168 > >Davis CA 95616-8638 USA > >http://wegelin.ucdavis.edu/ > >[EMAIL PROTECTED] > > > >______________________________________________ > >[EMAIL PROTECTED] mailing list > >https://www.stat.math.ethz.ch/mailman/listinfo/r-help > >PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html > > > > > ______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html