Re: [R] lme unequal random-effects variances varIdent pdMat Pinheiro Bates nlme

From: Jacob Wegelin <jawegelin_at_ucdavis.edu>
Date: Tue 13 Jul 2004 - 03:41:07 EST

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/
> >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 04:06:43 2004

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