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

[Spencer Graves]

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