R help - Jul 2004 - lme unequal random-effects variances varIdent pdMat Pinheiro Bates nlme

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

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
>  
>

How do I read a Systat file into R?  The  following email by Marc Schwartz
http://tolstoy.newcastle.edu.au/R/help/04/06/1005.html deals with reading Systat
on
a Linux machine.  I'm running R on Windows (precise version info below).  A
colleague
sent me data in systat, a *.SYD file.  Do I have to ask the colleague to re-save
the
data as Excel or text?

(I tried opening it in SPSS, and SPSS also does not seem to recognize an SYD
file.)

Thanks for any info.

Jake Wegelin
> version
         _
platform i386-pc-mingw32
arch     i386
os       mingw32
system   i386, mingw32
status
major    1
minor    9.1
year     2004
month    06
day      21
language R