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

[Jacob Wegelin]

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 at ucdavis.edu
> >
> >______________________________________________
> >R-help at 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
> >
> >
>