From: Jacob Wegelin <jawegelin_at_ucdavis.edu>

Date: Mon 12 Jul 2004 - 17:50:58 EST

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 Mon Jul 12 17:58:30 2004

Date: Mon 12 Jul 2004 - 17:50:58 EST

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 Received on Mon Jul 12 17:58:30 2004

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