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

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Tue 13 Jul 2004 - 03:18:59 EST

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