From: Spencer Graves <spencer.graves_at_pdf.com>

Date: Tue 13 Jul 2004 - 03:18:59 EST

R-help@stat.math.ethz.ch 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 Received on Tue Jul 13 03:28:09 2004

Date: Tue 13 Jul 2004 - 03:18:59 EST

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 explainwhen "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/
**>jawegelin@ucdavis.edu
**>
**>______________________________________________
**>R-help@stat.math.ethz.ch 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
**>
**>
*

R-help@stat.math.ethz.ch 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 Received on Tue Jul 13 03:28:09 2004

*
This archive was generated by hypermail 2.1.8
: Fri 18 Mar 2005 - 02:35:38 EST
*